A Well-Balanced Central-Upwind Scheme for the Thermal Rotating Shallow Water Equations

by   Alexander Kurganov, et al.

We develop a well-balanced central-upwind scheme for rotating shallow water model with horizontal temperature and/or density gradients—the thermal rotating shallow water (TRSW). The scheme is designed using the flux globalization approach: first, the source terms are incorporated into the fluxes, which results in a hyperbolic system with global fluxes; second, we apply the Riemann-problem-solver-free central-upwind scheme to the rewritten system. We ensure that the resulting method is well-balanced by switching off the numerical diffusion when the computed solution is near (at) thermo-geostrophic equilibria. The designed scheme is successfully tested on a series of numerical examples. Motivated by future applications to large-scale motions in the ocean and atmosphere, the model is considered on the tangent plane to a rotating planet both in mid-latitudes and at the Equator. The numerical scheme is shown to be capable of quite accurately maintaining the equilibrium states in the presence of nontrivial topography and rotation. Prior to numerical simulations, an analysis of the TRSW model based on the use of Lagrangian variables is presented, allowing one to obtain criteria of existence and uniqueness of the equilibrium state, of the wave-breaking and shock formation, and of instability development out of given initial conditions. The established criteria are confirmed in the conducted numerical experiments.



There are no comments yet.


page 1

page 2

page 3

page 4


An adaptive well-balanced positivity preserving scheme on quadtree grids for shallow water equations

We present an adaptive well-balanced positivity preserving scheme on qua...

Steady States and Well-balanced Schemes for Shallow Water Moment Equations with Topography

In this paper, we investigate steady states of shallow water moment equa...

Adaptive Central-Upwind Scheme on Triangular Grids for the Saint-Venant System

In this work, we develop a robust adaptive well-balanced and positivity-...

A Well Balanced Reconstruction with Bounded Velocities and Low-Oscillation Slow Shocks for the Shallow Water Equations

Many numerical schemes for hyperbolic systems require a piecewise polyno...

An adaptive central-upwind scheme on quadtree grids for variable density shallow water equations

Minimizing computational cost is one of the major challenges in the mode...

Nonlinear mixed-dimension model for embedded tubular networks with application to root water uptake

We present a numerical scheme for the solution of nonlinear mixed-dimens...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

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., [50]. 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 [44], 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 [50].

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 [15] 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 [13]. A path-conservative approximate Riemann-problem-solver based scheme for the TSW system (1.2) with was developed in [40]. 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 [22], 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 [10]. 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 [3].

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 .

We note that equations (2.1) and (1.1d) can be immediately integrated. This results in


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 [49] in the context of standard RSW equations for the TRSW case.

In order to simplify equation (2.10), we introduce a “straightening” change of the space variable such that . In terms of the new variable, the mass conservation equation (2.3) reads as


and equation (2.10) can be rewritten in the following form:


where and

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 [44].

We note that equation (2.13) differs from the corresponding equation in the standard RSW studied in [49] by the presence of the factors.

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 [49].

Remark 2.1

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 [49]. In the physically interesting frontal case when is a compactly supported function and and have constant asymptotics at infinities, it was shown in [49] 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).

The PDE for the “fast” component of is obtained by subtracting (2.14) from (2.15), and reads as


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 [49]. 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 [49], 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 [49] 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 [49]

. The eigenvalues of the matrix


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

and thus


Finally, substituting (2.25) into (2.24) and using (2.12) result in


This is a generalized Ricatti equation, which can be analyzed following [18], as it was done in [49] for the 1-D RSW model. Breakdown and shock formation correspond to reaching infinite values in finite time. Compared to the corresponding equations in [49], 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 [4].

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:


while equation (2.3) does not change. As in the -plane case, equation (2.27) can be easily integrated:

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


The thermo-geostrophic balance of the initial conditions (2.7) with provides a trivial solution . If the imbalance (2.8) is small, we linearize (2.30) as in §2.4:

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., [50]. 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 [37], 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 [5] 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 [37] 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


where .

Similarly, we can use a slightly different quadrature to obtain the point values of at the cell interfaces:


It should be observed that the recursive formulae (3.6) and (3.7) require starting values. We first take , then compute using (3.7), and then set .

We also notice that the point values of (which are used in (3.6) and (3.7)) are obtained as in [29], namely, we take

which reduces to when

is continuous. Then, the bottom topography is approximated using a continuous piecewise linear interpolant

and then take