A Central-Upwind Scheme for Two-layer Shallow-water Flows with Friction and Entrainment along Channels

We present a new high-resolution, non-oscillatory semi-discrete central-upwind scheme for one-dimensional two-layer shallow-water flows with friction and entrainment along channels with arbitrary cross sections and bottom topography. These flows are described by a conditionally hyperbolic balance law with non-conservative products. A detailed description of the properties of the model is provided, including entropy inequalities and asymptotic approximations of the eigenvalues of the corresponding coefficient matrix. The scheme extends existing central-upwind semi-discrete numerical methods for hyperbolic conservation and balance laws and it satisfies two properties crucial for the accurate simulation of shallow-water flows: it {\it preserves the positivity} of the water depth for each layer, and it is {\it well balanced}, i.e., the source terms arising from the geometry of the channel are discretized so as to balance the non-linear hyperbolic flux gradients. Along with the description of the scheme and proofs of these two properties, we present several numerical experiments that demonstrate the robustness of the numerical algorithm.


Introduction
In this paper we present a new high-order numerical scheme for simulating two-layer shallow-water flows along channels with a bottom topography and varying width (see Figure 1). These flows are characterized by a large horizontal length scale relative to their depth and are commonly observed in nature -e.g., channel flows, straits, mountain passes; and their modeling and simulation have applications in flood control, coastal engineering, or environmental assessment. Applications of this model include the study of internal waves observed in rivers and possibly caused by shear flow instabilities [Apel et al., 1975]. Shallow-water flows are typically modeled by the Saint-Venant equations, a hyperbolic balance law that results from the depth averaging of the Euler equations. The model poses various mathematical and computational challenges: the effects of the channel geometry on the flow are described by source terms that need to be discretized consistently with the existence of equilibrium solutions and the onset and propagation of discontinuities -e.g., hydraulic jumps. The dynamics of the interface between the two layers render non-conservative products of the flow variables and their space derivatives whose correct discretization is crucial for determining the location and propagation of these discontinuities. The various challenges that the Saint-Venant equations pose for simulating shallow-water flows are well known and have been studied extensively. The non-conservative terms may lead to the loss of hyperbolicity, and the common occurrence of steady states in geophysical flows require well-balance numerical schemes capable of capturing and resolving accurately steady-state solutions of the PDE model. These challenges become increasingly difficult to address as the flows and the corresponding PDE models that describe them increase in complexity, and a significant effort over the past couple of decades have lead to the development of numerical schemes for simulating a wide variety of flows. Numerous schemes have been proposed for the simplest case of one-layer flows along channels with constant width: a positivity preserving kinetic scheme capable of preserving the steady state of rest is presented in [Perthame and Simeoni, 2001], and in [Audusse et al., 2004] a finite volume scheme with similar properties is devised using hydrostatic reconstruction to capture steady-state solutions. The authors of [Khan and Lai, 2014] proposed a discontinuous Galerkin method. And [Kurganov and Levy, 2002, Bollermann et al., 2011 introduce positivity preserving well-balance central-upwind schemes for these flows. The approach suggested to achieve wellbalance and positivity within the central-upwind framework is extended in [Balbás and Karni, 2009] to simulate flows along channels with varying cross-sections using a central scheme, a type of flows also addressed in [Garcia-Navarro and Vazquez-Cendon, 2000] using an upwind Roe-type scheme. And the authors of this paper proposed a new central-upwind scheme for one-layer flows along channels with arbitrary geometry in [Balbás and Hernandez-Duenas, 2014]. Several schemes have also been devised for simulating two-layer flows.
For flows along straight channels, the authors of [Abgrall and Karni, 2009] propose to address the conditional hyperbolicity of the model with a relaxation approach that simplifies greatly its eigen structure. A similar approach is presented in [Chiapolino and Saurel, 2018] where a strictly hyperbolic formulation with pressure relaxation is obtained by considering weak compressibility of the phases of the flow. And in [Castro-Díaz et al., 2011] friction between the layers is added to prevent the loss of hyperbolicity. The central-upwind scheme in [Kurganov and Petrova, 2009a] provides a robust approach for flows along straight channels that is simple to implement, and the work in [Balbás and Karni, 2013] extends the central framework to channels with rectangular sections of varying width. For flows along channels with arbitrary geometry, the authors of [Castro et al., 2004] extended with great success (and high impact in the field) the Q-scheme for hyperbolic systems with source terms previously introduced in [Castro et al., 2001]. Some of these schemes for shallowwater flows have been employed -and new ones created-for simulating, studying or recreating geophysical flows from the real world like gravitational currents or strait flows. These typically pose additional challenges such as the handling of non-symmetrical channel geometries (described by bathymetry data) or accounting for phenomena like entrainment. For some examples, we refer the reader to the works presented in [Adduce et al., 2012, Apel et al., 1975, Castro et al., 2004, Castro et al., 2007, Kim and LeVeque, 2008.
In this paper we follow on our previous works on central schemes for one-layer flows along channels with arbitrary geometry, [Balbás and Hernandez-Duenas, 2014], and two-layer flows along channels with varying rectangular cross-sections, [Balbás and Karni, 2013]. We propose a new high-order central-upwind scheme to compute two-layer shallow-water flows along channels with arbitrary geometry that incorporates the treatment of friction and entrainment terms. These new terms allow us to simulate more realistic flows and to assess the limitations of the model and the scheme. In order to understand and address the challenges posed by the model and its limitations, we present a detailed analysis of the hyperbolic PDE model, with special emphasis on the conditions that lead to the loss of hyperbolicity. To this end, we derive rigorous asymptotic approximations of the eigenvalues of the quasilinear form of the model and, for the sake of completeness in the analysis, we prove the existence of an entropy function and an entropy inequality that physically relevant weak solutions must satisfy.
In order to address the various numerical and computational challenges we propose a central-upwind scheme based on the semi-discrete central-upwind schemes for hyperbolic conservation laws of Kurganov, Noelle, and Petrova, [Kurganov et al., 2001], characterized by their simple implementation and robustness.
The proposed numerical scheme evolves the cell averages of the flow variables with second order accuracy, and their implementation requires four main ingredients: a non-oscillatory reconstruction of point values from cell averages that preserves the positivity of the water depth, an evolution routine to advance the solution in time, estimates of the largest and smallest eigenvalues of the system, and the discretization of source terms and non-conservative products that balance the hyperbolic fluxes so as to recognize steady states at rest and add accuracy to the computation of flows near non stationary steady states.
The paper is structured as follows. In §1 we present the model, its properties, and the challenges that these properties pose for computing numerical solutions. In §2 we describe the proposed numerical scheme and prove that it preserves the positivity of the water height, and it is well-balanced, i.e., it recognizes and preserves the steady state of rest. Numerical solutions for a variety of flow regimes are presented in §3, validating the scheme's accuracy and robustness and demonstrating its ability to simulate a wide range of flows.

The Model and its Properties
The model for the two-layer shallow-water flows in channels is taken from [Castro et al., 2004]. After adding source terms due to friction and entrainment, it has been re-written as Here h 1 , and h 2 denote the depth of bottom and top layers respectively, u 1 , and u 2 the cross-sectional velocities; g the acceleration of gravity, B(x) describes the bottom topography, and σ(x, z) the width of σ(x, z) dz are the cross-sectional wet areas in each layer; w 1 = B + h 1 denotes the total elevation of the internal layer and w 2 = B + h 1 + h 2 that of the external layer; and Q 1 = A 1 u 1 , and Q 2 = A 2 u 2 are the flow rates or discharges for the internal and external layers respectively. Furthermore, σ B (x) = σ(x, B(x)), σ 1 (x, t) = σ(x, w 1 (x, t)), σ 2 (x, t) = σ(x, w 2 (x, t)) denote the channel's width at the bottom topography, and at the internal and external layers respectively. We note that σ 1 and σ 2 depend both on space and time since they are evaluated at the internal and external layers.
The ratio of densities is denoted by r = ρ 2 /ρ 1 ≤ 1. The vertically integrated hydrostatic pressure of the upper layer is given by and treats the internal layer as a moving topography. The hydrostatic vertically integrated pressure in the internal layer has to account for the contribution of the pressure exerted by the upper layer on it, and it is given by The source terms I 1 , and I 2 correspond to the vertically integrated pressure terms due to width variation, and are given by For piecewise trapezoidal channels, the integrals involved in the computation of the cross sectional areas coincide with the corresponding approximations given by the trapezoidal rule.
The friction terms are calculated as for the internal and external layers respectively. Here n i and n b are the Manning roughness coefficients for the interface and bottom respectively. The hydraulic radius, R, is defined as the ratio between the wet area and the wetted perimeter, and it is given by This friction terms in equation (5) have been adapted for the two-layer case from the formula used in [Khan and Lai, 2014] for one-layer non-rectangular channels (Ch. V, p. 83). See [Kim andLeVeque, 2008, Castro-Díaz et al., 2011] for other approaches. We have adapted to our model the expression of the source term due to entrainment for flows along channels with constant width found in [Adduce et al., 2012]. It is given by where V e is the entrainment velocity. Furthermore, we set S e proportional to A 1 to suppress entrainment in the absence of the heavier fluid in the internal layer. The cross sectional area in the external layer (A 2 ) is not small in the numerical tests in Section 3.5 where entrainment is included.
Equations (1a) and (1c) express conservation of mass. The momentum equation (1d) treats the elevation of the internal layer w 1 as a moving topography, resulting in a balance law with non-conservative products.
The momentum equation (1b) of the internal layer has a momentum exchange of the external layer on it, represented by non-conservative products. We note that in the limit h 2 → 0 we recover the one layer shallow-water system. On the other hand, when σ = σ(x) is independent of height (straight vertical walls), the pressure terms in (2) and (3) become p 1 = g σh 2 1 2 + rgσh 1 h 2 and p 2 = g σh 2 2 2 , and the source terms in (4) become resulting in the model presented in [Balbás and Karni, 2013] in the absence of friction and entrainment.
Source terms in hyperbolic balance laws that do not depend on the derivatives of the solution variables do not modify the Rankine-Hugoniot jump conditions of weak solutions. However, the source terms in system (1) include rgh 2 ∂ x A 1 and gh 2 σ 1 ∂ x w 1 with non-conservative products that depend on derivatives of the wet areas. The jump conditions are modified by such terms leading to theoretical and numerical challenges [Abgrall and Karni, 2010]. A theory of non-conservative products can give us a notion of weak solutions in those cases [Dal Maso et al., 1995]. Following the ideas in [Kurganov and Petrova, 2009b] implemented for 2D and 1D two-layer shallow-water flows (with no width variation), we rewrite system (1) in a more convenient way to treat the non-conservative products. To that end, we decompose the hydrostatic pressures as where w 2 = B + h 1 + rh 2 . Rewriting system (1) in terms of these, we obtain the alternative model We note that system (7) still has non-conservative products and it is equivalent to system (1) only when the associated solution is smooth. Therefore, both systems (1) and (7) present similar theoretical challenges.
However, the non-conservative products in system (7) have w 2 and w 2 as factors, which are approximately constant in many situations where the rigid-lid approximation is valid, such as the case of internal waves. In the ideal case where such factors are constant, however, those non-conservative products become conservative, and system (7) therefore minimizes the numerical challenges when simulating two-layer flows. If the spatial derivatives of A 1 and A 2 in the non-conservative products are treated as flux gradients, then any consistent second order reconstruction of w 2 and w 2 leads to robust numerical schemes. Details are given in Section 2.
Other efforts to add stability to the two-layer shallow-water equations includes the three-layer approximation, as described in [Chertock et al., 2013].

The Quasilinear Form
For one-layer shallow-water systems, one can find explicit expressions involving a speed of sound quantity.
For two-layer flows, however, the eigenvalues have no explicit formulas but can be approximated by explicit expressions under suitable conditions. Nonetheless, the definition of speed of sound for each layer will still be helpful. For the external layer, this quantity is analogous to the one layer case, and it is given by which, in practice, makes the internal layer a moving topography. The speed of sound for the internal layer depends on the width of both layers and it is given by where = 1 − r.
Similarly to channels with straight walls, the system can be written in quasi-linear form. The details are in the following proposition.
Proposition 1. Either system, (1) or (7) can be written in quasilinear form, and where Proof. The derivation of the quasilinear form follows from the fundamental theorem of calculus and the relations

Asymptotic Approximations of the Internal Eigenvalues as → 0
The characteristic polynomial of the matrix M in (10) is Due to the last term, the eigenvalues have no explicit simple formulas. The following proposition provides us with a first order approximation for the internal eigenvalues.
Proposition 2. The two eigenvalues associated with the internal layer satisfy is the convective velocity.
Proof. Expanding the characteristic polynomial we get Dividing by gA1 σ1 c 2 2 , we get The conditions for vanishing eigenvalues give us the composite Froude number where c 1 and c 2 are given, respectively, by equations (9) and (8). We note that c 2 1 gA1/σ1 → σ1 σ2 as → 0. So, the first two terms in the above equation are order 1 and The roots of the above quadratic polynomial are those given in (15).

Asymptotic Approximations of the External Eigenvalues as
Approximate explicit formulas of the external eigenvalues can obtained when u 2 − u 1 and are both small.
In that case, the fluid behaves like a one layer shallow-water flow. The details are given in the following proposition.
Let us assume that both parameters go to zero at the same rate. The external eigenvalues satisfy where are the vertically averaged velocity and the speed of sound of a one-layer shallow-water system with total cross-sectional area A 1 + A 2 and surface width σ 2 .
Proof. In what follows, we assume the parameters A 1 , A 2 , u 1 , σ 1 , σ 2 are all fixed. The only parameter that varies is u 2 = u 1 + δ 1 , where δ 1 = δc. We look for an expansion of the form with λ o independent of and δ 1 , and λ 1 = O(δ 1 ). Then Collecting order 1 terms we get Collecting order δ = O( ) terms, we get with solution Substituting λ o and λ 1 we get as desired.
The approximation for the internal (15) and external (17) eigenvalues indicates that the system is conditionally hyperbolic with an approximated condition given by As documented in [Abgrall and Karni, 2009], loss of hyperbolicity is associated with shear layer Kelvin-Helmholtz instabilities of the internal layer, and in such case the proposed balance law (7) is no longer valid to describe the flow.

Eigenvalue Bounds
Assuming we are in the hyperbolic regime where all the eigenvalues are real, and denote by γ ± 1 , γ ± 2 the roots given by The following proposition gives us the eigenvalue bounds. Such bounds will be used to define the one-sided local speeds, which in turn are used in the numerical scheme of Section 2.
Proposition 4. Let us assume that the eigenvalues of the coefficient matrix A in equation (10) are all real.
The proof of the previous proposition can be adapted to the present case of channels with arbitrary geometry from the original proof in [Abgrall and Karni, 2009]. We note that the values of γ ± k , k = 1, 2 are always real, regardless of the hyperbolicity condition.  (17) and (15) and the bounds in Proposition 4 are shown for different values of δ = from 0 to 0.5. Here A 1 = 1.5, A 2 = 2, u 1 = 1, σ 1 = 1.4, and σ 2 = 2 We have verified the approximations in equations (17) and (15) numerically for specific values and have found that the approximations are excellent even if and δ are not too small. Figure 2 shows the eigenvalues, the approximations (17) and (15)

Steady-state Solutions
System (1) admits non-trivial steady-state solutions. The following proposition characterizes conditions satisfied by general smooth steady states. We also describe the internal waves, consisting of a moving internal layer at equilibrium and an external layer at rest, as well as steady states at rest for both layers.
Proposition 5. In the absence of friction and entrainment (n i = n b = 0, V e = 0), smooth steady-state and steady states of rest that satisfy Proof. We use the relation in (12) to obtain Taking the difference between the flux gradients and the source term in the external layer we have For the internal layer, we notice that and using which concludes the proof for the general smooth steady states.

Entropy Functions
The existence of entropy functions and entropy inequalities have shown to be helpful in choosing the correct weak solution [Oleinik, 1957, Audusse et al., 2004. A numerical scheme that satisfies a fully discrete entropy inequality can be found in [Bouchut and de Luna, 2008]. To conclude this section and for the sake of completeness in the analysis of the hyperbolic model, in the following proposition we prove that system (1) admits one entropy function and entropy inequality that accounts for the contribution of both layers.
Proposition 6. In the absence of entrainment (V e = 0), system (1) is endowed with an entropy function , and its physically relevant solutions satisfy the entropy inequality Proof. Adding friction or viscous terms is usually the first step in deriving entropy inequalities. However, system (1) already admits friction terms. Without loss of generality, we assume that the Manning coefficients are positive. Using the fundamental theorem of calculus we get the relation Multiplying the momentum equation of the external layer by u 2 , the resulting equation can be rewritten as Subtracting the two equations and using ∂ t A 1 = σ 1 ∂ t w 1 , we get On the other hand, the pressure in the internal layer satisfies Multiplying the momentum equation of the internal layer by u 1 , the resulting equation can be rewritten as Subtracting the two equations we get Finally, using equations (25) and (27) and defining E = E 1 + rE 2 , we get which concludes the proof.

Numerical Scheme
In this section we describe a central-upwind scheme to approximate the solutions of system (7). The scheme must satisfy the well-balance and positivity-preserving properties. See [Balbás and Karni, 2013, Kurganov and Petrova, 2009b, Kurganov and Tadmor, 2000] for more details on central and central-upwind schemes. We extend the scheme here to two-layer flows along channels with arbitrary cross sections. The non-conservative products on the right hand side of (7) and the more complex geometry of the channels considered here makes the extension of these schemes more challenging than those presented in the above references.
The system now takes the form

Discretization of the Channel's Geometry
In order to arrive at a second order scheme for (28) -(29), we shall first address the discretization of the channel geometry given by the width function σ(x, z) and the bottom topography B(x). Since our approximate solution of the system (28) -(29) will be realized as the cell averages W j of the flow quantities over the grid cells I j := [x j− 1 2 , x j+ 1 2 ], where ∆x is the spatial scale in the direction of the flow, x j± 1 2 = x j ± ∆x 2 and x j is the center of the grid cell, we begin by sampling the bottom topography and width of the channel at the points x j+ 1 2 . From these samples, we define the cell average of the bottom topography as And similarly, we define the average width functions Using these, we build a piecewise trapezoidal approximation of the cross-sections of the channel at each cell where ∆z is a fixed spatial scale in the vertical direction (flow height).
These piecewise (bi-)linear functions allow us to approximate the spatial derivatives of B(x) and σ(x, z) with centered divided differences and are consistent with the second order accuracy sought for the scheme (see left panel of Figure 3). Unless otherwise noted, in the numerical experiments we use a resolution of ∆z = 0.01.

The Semi-discrete Central-upwind Scheme
With the above partition of the solution grid and discretization of the channel geometry, we denote by v j (t) the computed cell average of W(x, t) over the cell I j , v j (t) = 1 ∆x Integrating equation (7) over each cell I j , we obtain the semi-discrete formulation which is approximated by The flux at the cell interfaces, F(W(x j± 1 2 ), t), is approximated by the numerical flux H j± 1 2 (t) in [Kurganov et al., 2001] Here, W ± j± 1 2 (t) is recovered from the cell averages via a non-oscillatory piecewise polynomial reconstruction Furthermore, the one-sided local speeds in this scheme are approximated using the eigenvalues' bounds given by (21) is automatically satisfied given the expressions in equation (21), which is important for the positivitypreserving property.

Positivity Preserving Non-oscillatory Reconstruction
In order to recover the interface point values W ∓ j± 1 2 (t) from the cell averages v j (t), we seek a piecewise polynomial reconstruction at cell j given by with the limited slopes W j calculated as, [van Leer, 1997], This non-oscillatory reconstruction is applied directly to the discharges Q 1 and Q 2 . However, in order to ensure the positivity of the water layers and recognize steady states of rest (23), the point values of the areas A 1 and A 2 are obtained by first reconstructing the total elevations of the internal and external layers, w 1 = B + h 1 and w 2 = w 1 + h 2 respectively. These, unlike h 1 and h 2 , will remain constant for the steady state of rest, so their reconstruction will also render constant point values within and across grid cells, and help us achieve well balance as well as prevent the onset of spurious oscillations.
In order to reconstruct the interface point values of A ± 1,j∓ 1 2 and A ± 2,j∓ 1 2 from their cell averages A 1,j and A 2,j , we proceed as follows: • Determine the values z l1 and z l2 along the partition {z l } of the vertical axis such that and • Find the height δw 1,j of the next trapezoid in the discretization of the j th cross-section of the channel so that and set w 1,j := z l1 + δw 1,j .
Note that the left hand side of (43) and (44) are, respectively, increasing functions of δw 1 and δw 2 , and because of the linearity of σ j (z), the equations are quadratic. Therefore, δw 1,j and δw 2,j correspond, each, to one of the solutions of a quadratic equation. As an alternative, we may also determine their values using bisection so as to pick the right solution, which works well here because the cross sectional area is an increasing function of height.
σ(x, z) where σ j∓ 1 2 (z) = σ(x j∓ 1 2 , z l ) + is a piecewise trapezoidal approximation of the channel cross-section at the cell interfaces x j∓ 1 2 . These guarantee the positivity of the water layers and will remain constant across cell interfaces for steady states of rest, where the total elevation of the water layers, w 1 = h 1 + B and w 2 = w 1 + h 2 , remain constant.

Regularization of Flow Velocity and Discharge for Small A 1 , A 2
Although the positivity of depth in each layer is preserved, these point-values may still be very small and may lead to large values of the velocity of the flow, u 1 , u 2 . We can prevent this with the regularization technique suggested in [Kurganov and Petrova, 2007], for k = 1, 2. The discharges Q ± k,j± 1 2 are then recalculated as to ensure conservation. The value of δ A was empirically determined, usually choosing δ A = 10 −12 in this paper.

Well Balance
Flows described by the balance law (7) admit steady-state solutions of rest satisfying (23). Numerical schemes for approximating the solutions of (7) that are capable of capturing these equilibrium solutions are said to satisfy the well balance property.
The authors of [Kurganov and Levy, 2002] and [Kurganov and Petrova, 2007], where shallow-water models for flows along channels with constant width is discussed, propose to reconstruct the variables w 1 and w 2 , which remain constant for steady states of rest, and not the heights h 1 and h 2 , that, despite the water layers remaining flat, may vary across cells due to changes in the bottom topography. While our model (7) describes the evolution of wet areas, not heights, Q ± in the case of a steady state at rest. In order to obtain such discretization of the source terms, we begin by writing the numerical flux for Q 1 under these rest conditions, and seek a discretization of the cell average of the corresponding source term that matches this numerical flux S Q1 (x, W, W x ) = g w 2 ∂ x A 1 . A consistent discretization is given in the following proposition.
Proposition 7. Let us consider the following discretization of the source terms in (29) g ∆x g ∆x Then the numerical scheme given by equation (35) satisfies the well-balance property. That is, it recognizes steady states of rest.
The averages in the computations A k,j± 1 2 are done consistently with the discretization of the numerical fluxes in equation (36), which is relevant for the cases where w 2 ,ŵ 2 are constant.

Evolution
Once the interface values, the numerical fluxes and the average of the source terms have been calculated, the ODE system (35) is integrated in time using the second order Strong Stability Preserving Runge-Kutta scheme [Gottlieb et al., 2001], (56a) with the Runge-Kutta fluxes with S j (t) described by equation (29), and discretized by equation (55). The time step ∆t is determined so as to satisfy the CFL restriction where k = 1, 2 denote the corresponding quantities for the internal and external layers respectively. Here, τ e is given by and We note that in channels with straight walls, = 1 by construction. In the case of general channels, that quantity is only approximately 1, so the extra condition is not too restrictive. Furthermore, τ e is non-positive and such extra restriction is due to entrainment. We also note that we are considering entrainment only when the external layer A 2 is away from zero, and the above definition is valid. The quantity in (60) is a frequency scale given by the friction term. The condition in (58) guarantees that the friction is resolved in the sense that we include at least 10 time steps during each time period of length 1/τ f .
One can show that the CFL restriction in equation (58) also guarantees the positivity-preserving property.
That is, the positivity of A 1 and A 2 are preserved in the numerical evolution. The details are given in the following proposition.
Proposition 8. Consider the scheme (35) -(36) with the reconstruction algorithm described in §2.3 and the discretization of the source term given by proposition (7). If the cell averages A k (t), k = 1, 2 are such that then the cell averagesĀ(t + ∆t) as evolved with time evolution given by equation (56), under the CFL limitation (58), and with one-sided local speeds given by equation (38) yields A k,j (t + ∆t) ≥ 0 ∀j, k = 1, 2.
The proof can be done for each layer as in [Balbás and Hernandez-Duenas, 2014], and by adding the contribution of entrainment. It will not be repeated here.

Numerical Setup and Boundary Conditions
In this section we present numerical solutions for a range of problems illustrating the properties of our At the left (right) boundary, outflow occurs when γ − 1 < 0, γ − 2 < 0 (γ + 1 > 0, γ + 2 > 0) respectively. Inflow is considered to occur when such conditions are not satisfied. Unless otherwise noted, we extrapolate the variables w 1 , w 2 , Q 1 , Q 2 at outflow, and prescribe w 1 , w 2 , Q 1 , Q 2 at inflow (with the values given by the initial conditions).

Riemann Problem
System (7) consists of balance laws with non-conservative products. This adds both theoretical and numerical challenges such as in calculating physically relevant weak solutions. Non-conservative products may change the jump conditions in weak solutions. In [Dal Maso et al., 1995], a definition of weak solutions based on the theory of non-conservative products is provided and more on the theory of paths can be found in [LeFloch, 2002, LeFloch, 2004.
In this example we test the convergence and non oscillatory properties of the proposed scheme. We consider a Riemann problem consisting of a two-layer flow along a channel with constant width, σ(x, z) = 1, and a flat bottom B(x) = 0. Friction and entrainment are ignored so that n i = n b = 0, and V e = 0, and the We impose the boundary conditions as specified at the beginning of Section 3. The two-layers move originally at the same speed, u 1 (x, 0) = u 2 (x, 0) = 2.5 in both sides of the membrane, and the only difference is a jump on the location of the interface between the two-layers at t = 0 (h 1 increases and h 2 decreases by 0.05 while the flow remains flat at the top).
Although the geometry of the channel is trivial and does not pose any challenges for the numerical scheme, when the membrane is removed, the initial conditions develop strong shock discontinuities. If solved with non non-oscillatory numerical schemes or over coarse grids, these discontinuities will lead to the onset of spurious oscillations. The numerical results shown in Figure 4 demonstrate the robustness of our scheme and shows the convergence to the discontinuous solution as the grid is refined. The left panels show the internal layer w 1 (top left) and the external layer w 2 (bottom left). The corresponding velocities are shown in the right panels. In each panel, different resolutions are used. The red dotted line shows the numerical results using 1, 000 gridpoints, while the blue dashed line represents the solution computed with a finer grid of 2, 000 gridpoints. Due to the presence of non-conservative products, a Riemann problem is particularly challenging in two-layer flows. Although one can observe oscillations when a resolution of 1, 000 gridpoints is used, such oscillations are significantly reduced for the reference solution that is obtained with a very fine resolution of 10,000 gridpoints. The structure of the solution can be clearly identified. There are four shockwaves, one of them with negative speed of propagation.

Perturbation from a Steady State of Rest in a General Channel with Discontinuous Topography
A perturbation to a steady state of rest in a channel with general walls and discontinuous topography is applied in this section. The wall's width defined over the domain [0, 1] is given by The ratio of densities between layers is r = 0.98. Here V e = 0 but friction is included with Manning coefficients n i = n b = 0.009 s m −1/3 . Such values were obtained from [Khan and Lai, 2014]. The topography is given by The initial conditions are given by u 1 (x, 0) = u 2 (x, 0) = 0, w 1 (x, 0) = 0.7, and w 2 (x, 0) = The left panel of Figure 5 shows the 3D view of the channel described by equation (62)  The well-balance property is particularly important in two-layer flows. Any unbalance in the external layer affects the internal one and viceversa. Figure 6 (without the initial perturbation) shows the numerical results when the numerical scheme does not satisfy the well-balance property. In particular, here we use the same scheme but implement the reconstruction of A 1 and A 2 directly, instead of the procedure described in Waves due to gravity may form in the internal layer, generating internal waves in two-layer flows. Internal waves may appear in rivers, and shear flow instabilities may be one of the mechanisms leading to such waves [Apel et al., 1975]. Here we test the ability of the numerical scheme to capture steady states that are not at rest for the internal layer while the free surface is at rest. We consider a channel extended over the interval x ∈ [0, 2] with the topography and channel's width given, respectively, by and For the purpose of this numerical test, we ignore entrainment (V e = 0). We also ignore friction (n i = n b = 0) to verify that the conditions for steady states in Proposition (5) are met.
The discretization of the source terms (55) guarantees the balance between source terms and non linear fluxes for steady states of rest, but that balance is not guaranteed in steady states characterized by the conditions (22), where the equilibrium in the internal layer is achieved through the effects of hydrostatic pressure exerted by the external layer. Perturbations in any of the layers affect each other. The results obtained for this example demonstrate that the proposed numerical schemes computes correctly these non trivial steady states.
It is also worth noting that in these examples, despite the loss of hyperbolicity in a subregion, the approximations (21) that we employ to bound the eigenvalues allow us to capture the steady-state flows and resolve accurately their internal waves. Other steady states with large enough velocity u 1 , however, may lose hyperbolicity in the entire domain. In such cases, the proposed model is not suitable to describe those type of flows and should not be used at all. we specify w 1 = 0.9, u 1 = 0.3, w 2 = 1.5, u 2 = 0, and those same values are used for the initial conditions. We confirm the solution is a steady state with the plots of Q 1 , Q 2 , E 1 , E 2 in the bottom middle and bottom right panels. The discharges Q 1 , Q 2 and the energy E 1 are approximately constant. The energy E 2 seems to be converging to a piecewise constant profile. The topography, and the elevations of the internal and external layers are shown in the middle top panel, while the velocity is displayed in the top right panel. We clearly see that w 2 is constant, indicating that the external layer is at rest. The 3D view of this flow is shown in the left panel. In the middle top panel we also show with red asterisks the region where hyperbolicity is lost.
We notice that hyperbolicity is lost where u 1 is highest. This is consistent with the fact that hyperbolicity is lost when u 1 and u 2 are far from each other.
We now continue the simulation of the steady state reached in Figure 7 with a reconstruction that does not follow the procedure in sections 2.3 and 2.4 where w 1 and w 2 are reconstructed frist. Instead, both A 1 and A 2 are reconstructed directly with the non-oscillatory procedure in Section 2.3. As a consequence, this    . From top to bottom: the time evolution at times t = 0, 0.2, 1, 3, 10 is shown for w 1 and w 2 in the left column, and u 1 and u 2 in the right column respectively.
We further test the numerical scheme by adding a perturbation of δ w = 0.2 in the internal layer in the interval [0.1, 0.2]. The time evolution of this perturbation at times t = 0, 0.2, 1, 3, 10 is shown in Figure 9 from top to bottom for w 1 , w 2 (left column) and u 1 , u 2 (right column). We observe that the perturbation in the internal layer propagates throughout the domain and eventually leaves it, recovering its initial profile.
On the other hand, the external layer w 2 is always almost flat. The formulation in equation (7) helps maintaining the balance in this internal wave flow and computing the free surface correctly.  Figure 10: Positivity numerical test with initial conditions given by (66), topography given by (64) and wall's width by (65). Left column: 3D plot at times t = 0 and t = 1. Middle column: topography (black solid line), w 1 (blue solid line), and w 2 (red dashed line) at times t = 0, 0.25, 0.5, 50. Right column: the same as in the middle column for u 1 (blue solid line) and u 2 (red dashed line).

Lock Exchange
In this example we test the positivity preserving property of the numerical scheme presented in Section 2. The topography and wall's width are given by equations (64) and (65) from the previous numerical test.
The Manning coefficients for the interface and bottom are n i = n b = 0.009 s m −1/3 , and V e = 0.
The two-layers are initially at rest, u 1 (x, 0) = u 2 (x, 0) = 0, with the external layer occupying the left side of the channel and the internal layer the right side. That is, Here δ B = 10 −3 is used as a threshold. We note that this threshold is used to avoid loss of hyperbolicity.
The ratio of densities is r = 0.95. The boundary conditions used here are those specified at the beginning of Section 3. We expect the heavier fluid to push the external layer, resulting in a displacement of the internal layer below the external one. The time evolution at times t = 0, 0.25, 1, 50 are displayed in Figure 10. The elevations are shown in the middle column. As one can observe, preserving positivity is challenging and such property in the numerical scheme adds stability and accuracy to the approximated numerical solution. At t = 50, the flow has reached a steady state that is not at rest. The velocity in the right column indicates that the external layer moves to the right while the internal one goes in the opposite direction. The 3D view on the left column exhibits the complex geometry of the channel's wall and the fluid's evolution over time.  Figure 11: Time evolution of the solution with initial conditions given by equation (67). The topography, w 1 and w 2 are shown at times t = 0, 1, 1.5, 2. The initial conditions are specified in equation (67).

Currents
Gravity currents produced by lock exchanges have been studied in [Adduce et al., 2012]. Experimental data and numerical results are compared in [Adduce et al., 2012] using a shallow-water model with entrainment. It was shown that entrainment helps getting more realistic numerical results. In the numerical test considered in this section, we initiate with a fluid that is composed only of the external layer (lighter fluid). We then impose a discharge of heavier fluid in the internal layer. Specifically, we impose a discharge Q 1,left = 0.1, and w 1,left = 0.6, Q 2,left = 0, w 2,left = 1.5 at the left boundary. The right boundary is treated as specified at the beginning of Section 3.
Here, r = 0.95, g = 9.81 m 2 s −1 , and n i = n b = 0.009 s m −1/3 . Following [Adduce et al., 2012], the velocity entrainment is given by where G is the composite Froude number in equation (16), and k = 0.1. Furthermore, the Manning coefficients for the interface and bottom are n i = n b = 0.009 s m −1/3 respectively.
One of the challenges in this numerical test is that the internal layer has wet-dry states since the depth h 1 is initially small. The discharge on the left produces a wet-dry front propagating to the right, inducing a current in the internal layer. Figure 11 shows the topography (black dashed line), w 1 (solid blue line) and w 2 (red dashed line) at times t = 0, 1, 1.5, 2. One can clearly identify the wet-dry front in the internal layer. In addition to the above variables, the bottom right panel also includes the elevation of the internal layer computed without entrainment. Consistently, one can observe a higher elevation near the front for the computations implemented with entrainment.