Rotating shallow water (RSW) model is a classical modeling tool in geophysical fluid dynamics, which helps to understand a variety of major dynamical phenomena in the atmosphere and oceans by simple analytical and computational means; see, e.g., . It has, however, a soft spot: the absence of horizontal gradients of temperature and/or density. The RSW model can be obtained from the full “primitive” hydrostatic equations of the geophysical fluid dynamics by vertical averaging under the hypothesis of columnar motion. As it is quite common in geophysical modeling, we will be considering the fluid on the tangent to the rotating planet plane, and take into account only the normal to the plane component of the angular velocity of planet’s rotation (so-called “traditional” approximation). The assumption of horizontal homogeneity of temperature/density accompanies the standard derivation. Yet, this assumption can be relaxed, which does not substantially alter the derivation and leads to the so-called thermal shallow water (TSW) model , which we will call thermal rotating shallow water (TRSW) equations in the presence of rotation. This model was repeatedly rediscovered and used in the literature both in the atmospheric and oceanic context, in particular for studies of the mixed layer in early papers [16, 32, 35, 38, 39, 48]. It was recently applied to planetary atmospheres [14, 45]. Structurally, while the classical shallow water equations are equivalent to those of isentropic gas dynamics with pressure depending only on density, the TSW model corresponds to the dynamics of a gas with a specific equation of state, depending both on density and temperature; see .
One of the principal applications of the TRSW equation is to model atmospheric and oceanic temperature fronts; see [16, 47]. This is why a numerical method capable of accurately resolving sharp temperature and pressure fronts is desirable. Development of such method will be in the focus of this paper.
We consider a one-dimensional (1-D) TRSW equations with nonflat bottom topography, Coriolis force, buoyancy, and without dependence on the zonal coordinate . The studied system (see [21, 22]) reads as
where is a meridional coordinate, is time, is the Lagrangian (material) derivative, denotes the thickness of the fluid layer, represents the bottom topography, and stand for the zonal and meridional velocities, respectively, and will be called buoyancy. In the oceanographic applications, , where is the density difference between the density and a reference density . Notice that density variations in the ocean are proportional to temperature variations, whence the name of the model. In the atmospheric applications, and should be replaced with the potential temperature and , while . In both cases, is the acceleration due to gravity. Finally, rotation enters via , the Coriolis parameter. In the simplest -plane approximation, where the effects of curvature are neglected, is assumed to be constant, that is, . If the curvature is taken into account to the first order in the beta-plane approximation, becomes , which reduces to near the Equator, where the tangent plane is parallel to the axis of rotation.
As one can easily show, the system (1.1) can be rewritten in the form of the following system of balance laws:
where and are the zonal and meridional discharges or, in other words, the densities of zonal and meridional momenta, respectively.
We notice that if the Coriolis force is not considered () and so that , then the system (1.2) is reduced to the classical Saint-Venant system of shallow water equations, which was originally derived in  and is still widely used in modeling water flows in rivers, canals, lakes, reservoirs and coastal areas. The 1-D Saint-Venant system reads as:
It is well-known that steady-state solutions are of great particular importance because they minimize, if they are stable, the energy, and many practically relevant waves can be viewed as small perturbations of those equilibria. The system, whenever there is an energy sink of any nature, tends to reach a steady state. This process is called adjustment. The steady-state solutions of the Saint-Venant system (1.3) satisfy the time-independent system:
It is easy to show that the system (1.4) admits a family of smooth steady-state solutions:
Among which, one of the most practical relevant steady states is the so-called “lake at rest” equilibria, which is obtained from (1.5) when :
Capturing these steady states or their small perturbations (quasi-steady flows) accurately is a highly nontrivial task, because it may lead to spurious oscillations when the shock capturing numerical methods are directly applied. Although such unphysical oscillations may be attenuated once a very fine computational grid is used, this may not be affordable in practice where coarse grids have to be used. Therefore, in the past decades, there have been many attempts to develop well-balanced numerical methods that are capable of exactly preserving the aforementioned steady-states solutions at the discrete level; see, e.g., [1, 2, 3, 7, 8, 9, 31, 25, 29, 33, 46] and references therein.
Similar to the Saint-Venant system (1.3), there are some special steady states of the 1-D TRSW system (1.2). However, due to the presence of Coriolis force and density/temperature variations, the structure of the steady state solutions becomes more complex. In fact, the steady-state solutions of (1.2) can be obtained by solving the system:
As one can easily see, in the absence of Coriolis force , the system admits several particular steady-state solutions, two of which are the following “lake at rest” ones:
A well-balanced central-upwind scheme supplemented with interface tracking method for preserving the steady states (1.6) and (1.7) was introduced in . A path-conservative approximate Riemann-problem-solver based scheme for the TSW system (1.2) with was developed in . This scheme is well-balanced, positivity preserving and entropy dissipative in the case of flat or continuous bottom topography . For several other recently developed well-balanced schemes for the TSW model, we refer the reader to [17, 23, 24, 42].
If both the Coriolis force and buoyancy are taken into account, the situation is even more complicated and it is quite challenging to design an accurate and robust well-balanced scheme for the studied system. In this paper, we follow the idea from [8, 10, 12], and incorporate the source term in (1.2c) into the corresponding flux term and rewrite (1.2) in the following equivalent form:
is a global equilibrium variable and
It is easy to see that the system (1.8) admits the equilibria, which will be called thermo-geostrophic , as they are due both to effects of rotation and temperature variation, and can be written using the global variables in an extremely simple form:
On the other hand, (1.8) is a hyperbolic system with a global flux and therefore it might be quite challenging to design a Riemann-problem-solver-based upwind scheme for solving it numerically.
In this paper, our main goal is to develop a robust, stable and highly accurate finite-volume method for the TRSW system (1.8). We note that in this system, the flux is global as it depends on a global variable . This makes it hard to design upwind methods for (1.8). We therefore derive a Riemann-problem-solver-free second-order central-upwind scheme. Central-upwind schemes have been originally developed for hyperbolic systems of conservation laws [28, 26, 30] and then extended and applied to hyperbolic systems of balance laws arising in modeling shallow water flows; see, e.g., [6, 13, 31, 25, 27, 29]. Here, we proceed along the lines of previous works in [8, 10, 11, 12] and develop a central-upwind scheme, which is well-balanced in the sense that it exactly preserves the geostrophic equilibria (1.11). This is achieved by performing a piecewise linear reconstruction of the equilibrium variables followed by well-balanced evolution, which uses modified central-upwind fluxes similar to those presented in . Preserving the positivity of and is another crucial property a good scheme must possess. This is enforced with the help of the so-called “draining time step” technique originally introduced in .
The paper is organized as follows. In §2, we present a general analysis of the TRSW system, first in the -plane, then in the equatorial beta-plane approximations. We obtain conditions of existence and uniqueness of the thermo-geostrophic equilibrium corresponding to arbitrary initial conditions, study the process of relaxation (adjustment) to such equilibrium, and get conditions of wave-breaking and shock formation. A well-balanced semi-discrete second-order central-upwind scheme for the TRSW system (1.8) is presented in §3. The proposed scheme is tested on a number of numerical examples reported in §4.
2 Analysis of the One-Dimensional TRSW Equations
In this section, we analyze the studied TRSW equations with the flat bottom topography, that is, with . We first obtain the results in the case of constant Coriolis parameter (§2.3–§2.5). To this end, we first rewrite the TRSW equations in the Lagrangian form (§2.1) and then analyze the existence of thermo-geostrophic equilibria for given initial conditions (§2.3), study the adjustment process (§2.4), and demonstrate the breakdown of smooth solutions (§2.5). We then extend these results to the equatorial case with , (§2.6), where the situation proves to be different in several aspects.
2.1 Lagrangian Formulation of the TRSW Equations
We first introduce Lagrangian positions of fluid parcels on the -axis, , with, correspondingly, and
. We then use the Euler-Lagrange duality in the description of the fluid motions, that is, the fact that at any time moment at each point in space we have a Lagrangian parcel, in order to rewrite equations (1.1b) and (1.1c) with as follows:
where we have used a dot notation for the Lagrangian time derivative. We notice that equation (2.1) expresses the Lagrangian conservation of the geostrophic momentum . The mass conservation equation (1.1a) in the Lagrangian formulation becomes
where is the initial thickness distribution and a prime notation is used for the partial derivative with respect to .
where and are the initial distributions of zonal velocity and buoyancy, respectively. For the sake of physical consistency, is assumed to be positive, that is, for all . Using (2.3) and (2.4) equation (2.2) can be rewritten in the following form:
It is often convenient to introduce the deviation of the fluid parcels from their initial positions: . Equation (2.5) then takes the following form:
2.2 Thermo-Geostrophic Equilibrium and Adjustment
It immediately follows from (2.6) that its stationary solution with no displacement of fluid parcels, that is, with exists if the initial data are in thermo-geostrophic equilibrium, that is, if
is satisfied for all . If , and do not obey (2.7), that is, if there is an initial imbalance
then the solution of the TRSW system evolves. One can show that a stationary solution minimizes the total energy of the system, which is the sum of the kinetic and potential energies:
As mentioned in §1, the system out of equilibrium tends to reach an equilibrium state. As the dissipation is absent in the studied TRSW model, this can be only achieved by evacuating energy through the emission of waves. This process will be called thermo-geostrophic adjustment. In mathematical terms, the adjustment process describes a solution of the Cauchy problem for (2.6) with radiation boundary conditions at spatial infinity.
2.3 Existence of the Adjusted State
The adjusted stationary state is described by the stationary part of equation (2.5), which can be rewritten in the following form:
Differentiating (2.9) with respect to results in
is a generalization of the Lagrangian pressure introduced in  in the context of standard RSW equations for the TRSW case.
and equation (2.10) can be rewritten in the following form:
is the potential vorticity. We would like to emphasize that the potential vorticity is a Lagrangian invariant in the 1-D TRSW equations, which is not the case in the case of two spatial dimensions; see .
Finally, we make two more changes of variables, and , and rewrite equation (2.13) in the canonical form:
Then, in the physically relevant case of monotone with constant asymptotics at infinities, which corresponds to a density/temperature front, the main configuration we are aiming at, existence and uniqueness of the decaying at infinities solution can be proved for nonnegative potential vorticity along the lines of .
The case of non-monotone requires a special consideration.
2.4 Linear Theory of Thermo-Geostrophic Adjustment
We now analyze the adjustment process in the linear approximation, assuming that the norm of the deviations of fluid parcels is small. The linearization of equation (2.6) with respect to results in
Let us split the dynamical variable into the “slow” and “fast” components: , where the slow component is defined as the time-average: . By averaging equation (2.14) in time, and assuming that both and remain bounded at all times, we obtain the following ODE for :
We then introduce a new variable and rewrite equation (2.15) in the canonical form:
Notice that the quantity reduces to the geostrophic potential vorticity in the corresponding equation obtained for the standard RSW model with a constant buoyancy studied in . In the physically interesting frontal case when is a compactly supported function and and have constant asymptotics at infinities, it was shown in  that a solution decaying at infinities exists, and it is unique for nonnegative geostrophic potential vorticity. Decay conditions are imposed by an obvious reason that fluid parcels should not move far away from the front, where they are already at equilibrium. Correspondingly, decaying solutions of (2.16) exist for (also see the nonlinear existence result obtained in §2.3).
If and are both constant, this is the 1-D Klein-Gordon equation describing gravity waves propagating at the surface of the shallow water layer. For nonconstant and , equation (2.17) describes the wave propagation over a variable background. For harmonic waves with , where c.c. stands for complex conjugation, equation (2.17) becomes
In the aforementioned frontal case, where and as , the far asymptotics of the solutions of (2.18) satisfy the following constant-coefficient ODEs:
and correspond to propagating inertia-gravity waves with the standard dispersion relations, obtained after applying the Fourier transform in space:
where is the wavenumber. Hence, the adjustment process corresponding to a solution of the Cauchy problem for (2.14) with zero initial and consists of the emission of inertia-gravity waves out of the initial unbalanced front, which evacuate an excess of energy and drive the system to an equilibrium state given by a solution of (2.16). The question, however, arises, whether some part of the wave signal could be trapped at the front and, thus, prevent the front from complete equilibration. In order to answer this question, we rewrite equation (2.18) as follows:
then multiply (2.21) by , where the star denotes the complex conjugation, and rewrite it in the following form:
Integrating this equation from to in , and assuming the decay of
far from the front, as we are looking for trapped modes, leads to the following estimate for the eigenfrequencies:
which shows that the frequencies are suprainertial, that is, , while in order to have trapped modes they should be sub-inertial, which can be shown following the lines of . It should be kept in mind, however, that the group velocity of near-inertial waves with frequencies close to is small, as follows from (2.20). This practically means that the portion of the “fast” component of the perturbation with will remain at the initial location for a very long time. We will see a manifestation of this fact in the results of numerical simulations reported below.
2.5 Breakdown of Smooth Solutions
The standard Saint-Venant system of shallow water equations is equivalent to the isentropic Euler equations of gas dynamics. Hence, shock formation is ubiquitous in this system. As shown in , including rotation into the Saint-Venant system does not prevent shock formation, but changes the breakdown conditions. In this section, we study the shock formation and properties of shocks in the TRSW equations, which are pertinent in the context of numerical simulations using shock-capturing finite-volume methods. In order to address these questions, we follow  and use the Lagrangian description in the case of constant initial , which can be always achieved by a “straightening” change of the independent spatial variable; see §2.3.
We rewrite the Lagrangian equations of motion in the following form:
where is given by (2.11) with the constant . In (2.22), and are dependent variables, while is not and needs to be determined from the conservation of the geostrophic momentum and potential vorticity :
It is easy to check that the system (2.22) can be rewritten in an explicit quasi-linear form:
compare with [49, equation (57)].
We now take (without loss of generality) and proceed along the lines of 
. The eigenvalues of the matrixare
and the corresponding left eigenvectors are. Hence, the Riemann invariants are and we have
Next, we differentiate (2.23) with respect to and obtain that the derivatives of the Riemann invariants, , satisfy the following PDEs:
Using the expressions for the Riemann invariants, we obtain , which implies that
This is a generalized Ricatti equation, which can be analyzed following , as it was done in  for the 1-D RSW model. Breakdown and shock formation correspond to reaching infinite values in finite time. Compared to the corresponding equations in , after obvious changes and in the latter, we observe the following two differences.
First, the first term on the right-hand side of (2.26), corresponding to vorticity, acquires an addition . Recall that if the initial vorticity is sufficiently negative, breakdown always takes place in the RSW equations. Here, in the TRSW model, it is the vorticity plus this new term (which depends on the initial distributions of buoyancy and thickness), should be sufficiently negative for the breakdown to take place.
Second, the breakdown is conditioned by the signs of the derivatives of Riemann invariants. As follows from their definitions, the overall sign of the derivatives depends not only on the signs of derivatives of and (as in the standard RSW model), but also on the sign of the derivative of , which makes a difference. It is worth emphasizing, however, that in practice it is difficult to discriminate the role of each contribution in the simulations, and specially designed initial conditions are required in order to do this; see .
2.6 Extension to the Equatorial Case
In this section, we discuss a possible extension of the results obtained in §2.3–§2.5 to the equatorial case with a variable Coriolis parameter , . We first reformulate the 1-D TRSW equations on the equatorial beta-plane in Lagrangian variables. Equations (2.1) and (2.2) take the following form:
and then by expressing , and in terms of their initial values, equation (2.28) reduces to
Already the inspection of (2.28) and (2.29) shows the fundamental difference from the -plane case: the dependence on meridional coordinate in the Coriolis parameter introduces higher powers of . This renders impossible the procedure used above in the -plane case in the demonstration of both existence and uniqueness of the adjusted state and breakdown. This procedure consisted of transforming, by differentiation, the original system of Lagrangian equations into a pair of equations for the - and -derivatives of (or a single equation for the space derivative in the stationary case), while eliminating the itself. Obviously, this is not possible on the beta-plane, which introduces serious technical difficulties (and even principal ones, like finding Riemann invariants for a system of three quasilinear equations). We therefore will not pursue these demonstrations below and will limit ourselves only by the linear analysis of the thermo-geostrophic adjustment on the equatorial beta-plane.
Introducing, as before, the deviations of the fluid parcels from their initial positions, we rewrite (2.29) as
split into slow and fast variables, and obtain the equation for the slow motion,
which is rewritten in terms of in the following form:
As in §2.4, a solution decaying at exists if the quantity in the square brackets is nonnegative. Notice, however, that the novelty of this expression with respect to its -plane counterpart is that it has an extra possibility to become negative if the initial flow is oriented westward (negative ) and sufficiently strong, which makes the factor negative. We will see that the sign of this factor, which is related to the absolute zonal momentum density on the beta-plane, also plays an important role in the dynamics of fast motions.
Again, as in §2.4, we obtain for the fast component
and then performing the Fourier transform in time results in
In the front/jet configurations, which are of primary interest, where and have constant asymptotics and tends to zero at , equation (2.31) becomes at the far left and far right sides of the front, respectively:
The crucial difference between these equations and (2.19) is appearance of as a coefficient. After rescaling with the corresponding equatorial deformation radii , and by the corresponding equatorial inertial periods , (2.32) takes a canonical form
where . If the natural for the equatorial region decay boundary conditions are imposed, equation (2.33) can be solved in terms of Gauss-Hermite (parabolic cylinder) functions obeying the equation
where is a number of zeroes of in . Hence, even in the absence of initial inhomogeneities in the buoyancy, thickness and zonal velocity, the fast component does not represent freely propagating inertia-gravity waves, but the waves are trapped at the Equator. The resulting eigenfrequencies correspond to the infinite zonal wavelength limit of the classical spectrum of the equatorial waves; see, e.g., . An important conclusion, following from this analysis, is that the fast component cannot be evacuated, like in the -plane case, but remains trapped at the Equator.
Another peculiarity of the equatorial adjustment is a possible appearance of additional trapped modes and even of instability for strong enough westward jets. Indeed, along the same lines as in §2.4, equation (2.31) leads to the following integral estimate for the eigenfrequencies:
which shows that not only the lowest nondimensional (according to the scaling above) eigenfrequency in (2.31) can be lower than the minimal frequency of the equatorial waves, meaning that these eigenmodes can be trapped inside the jet, but that the eigenfrequency squared can become negative (if is sufficiently negative), meaning imaginary eigenfrequencies and thus linear instability. Such instability is known in the RSW equations in the equatorial beta-plane (see , where it was analyzed along the same lines) as a symmetric inertial instability. The present analysis shows that it exists in the TRSW model as well. Obviously, the exponential growth of unstable modes corresponding to imaginary eigenfrequencies rapidly invalidates the linear approximation used in this section. However, it is known from other studies  that the growth of symmetric inertial instability modes in shallow water models leads to their breakdown and shock formation in the negative-vorticity (anticyclonic) part of the jet. We therefore expect a similar scenario here. The analysis of this instability in  shows that it appears at Rossby numbers of the order unity and Burger numbers below . Here, we define the Rossby number and Burger number as follows:
where is the maximum velocity of the jet, is its typical width, is the mean fluid depth, and is the mean buoyancy.
3 Well-Balanced Semi-Discrete Central-Upwind Scheme
In this section, we describe a semi-discrete second-order central-upwind scheme for the 1-D TRSW system (1.8
). To this end, we rewrite this system in a vector form as
We then derive a semi-discretization of (3.1)–(3.2) as follows. We divide the computational domain into a set of uniform cells , which are centered at . We denote the cell averages of the numerical solutions at time by and then integrate the system (3.1), (3.2) in space to obtain the following system of ODEs:
Here, are numerical fluxes, which typically depends on the reconstructed left- and right-sided point values of at the cell interfaces , and
are the approximations of the cell averages of the source term. When is taken to be a constant, we simply take . However, when a more realistic case of is considered, we use Simpson rule to obtain
where and are the one-sided values of at the cell interface; these values will be defined in §3.1 below.
Details on the computations of are provided in §3.2. For the sake of brevity, we will omit the dependence of all of the indexed finite-volume quantities on in the rest of this paper.
3.1 Well-Balanced Reconstruction
It is quite well-known that in order to derive a well-balanced scheme one has to perform piecewise polynomial reconstruction of equilibrium variables rather than the conservative ones; see, e.g., [8, 9, 10, 11, 12, 31, 25, 29]. We therefore reconstruct the equilibrium variables .
To this end, we first compute the values at the cell centers as follows. If the cell averages are available at a certain time level , then according to (1.9), one obtains
where is a total number of cells and can be computed using (1.10):
We notice that formula (3.4) can be rewritten in the following recursive way:
and then we apply the following second-order quadrature to the integral in (3.5), to obtain
Similarly, we can use a slightly different quadrature to obtain the point values of at the cell interfaces: