A robust well-balanced scheme for multi-layer shallow water equations

. The numerical resolution of the multi-layer shallow water system encounters two additional diﬃculties with respect to the one-layer system. The ﬁrst is that the system involves nonconservative terms, and the second is that it is not always hyperbolic. A splitting scheme has been proposed by Bouchut and Morales, that enables to ensure a discrete entropy inequality and the well-balanced property, without any theoretical diﬃculty related to the loss of hyperbolicity. However, this scheme has been shown to of-ten give wrong solutions. We introduce here a variant of the splitting scheme, that has the overall property of being conservative in the total momentum. It is based on a source-centered hydrostatic scheme for the one-layer shallow water system, a variant of the hydrostatic scheme. The ﬁnal method enables to treat an arbitrary number m of layers, with arbitrary densities ρ 1 , . . . , ρ m , and arbitrary topography. It has no restriction concerning complex eigenvalues, it is well-balanced and it is able to treat vacuum, it satisﬁes a semi-discrete entropy inequality. The scheme is fast to execute, as is the one-layer hydrostatic method. scheme, hydrostatic reconstruction, discrete entropy inequality.


Introduction
The multi-layer shallow water model is used to describe incompressible flows in the shallow water regime, in the situation where several layers with different densities can be identified. With rotating Coriolis force, the model can be used to describe large scale atmosphere or ocean flows. In the setting of this model, the following assumptions are needed: shallowness of the layers with respect to horizontal scales, small viscosity, hydrostatics, small slope for the topography, almost uniform velocities across each layer, uniform densities for each layer.
Suitable numerical methods for solving one-layer shallow water models have been developed in the last fifteen years, and in particular for getting the so called well-balanced property, that enables to exactly Key words and phrases. multi-layer shallow water, nonconservative system, nonhyperbolic system, well-balanced scheme, splitting scheme, hydrostatic reconstruction, discrete entropy inequality.
resolve steady states at rest. Well-recognized methods are those of [9] and [4] (that are indeed identical at first order), in the context of the Roe solver, and the hydrostatic reconstruction method of [3].
The two-layer shallow water model have been studied more recently, and involves much more difficulties, because of nonconservative terms and nonhyperbolicity. The only works that deal with it are [9], [2], [1], [6].
Here we introduce a variant of the splitting method of [6], that enables to correct the problem of wrong solutions with unexpected discontinuities found in [6]. The scheme involves a source-centered hydrostatic reconstruction scheme for the one-layer system, a variant of the hydrostatic reconstruction scheme. In contrast to the previous works, this new scheme allows to treat nonhyperbolicity, arbitrary number of layers, and arbitrary densities in a very robust manner, it is wellbalanced and satisfies a semi-discrete entropy inequality.

One-layer shallow water
Before writing the multi-layer model, let us first recall the main features of the one-layer shallow water model.
2.1. The system. The equations of the one-dimensional one-layer shallow water problem are: where h(t, x) ≥ 0 is the water height, u(t, x) ∈ R is the velocity, z(x) is the topography, and g > 0 (see Figure 1). The main properties are z(x) : topography h(t, x) : water level steady state at rest Figure 1. One layer that • The system has a convex entropy η = hu 2 /2 + gh 2 /2 + hgz with entropy flux G = ( η + gh 2 /2)u.
The steady states at rest are those for which u = 0 and h + z = cst. • The system is hyperbolic, except at resonant points (crossing of eigenvalues), defined by u = ± √ gh.
2.2. Numerical difficulties. The desirable properties of a numerical method for solving the one-layer shallow water system are of: • keeping the water height h nonnegative, • being able to compute dry areas where h = 0, • preserving the total mass, • maintaining steady states at rest (well-balanced property), • satisfying a discrete entropy inequality, • producing stable computations (no oscillations) for all data, including transcritical cases.

2.3.
Hydrostatic reconstruction scheme. The hydrostatic reconstruction scheme proposed in [3] satisfies all of the above properties (except the entropy inequality which is only semi-discrete), and is cheap computationally. With U = (h, hu) it can be written with classical notations (see [7]) where p(h) = gh 2 /2, F (U l , U r ) is a consistent numerical flux for the shallow water problem without source (z = cst), and the reconstructed states U * l , U * r are defined by

Multi-layer shallow water system
The one-dimensional multi-layer shallow water system can be written where h j ≥ 0, j = 1, . . . , m are the fluid depths, u j are the velocities, and z(x) is the topography. The constants g, are respectively the gravity and the densities of the fluids. The fluids 1, . . . , m are labeled from top to bottom. The situation is represented on Figure 2 for two layers. This system admits a convex entropy, and z(x) : topography steady state at rest Figure 2. Two layers thus we are looking for entropy solutions, satisfying (9) The system has the steady states at rest (10) Notice that if ρ 1 < · · · < ρ m , this reduces to while if ρ 1 = · · · = ρ m , (10) reduces to As for the one-layer shallow water system (m = 1), the numerical difficulties related to this system are positivity of the depths h j , with the possibility of treating out-or in-cropping ("drying" of one ore more layers), the exact preservation of the steady states at rest (wellbalanced property), and the property to have a discrete entropy inequality. Overall, the multi-layer system has extra difficulties which are the non-conservativity of the system (even for smooth topography z) that implies the lack of appropriate Rankine-Hugoniot relations, and the possibility of having complex eigenvalues (the system is not everywhere hyperbolic). This last property raises the question of the mathematical well-posedness of the Cauchy problem. The answer is not clear, but some hope is possible, see in particular [12] for a study on a related system. In physical terms the loss of hyperbolicity corresponds to the onset of the well-known Kelvin -Helmholtz (KH) instability, i.e. instability due to the strong velocity shear between adjacent fluid layers. It was shown in [11] that criteria of hyperbolicity loss and of KH instability coincide in two-layer shallow water with a rigid lid. This argument can be easily extended to the (present) free-surface case. The same line of argument is applicaple to the sheared continously stratified hydrostatic flows, as was shown recently in [10], where the so-called Howard -Miles criterion of instability coincides with that of hyperbolicity loss. The shear instabilities persist even in the presence of dissipation, i.e. friction between the layers. In the multi-layer model, the mixing is replaced by numerical dissipation mechanisms, however of unknown nature since the system is not hyperbolic. One can hope that such dissipative modelling gives a rough description of the missing mixing, in the same spirit as the idea that a shock replaces multivalued solutions in the case of wave-breaking in fluid dynamics. Thus, the mathematical difficulty of having a mixed system corresponds to well-established physical phenomenon, which reinforces the hope of properly dealing with the Cauchy problem.
Several attempts have been made in order to solve the multi-layer shallow water system. In [9] and subsequent papers of the Malaga-Sevilla school, the two-layer case is treated by a Roe type method. A special treatment is done in order to recover positivity, and a special treatment is performed for complex eigenvalues, making the scheme inconsistent in this case. A relaxation method is proposed in [1], with similar properties, but which is not able to treat outcropping, nor complex eigenvalues. In [2], the m-layer system is treated in the case ρ 1 = . . . ρ m , without restrictions on the eigenvalues and including drying. However, topography is not included, and the scheme is not extensible to arbitrary densities. None of these solvers is able to treat the general case with complex eigenvalues in a consistent way.

Splitting method
In [6] we introduced the operator splitting method for the multilayer shallow water system, that leads to solving successively (7) for each j, with time independent h k for k = j. This means that we solve successively m one-layer shallow water systems for U j = (h j , h j u j ), with topography Applying a one-layer solver with numerical fluxes F l , F r gives The advantages of the method are that • The scheme is well-balanced if the resolution of each subsystem is performed with a well-balanced scheme. This is true even in the case of constant densities ρ 1 = · · · = ρ m . • The scheme is nonnegative in heights and entropy satisfying if the resolution of each subsystem is performed with an entropy satisfying scheme, because the entropy of the system is an entropy for each of the subsystems, as explained in [6]. • No estimate of the (possibly complex) eigenvalues of the full system is necessary. In particular the scheme is always consistent, and one can use the CFL condition associated to each subsystem.

Sum scheme.
A variant of the splitting scheme is the "sum scheme", that can be defined as follows for two operators. Consider an ordinary differential equation Assume that we have two schemes for solving dU/dt + A(U) = 0 and dU/dt + B(U) = 0 respectively, given by Then the splitting scheme can be written (17) while the sum scheme is simply In [6], the splitting scheme (17) was used. The application of the sum scheme gives the same properties (entropy inequality, well-balanced property, nonnegativity...), and give very similar results. Note that in our context, applying the splitting scheme means that we solve (14) successively for each j, while applying the sum scheme means that we solve (14) simultaneously for all j.

4.2.
The problem of momentum conservation. The application of the splitting method of [6], using the one-layer hydrostatic solver, leads to wrong solutions, as seen on Figure 4. It happens also with the sum scheme. In [6] we thought the problem was due to only a semidiscrete entropy inequality. Here we show that the problem is in fact due to non-conservativity of total momentum. By summing up the momenta equations in (7) with weights ρ j , we deduce the total momentum balance that follows from the identity The equation (19) is conservative as soon as z is smooth. Since this conservative equation is physically meaningful, one would like to use a scheme that is also conservative in total momentum. At the discrete level, this means that one would like a discrete analogue of (20), telling that the discrete nonconservative terms add up to give a conservative contribution. A natural way of achieving this is by relying on the quadratic identity where indices l and r refer to left and right values with respect to a given interface between two cells, and h 1 and h 2 denote two different layers. Applying this several times leads to which is the desired discrete analogue of (20) (we denote ∆z = z r − z l ). Since the right-hand side of (22) is in conservative form (in the limit ∆z → 0, which corresponds to the assumption of continuous z(x)), it is enough in order to get the discrete conservation that the nonconservative term in the momentum equation (7) is discretized for each j as where ∆z j = z j r − z j l and is the topography seen by the j's layer in the splitting scheme. However, in order not to have different times involved, such as t n+1/2 in (17) (that would break the conservation), we need to apply the sum scheme (18), that involves only the data at time n. One can observe that the one-layer hydrostatic reconstruction method (3)-(6) gives a nonconservative contribution F l − F r = (0, g h l +hr 2 ∆z), explaining the failure of conservation of total momentum.
Note that the centered discretization (23) is consistent with the following Rankine-Hugoniot condition for the j's momentum where as usual [. . . ] denotes the jump of a quantity through a spacetime discontinuity, s denotes the speed of the discontinuity, and indices l and r refer to the values on each side. These Rankine-Hugoniot conditions are also those corresponding to the scheme of [9]. Thus we expect to get the same solution, even if no rigorous proof of convergence to generalized Rankine-Hugoniot conditions exists for nonconservative systems, see [8].

Source-centered hydrostatic reconstruction
As explained in the previous section, the sum scheme for the multilayer shallow water system satisfies the discrete conservation of total momentum whenever the one-layer solver used for U = (h, hu) with left/right numerical fluxes F l , F r , verifies the source-centered identity We can observe that this property (28) is true for Roe type solvers like those of [4], [9], because in this context the source discretization can be chosen arbitrarily. However this property is not true for the hydrostatic reconstruction solver, that we would like to keep because of its robustness. We provide below a correction to it in order to enforce the source-centered property, but without breaking the semi-discrete entropy inequality.

Dispersive correction.
In order to provide some intuition, at the continuous level, the idea is to find a dispersive correction to the shallow water system, that is neutral at the level of entropy. A parametrization of possible modified systems is as follows, with ∂ t z = 0, and µ(t, x), w(t, x) arbitrary functions. Then the entropy η = hu 2 /2 + gh 2 /2 + hgz satisfies showing that the correction is conservative in entropy. It also keeps steady states at rest.
The aim is to find a correction µ∂ x w of order ∆x (it will be indeed of order ∆z), in such a way that the nonconservative part of the momentum equation −µ(∂ x w)g∂ x (h+ z) exactly balances the noncentered part of the numerical source coming from the hydrostatic reconstruction scheme.
At the discrete level, omitting the standard shallow water terms, a dispersive correction can be performed with U = (h, hu) as and we take the momentum formulas where and X + ≡ max(0, X), X − ≡ min(0, X). The scalars δ i+1/2− , δ i+1/2+ parametrize the nonconservative part δ i+1/2− + δ i+1/2+ of the momentum correction. Writing (31) by interface, we deduce sufficient nonnegativity conditions by interface for h, Let us now write down discrete entropy inequalities. We decompose the entropy η = hu 2 /2 + gh 2 /2 + hgz as η(U, z) = K(U) + L(U, z) with Let us define J i+1/2 = (J 0 i+1/2 , J 1 i+1/2 ), and Then (37) say that U ′ i and U ′′ i have nonnegative first component. Writing A similar computation leads to (42) Therefore, one gets Using (34) and the fact that K is quadratic in hu at fixed h, (43) yields Next, we treat L(U, z), which is easier since it does only depend on h and z. We have for some ϑ L i+1/2 to be determined. Adding up (45) and (46) gives (47) Now, assume for a moment that the errors in (47) are both nonpositive. Then, from (36) we get showing that a discrete entropy inequality is satisfied with numerical entropy flux But since requiring nonpositive errors in (47) is too restrictive, we shall only look for a semi-discrete entropy inequality, which is the limit of (48) divided by ∆t when ∆t → 0. For this to hold it is enough that the errors in (47) divided by ∆t give nonpositive contributions in the limit ∆t → 0. These errors coming from (45) and (46), this yields the sufficient conditions or in other words (51) In order to be able to find some ϑ L i+1/2 satisfying (51), it is necessary and sufficient that the lower bound is less than the upper bound, which is Then one can choose for ϑ L i+1/2 the half sum of the lower and upper bounds in (51), Noticing that the CFL conditions for nonnegativity (37) are satisfied in the semi-discrete limit, we conclude the following.
Proposition 1. The scheme (31)-(32) with numerical fluxes J l/r (U l , U r , z l , z r ) satisfying J l/r = (J 0 , J 1 l/r ), and is semi-discrete entropy satisfying for the entropy η = hu 2 /2 + gh 2 /2 + hgz as soon as Then the entropy flux can be taken as with J l/r = (J 0 , J 1 l/r ) defined by (54), (55). We have to take δ l , δ r so that J 1 r − J 1 l = −(δ l + δ r ) will kill the source error between the hydrostatic source and the desired centered source, that gives (28). Then we define J 0 in order to have equality in (56). Notice that even if the construction (58)-(59) looks possible for any scheme instead of F HR l/r , the issue is in the possibility to compute J 0 from (56), which somehow involves the property that h r + z r − h l − z l is in factor in δ l and δ r , and thus must be factorizable from the right-hand side of (59). Indeed, the well-balanced property ensures that whenever u l = u r = 0 and h r +z r −h l −z l = 0, the right-hand side of (59) vanishes. But this does not imply that h r + z r − h l − z l can be factorized, unless (F HR r − F HR l ) 1 is independent of u l , u r . This is wrong in general, but this is true for the hydrostatic solver. Indeed one can check that One can observe that Thus one would like to take δ l = 1 2 (1 + θ)g(h r − h l + ∆z) κ, δ r = 1 2 (1 − θ)g(h r − h l + ∆z) κ, for some −1 ≤ θ ≤ 1, and J 0 = 1 2 (1 + θ)u l + (1 − θ)u r κ. However, in order for the CFL condition (37) to be finite, we need that |J 0 | ≤ C min(h l , h r ). But one can check that this is not true, and only (63) is valid.
Therefore, our choice is the following. We define and take so that the semi-discrete entropy condition (56) is satisfied. The value of θ is not very important, and we have chosen an upwind formula, Now, denote by a HR (U l , U r , ∆z) the speed involved in the CFL condition for nonnegativity of the hydrostatic reconstruction scheme, and by a J (U l , U r , ∆z) the one associated to the correction scheme (31)-(32). Indeed, according to (37), one can take We recall that according to [3], a HR (U l , U r , ∆z) = a F (U * l , U * r ) where a F stands for the speed associated to the homogeneous solver F . Then, the CFl condition for nonnegativity for the scheme (58) is ∆t with a SCHR (U l , U r , ∆z) = a HR (U l , U r , ∆z) + a J (U l , U r , ∆z). Proof. Recall that the hydrostatic reconstruction scheme satisfies (i)-(iv). Then, according to (63), (64), (65) one has which yields (i) and (iii). Property (ii) follows from the above construction (by applying Proposition 1), as well as property (iv). Finally, for (v), we observe that δ l + δ r = g(h r − h l + ∆z)κ. Thus, whenever κ = κ, the identity (60) gives (59), proving (28). Therefore, the only restriction is that κ = κ, which according to (64), (63), holds as soon as This inequality is what we mean by "away from sharp variations of h close to vacuum".
In practice, this scheme gives numerical results for the one-layer system that are almost identical to the usual hydrostatic reconstruction scheme, except that it needs slightly smaller timesteps because of the restricted CFL condition. This can be improved by diminishing a bit the constant 5/2 in (64) (take for example 3/2), at the cost of losing a bit of source-centering for sharp variations of h close to vacuum.
In order to treat the one-layer shallow water system with transverse velocity v solving ∂ t (hv) + ∂ x (huv) = 0, (71) one has to apply also the correction J . This means that we have to complete (31) with the passive advection scheme with (73) Then the discrete entropy inequality associated to the entropy hu 2 /2 + hv 2 /2 + gh 2 /2 + hgz holds with an additional numerical entropy flux Equation (71), completed with (1), corresponds to the one and half dimensional one-layer shallow water model. For treating further zeroorder right-hand sides in the momentum equation (such as Coriolis force) in a well-balanced setting, one can apply the apparent topography method of [5], see [7] for details.

5.3.
Sum scheme for the multi-layer system. From Theorem 1 we get directly Theorem 2. For the multi-layer shallow water system, the sum scheme (14), (13) using the source-centered hydrostatic reconstruction scheme on each layer is (i) consistent for smooth solutions, (ii) semi-discrete entropy satisfying, (iii) well-balanced, (iv) conservative in h j for each layer j, and nonnegative under the CFL conditions (68), (69), (67) for each j, (v) conservative in total momentum away from sharp variations of h j close to vacuum.
It is noticeable to remark that for a sum scheme if each scheme A j is stable (in the sense of invariant domain or fully discrete entropy inequality) under a CFL condition ∆t a j (U n ) ≤ 1, then the sum scheme is stable under the CFL condition ∆t m j=1 a j (U n ) ≤ 1.
However, since in our case the height h j evolves only through the term A j (and not through A k for k = j), the nonnegativity is true under the weak CFL condition ∆t max j=1,...,m as stated in (iv), while by nature the semi-discrete entropy inequality does not need any CFL condition (because it involves only the limit ∆t → 0). Therefore, even if one would expect that for a really stable sum scheme, the CFL condition (76) is needed, in practice all our computations with a large number of layers show that (77) is sufficient, at the cost of putting a factor 1/2 in the right-hand side, or a factor 1/4 in two dimensions. Note that the one and half dimensional multi-layer model can be treated similarly by applying the sum scheme, leading to additional resolutions of (71) for each layer. Then, fully two-dimensional finite volume simulations can be performed as usual by applying interface one and half dimensional solvers.

5.4.
Upwinding dispersive correction. The SCHR scheme obtained in the previous paragraph for the multi-layer shallow water system shows good numerical results, except for "upwind data" for which all eigenvalues have the same sign. For these we observe some important oscillations, as shown on Figure 3. These oscillations do not increase when refining, thus are not instabilities, but rather dispersive errors. An explanation of this phenomenon is that in the upwind situation, the scheme differs from the standard upwind scheme evaluating all the data on the upwind side, because the apparent topography z j = z + k>j h k + k<j ρ k ρ j h k is treated as a steady variable (instead of upwind).
A correction can be done at the level of the one-layer solver, by changing the definition (64) of κ into with κ = κ + (∆z) + − κ min 1, 1 4 The idea is that for Froude number F r = u/ √ gh small, κ = κ and the scheme is as previously, while for large Froude number, κ = (∆z) + (F r positive, i.e. positive eigenvalues), or κ = (−∆z) + (F r negative, i.e. negative eigenvalues).
The reason for putting these values is as follows. Consider the case of positive eigenvalues. Then the HR scheme (5)-(6) gives the numerical fluxes, if ∆z ≥ 0, while if ∆z ≤ 0, These values have to be compared to the "ideal" upwind values Taking into account (60), this yields the errors In (83) we can neglect g(h r − h l + ∆z) κ since it is of higher order (recall that 0 ≤ κ ≤ |∆z|/2). Now it remains to take into account the contribution of the correction fluxes (54)-(55) into (58). Since in (65) δ l and δ r can also be neglected, it remains J 0 = κu l , and (55) gives J 1 = κu 2 l . Thus (J 0 , J 1 ) = κ(u l , u 2 l ). In order for this to compensate (83), we need to have κ = (∆z) + , which justifies (79). The case of negative eigenvalues is similar.

Numerical tests
We only intend to show here that our approach works on the basic tests proposed in [6]. More involved tests (non-hyperbolic regime, more than two layers, second-order, two dimensions, Coriolis force) will be reported in other places.
In all our tests we take g = 9.81 with two layers, the one-layer UDCHR scheme of Subsection 5.4 being used into the sum scheme. The homogeneous solver F of (5) is taken as the Suliciu relaxation solver (the HLL solver gives similar results). A second-order two-step scheme in time is applied (Heun scheme). The weak CFL condition (77) is implemented with right-hand side 1/2. We use the first-order scheme in space, but nevertheless for Tests 1, 2, 4, 5, a reference solution computed at second-order with a very fine grid is plotted. Except for Test 1, the SCHR scheme gives similar results as the UDCHR scheme. 6.1. Test 1. It is an upwind Riemann problem with ρ 1 /ρ 2 = 0.98, z = 0, u 1 = u 2 = 2.5, x ∈ (0, 1) and the discontinuity is taken at x = 0.5. The final time is t = 0.05 and we use 100 points. Figure 3 shows the interface z +h 2 . For comparison the SCHR scheme is plotted also (dispersive oscillations).

Test 2.
It is a centered Riemann problem, ρ 1 /ρ 2 = 0.7, z = 0, u 1 = u 2 = 0, x ∈ (0, 10) and the discontinuity is taken at x = 5. The final time is t = 1 and we use 500 points. Figure 4 shows the interface. The result with the HR scheme is plotted also (unphysical discontinuity). The UDCHR scheme gives the same solution as the Roe method, indeed this test was proposed in [9].
h 1 + h 2 = 1, x ∈ (0, 10), the final time is t = 1 and we use 1000 points. Figures 5 and 6 show the benefit of using the UDCHR scheme.
x ∈ (0, 1) and the discontinuity is taken at x = 0.5. We use Neumann boundary conditions, the final time is t = 0.32 and we use 100 points. The interface is plotted on Figure 7.

Conclusion
We have proposed a simple uncoupled method to solve the multilayer shallow water system. Its features are: • It keeps the water heights nonnegative, and it is able to compute with dry areas, • It is well-balanced, • It satisfies a semi-discrete entropy inequality, • It works for an arbitrary number of layers, arbitrary densities, and arbitrary topography, • It can deal with complex eigenvalues consistently, indeed the eigenvalues of the system are never computed, • It runs very fast since it needs only to run a one-layer solver for each layer.   We have performed simple tests in order to validate the method, that gave good results, showing the exceptional robustness of the scheme. More involved tests will be presented in other publications, including [13].