The Momentum Conserving Scheme for Two-Layer Shallow Flows

: This paper confronts the numerical simulation of steady ﬂows of ﬂuid layers through channels of varying bed and width. The ﬂuid consists of two immiscible ﬂuid layers with constant density, and it is assumed to be of a one-dimensional shallow ﬂow. The governing equation is a coupled system of two-layer shallow water models. In this paper, we apply a direct extension of the momentum conserving scheme previously used for solving the one layer shallow water equations. Computations of various steady-state solutions are used to demonstrate the performance of the proposed numerical scheme. Under the inﬂuence of a given ﬂow rate, the numerical steady interface is generated in a channel topography with a hump. The results obtained conﬁrm the analytic steady interface of the two-layer rigid-lid model. Furthermore, the same scheme was used with an additional artiﬁcial damping to simulate the maximal exchange ﬂow in channels of varying width. The numerical steady interface agreed well with the analytical steady solutions.


Introduction
The presence of internal waves in Indonesian waters such as the Lombok Strait and the Andaman Sea inspired this research. The internal waves in the Lombok Strait form when a layer of warm low-salinity water from the Pacific Ocean enters the strait and flows into the Indian Ocean as part of the Indonesian throughflow (ITF); the complex bathymetry of the strait with a sill also contributes to internal wave generation in Lombok Strait. These internal waves, once generated, travel through the strait; after leaving the strait's mouth, the wave propagates radially across the ocean, to the north of the Flores Sea, and to the south of the Indian Ocean.
In order to construct a numerical model for simulating the generation and propagation of internal waves in natural straits, we must begin with a simpler numerical model that captures the most important behavior of exchange flows through channels connecting two basins with different hydrological characteristics, such as the Lombok Strait. Several researchers [1][2][3] and many others discuss the generation and propagation of internal waves in the Strait of Gibraltar. They adopted a one-dimensional two-layer shallow flow model, which can be applied to channels of varying topographies and widths, to model the influence of external flows in the strait. This one-dimensional approach is considered to be quite feasible for modelling the internal hydraulics of the strait and its effect on exchange flows [3][4][5]. While internal wave propagation in the radial direction across the ocean requires a higher dimensional approach (two-dimension or even three-dimension). Inspired by their research, in this article, we will focus on developing a numerical scheme for a one-dimensional two-layer shallow flow model as well as testing the ability of the numerical scheme to simulate the steady interface of the two-layer fluid flows.
The goal of this research study is to investigate the dynamics of flow through a one-dimensional channel connecting two reservoirs. Our method employs numerical simulation of a steady interface in a channel with varying bed and width. In this paper, we present a numerical model capable of simulating steady interfaces in a variety of situations. In contrast to the one-layer shallow water equation, the two-layer shallow water equation model is not always hyperbolic [1,[6][7][8][9]. This presents its own challenges for the proposed numerical scheme. So far, the majority of numerical schemes that have been developed are associated with the finite volume or the finite difference methods. A numerical scheme based on a Riemann solver is proposed and developed by [3,9]. A finite difference scheme for solving the two-layer rigid-lid model is developed in [1], whereas a phase resolving numerical model is described in [2,5]. References [2,5] developed a numerical model capable of describing dispersive waves. In this article, we will develop a numerical scheme for solving the two-layer shallow water equations. The scheme is an extension of the momentum conserving staggered-grid scheme (MCS) for the one-layer shallow water equations [10] (see also [11]). It turns out that the newly developed scheme can be applied directly to solve two-layer shallow water models. The scheme is then validated by simulating different steady two-layer flows, including those with hydraulic jumps. We will simulate the formation of a steady interface in the channel due to the presence of a bottom sill. The numerical scheme will be validated further by simulating the maximal exchange flow in a channel of varying width.
Thus, the outline of this paper is as follows: discussions about the two-layer shallow water models followed with some review about the steady hydraulics of two-layer rigid-lid models are presented in Section 2; formulations of the momentum conserving schemes for solving the two-layer shallow water models are presented in Section 3; and various steady-state problems are simulated in Sections 4 and 5 to validate the proposed scheme. Finally, the conclusion is provided in Section 6.

Mathematical Model
We begin with a discussion on the governing equations for the two-layer fluid in a one-dimensional channel. Here, the channel cross-section is assumed to be rectangular, with varying bed, and varying width. Assume that the fluid consists of two immiscible layers of different densities. The upper fluid has density ρ 1 , and the lower fluid has density ρ 2 , with ρ 2 > ρ 1 . Let z = η 1 (x, t) represent free surface deformation, z = η 2 (x, t) denotes the interface between two fluids, and z = d(x) denots the channel bottom topography, as shown in Figure 1. Under this setting, the thickness of the upper and lower layers is h 1 = η 1 − η 2 and h 2 = η 2 − d, respectively. We denote H 0 as the reference height, where H 0 = h 1 + h 2 + d. We employ a one-dimensional approach, which requires us to assume that the flow is uniform across the channel cross section. Moreover, the channel under consideration is symmetrical about the x-axis, with channel bed d(x) and width b(x) assumed to be slowly varying. The fluid particle velocities in the upper and lower layers is denoted as u 1 (x, t) and u 2 (x, t), respectively. Let A i (x, t) represent the cross-sectional area of fluid in each layer, whereas Q i (x, t) = A i (x, t)u i (x, t) denotes the discharge. We restrict our discussion to rectangular channels, so the cross-section area is sim- being the channel width. According to [7,12,13], the governing equations of the two-layer flows under hydrostatic pressure is given by the following.
∂A 2 ∂t Equations (1) and (3) are mass conservation for the upper and lower layer, respectively. Equations (2) and (4) are the momentum balances for each layer. In this paper, we will propose a numerical scheme for solving the two-layer model (1)-(4). However, we first provide a brief review here on the formulation of steady two-layer solutions.

Steady Two-Layer Solutions
Farmer and Armi [4], Armi [6] conducted an in-depth study of the steady-flow of the two-layer rigid-lid model. The results of their study have been widely used to validate numerical schemes [1,3,7,9]. In the following description, we will present a brief review of the Armi Farmer steady solution, which will be used to validate our numerical scheme.
In the steady situation, mass conservation (1) and (3) reduces toQ 1 = A 1 b and Q 2 = A 2 b, respectively, whereas the momentum Equations (2) and (4) lead to the specific energy for each layer as follows.
In the above description, notations with hats represent steady variables. By subtracting the energy equations in (5) and (6) after dividing with g ρ 2 H 0 and after some algebra, one obtains the following:Ê where r = ρ 1 ρ 2 and g ≡ g(1 − r), which represents reduced gravity. The energy difference (7) is also known as the Bernoulli equation [4].
Let the Froude number for each layer be denoted as follows.
As discussed in [6], the two-layer model (1)-(4) reaches a critical condition when G 2 = 1 in which G 2 the composite Froude number for two-layer flows with free surface is given by the following.
For layers with r ≈ 1, the composite Froude number can be approximated to G 2 ≈ F 2 1 + F 2 2 , and the criticality condition becomes the following.
If G 2 < 1, the flow is called subcritical, whereas if G 2 > 1 the flow is called supercritical. When G 2 = 1, the flow is said to be critical, and its location is called a control. Unlike the one-layer shallow water equation which is always hyperbolic, the two-layer Equations (1)-(4) are not always hyperbolic. The hyperbolic property of the two-layer equations depends on a parameter called the shear parameter, The model will be hyperbolic if S < 1, in this case, the characteristic velocities are real numbers. However, if S > 1 the characteristic velocities becomes complex numbers, and the model changes from hyperbolic to elliptic [14]. The stability condition (10) is related to the Kelvin-Helmholtz instability. If this instability occurs, this means that the two previously immiscible fluids begin to mix. The energy will be reduced as a result of this mixing process. Numerically, this instability will cause spurious error in the calculation, causing the scheme to become unstable [1,8,9,12]. To control this instability, some dissipation must be added to the model; in this study, the artificial damping terms ν∂ xx u 1 and ν∂ xx u 2 is added to the momentum Equations (2) and (4), respectively.
The most common approximation used is the rigid-lid approximation. In this approximation, free surface deformation is neglected or η 1 is constant over time. When dealing with a rigid lid, the following relation must be added: with H 0 the reference height.

Formulation Using Normalized Variables
Furthermore, we define the following non-dimensional variables of free surface, interface, fluid thicknesses, channel width, bottom elevation, and discharges: with b 0 denoting the reference width. The dimensionless variables are denoted with bars, whereas the physical variables do not have bars. From the Froude number formula (8), we have the following useful relationship.
By using the relation (13), we can rewrite the rigid-lid Equation (11) in terms of Froude number as follows.
Under the rigid-lid assumption, we have −h 2 −d =h 1 − 1, so the Bernoulli Equation (7) can be rewritten as follows.
If the channel has constant widthb = 1, the formula (16) can be written in terms of the Froude numbers and discharge ratioQ r = |Q 1 /Q 2 | as follows: where

∆Ê.
Similarly to the one-layer shallow water model where steady hydraulics can be described as a trajectory of the energy curve, the steady interface of the two-layer model (assuming existence is guaranteed) can be described as a trajectory on the energy difference curve (17). This aspect will be discussed in detail by using specific examples in Sections 4 and 5.

Numerical Model
In this section, we discuss the numerical scheme for solving the two-layer model (1)-(4). The scheme is based on the momentum conserving scheme for one-layer shallow water flows originally proposed by [15]. The scheme can be used to solve problems with rapidly varying flows, such as those involving hydraulic jumps and bores. This scheme has been implemented and extended to various shallow water flow problems, such as in [10,16,17] where it was referred to as the MCS scheme, an abbreviation for the momentum conserving staggered-grid scheme. The extension of the MCS scheme to the two-layer model has been successfully used for the study of internal waves in [18]. In this article, the scheme, which is limited to the hydrostatic model, will be used for discussion with a focus on simulating various steady flows in a channel with varying bed and width. The formulation of the MCS scheme for the two-layer hydrostatic model (1)-(4) will be described below.
Consider the spatial domain x ∈ [α, β] and time domain t ∈ [0, T] with α, β, T ∈ R. The spatial domain is discretized uniformly with a spatial step size ∆x/2 to yield partition points: domain is also discretized uniformly using the time step ∆t. In this staggered configuration, we approximate A i (x j , t n ) ≡ A n i,j at full grid points x j . On the other hand, we approximate u i (x j+1/2 , t n ) ≡ u n i,j+1/2 at half grid points x j+1/2 , see Figure 2. The semi discrete equations for (1)-(4) are read as follows. dA n It should be noted that the numerical damping terms with coefficient ν have been added to Equations (18) and (21). When the problem being simulated is not fully hyperbolic, as suggested by several authors [1,8,12], these additional damping terms are required to maintain the stability of the numerical scheme. This issue will be discussed in more detail in Section 5.
Using the analogy of the momentum conserving approximation of the shallow water equations for one-layer model as described in [10], we use the following approximation for the advection terms: and the first-order upwind approximation for horizontal velocity is the following.
Equations (18)-(21) are semi-discrete schemes. A fully discrete scheme that we used here is obtained by implementing forward time integration to (18)-(21), yielding the following.
Thus, the MCS-schemes for the two-layer shallow water equations are (27)-(30). This scheme is explicit and is based on the leapfrog method, which is discussed in [11] for the case of a one-layer shallow water flow. This staggered grid arrangement has the advantage of being free of numerical damping error. Various test cases will be simulated in Sections 4 and 5 to validate the proposed scheme. It should be noted that the artificial damping term is only used when the shear stability condition is violated, i.e., when S > 1.
For the linearized form of the two-layer model (1)-(4), there are two phase speeds: gH 0 corresponds to the fast mode, and g h 1 h 2 H 0 corresponds to the slow mode [18]. Thus, for all simulations performed here, ∆t is chosen so that the Courant number is gH 0 ∆t ∆x ≈ 0.5. Furthermore, a wet-dry procedure should be used to avoid numerical instability in simulations involving dry areas. In this case, a simple wet-dry procedure is used for each layer; for example, in the lower layer, u 2 is calculated by using (30) only when the corresponding cell is wet. While the cell is considered wet if A 2,j+1/2 > A thres , that is, when the wet cross-sectional area exceeds the predetermined value of A thres . The same is also holds for the upper layer.
Next, we will first show that the momentum conserving staggered grid scheme admits the well-balanced property. Well balanced here means that the scheme can maintain a steady state; if it was originally steady, it will stay steady. The following discussion is a generalization of the well-balanced property of the one-layer shallow water flow that has been derived in [19].
Since the channel width and depth does not vary in time, after dividing by b j = 0, for This means that η 1,j and η 2,j are constant at any time t ≥ t n . Furthermore, as η i,j , i = 1, 2 are both flat, Equations (19) and (21) become the following.
This means that u i,j+1/2 is constant at any time t ≥ t n . As we know u n i,j+1/2 = 0, then u i,j+1/2 is zero for any t ≥ t n . Hence, we have shown that surface and interface stay flat, and the two-fluid system remain motionless at any time t ≥ t n ; hence, the proof is complete.
Theorem 2 (Well balanced fully discrete). The fully discrete scheme (27)-(30) preserves the steady state of a lake at rest.
By using (22), the above equations can be expressed in terms η 1 and η 2 as follows.
Thus, for any time t ≥ t n , we obtain u i,j+1/2 = 0, and η i,j is constant for i = 1, 2. The proof is complete.

Steady State of a Lake at Rest
The first test is the simulation of steady state at rest. The spatial domain is [0, 12], with the spatial step size ∆x = 0.2. This simulation uses a channel with randomly generated topography d(x) and width b(x) as shown in Figure 3. The initial condition is u 1 (x, 0) = 0 and u 2 (x, 0) = 0 and η 1 (x, 0) = 1 and η 2 (x, 0) = 0.5, which represent flat surface and flat interface at rest. Both left and right boundaries are hard wall boundaries u i (0, t) = 0 and u i (12, t) = 0 for i = 1, 2. During the calculations it was observed that the free surface and interface remain flat and that velocities u 1 and u 2 are nearly zero (see Figure 4 (left)). As shown in Figure 4 (right), non-zero discharges Q 1 and Q 2 appear during computation, but they are still in the order of the round-off error; their values also do not increase over time. This test provides a computational evidence of the well-balanced property of numerical scheme (27)-(30) in the case of a channel with varying bed and width.

Simulation of Hyperbolic Case
In this section, we present the results of numerical calculations using the momentum conserving scheme (27)-(30) and simulate the formation of a steady interface in a two-layer flow model under the influence of bottom channel with a sill. All of the test cases presented here are simulated without the use of artificial damping, so we set ν = 0; additionally, the problems are hyperbolic with shear parameters less than one.

Smooth Subcritical Flow
This test case is taken from [7], focusing on simulating a steady interface in a channel with a bottom sill. The computational domain is x ∈ [−3, 3], with a spatial step size of ∆x = 0.05 and gravity acceleration of g = 10. The channel width is assumed to be one, and the bottom topography consists of a single hump represented by a smooth exponential function.
The initial condition is given by u i (x, 0) = 0, η 1 (x, 0) = 0 and η 2 (x, 0) = −0.5. The boundary conditions are the discharge Q 1 (−3, t) = 0.15 and Q 2 (−3, t) = −0.15 at the inflow, and the layer depths are h 1 (3, t) = 0.5 and h 2 (3, t) = 1.5 at the outflow. The steady situation of subcritical flow is depicted in Figure 5 (left). As shown in that figure, the steady interface forms a negative hump, while the free surface remains flat. The lower layer flow moves from left to right (u 2 positive), while the upper layer flow moves from right to left (u 1 negative). The discharge and energy curves of both layers are depicted in Figure 5 (right); their constant values indicate that the steady state has been reached. We also plot the energy difference and shear instability; the shear instability is less than one, indicating that the problem is still hyperbolic.

Smooth Transcritical Flow
In this section, we conducted simulations of transcritical flows using the two-layer model. This is the test case presented in [7], also known as the maximal two-layer flow of the rigid-lid model by [4]. To simulate the transcritical flow, we performed two simulations: one where both layers flow in the same direction (parallel flow) and the other where the two layers flow in opposite directions (exchange flow).

Parallel Flow
The computational domain is x ∈ [−3, 3], with spatial step size ∆x = 0.05, and gravity acceleration of g = 10. The channel width is taken to unity, whereas the bottom topography consists of a single hump represented by a smooth function.
The initial condition used is u i (x, 0) = 0, η 1 (x, 0) = 1.3338331, and η 2 (x, 0) = 1. In this case, the fluid in both layers are flowing in the same direction (parallel flow) by adopting the boundary conditions: positive discharges at the left inflow Q 1 (−3, t) = 0.09282893, Q 2 (−3, t) = 0.09282893, and fixed water thicknesses for both layers at the outflow are given by h 1 (3, t) = 1.3338331, where Q = 0.09282893 and h 2 (3, t) = 0.1616669. Starting from the initial condition, the numerical calculation shows the development of a steady interface, as shown in Figure 6 (left). That steady interface is of transcritical-type, which smoothly connects subcritical flow and supercritical flow on the left and right of the sill, respectively. As shown in that figure, the numerical steady interface agrees well with the analytical steady interface, with a slight deviation on the upstream part of the hump. When a steady state is reached, the discharge curves Q 1 and Q 2 and the energy curves E 1 and E 2 show constant values, as shown in Figure 6 (right). These flow parameters must indeed be constant for the steady analytical solution of the rigid-lid model. The fact that our simulation was able to produce a constant value for these parameters indicates that our two-layer MCS scheme can be used to simulate the rigid-lid case. In addition, we observed that the free surface is barely deformed during the simulation of this case. Furthermore, the shear parameter shown in Figure 6 (right) is less than one, indicating that the problem remains hyperbolic. We also ran another simulation for this smooth transcritical flow in which both layers flow in opposite directions with the same discharge. This is accomplished in the numerical simulation by using the boundary conditions as Q 1 (3, t) = −0.09282893, Q 2 (−3, t) = 0.09282893, and h 1 (−3, t) = 0.4311358 and h 2 (3, t) = 0.1616669. The numerical simulation yields a steady surface that is similar to the parallel flow case. This is to be expected because, as predicted by the rigid-lid analytic model, the steady interface is determined by the same values ofQ r and ∆Ê as in the parallel case.
Next, we will discuss the hydraulic aspect of this smooth transcritical steady flow. This description holds for both cases: parallel flow as well as exchange flow, with the discharge ratioQ r = |Q 1 /Q 2 | = 1. The channel topography (38) with constant channel widthb = 1 has the form of a smooth hump with a crestd = 0.33 (corresponds to the volume flow rate ofQ 2 (1−d) 3/2 = 0.21) and a footd = 0 (corresponds to the volume flow rate ofQ 2 (1−d) 3/2 = 0.11).
On the Froude number plane, we draw two rigid-lid contours: one corresponds to the foot, and another corresponds to the crest of the topographical hump. In numerical simulation, after reaching steady state ∆Ê = 1.5, we plot the Bernoulli curve with contour ∆Ê = 1.5. In Figure 7 (right), this Bernoulli curve intersects the rigid-lid contour = 0.11 at the point A; this condition represents the subcritical upstream part of the steady interface indicated by the point A. Moving forward, at the bottom crest, the Bernoulli and the rigid-lid contour = 0.21 intersect at B, exactly at a control point with G 2 = 1. Furthermore, the supercritical flow on the downstream part can be recognized by the point C. This downstream condition meets the fact that h 2 is small in the downstream area, implying that F 2 2 must be large.

Transcritical Flows with Jump
This is a simulation of steady interface in the case of transcritical flow with jump. Parameters used in this simulation are the same as the previous smooth transcritical case except that here we specify different boundary conditions. Discharges are specified at the inflow Q 1 (−3, t) = Q 2 (−3, t) = 0.09282893, and water heights for each layer are specified at the outflow h 1 (3, t) = 0.5794783 and h 2 (3, t) = 0.9205217. The resulting steady interface is depicted in Figure 8; when the flow proceeds over the sill, the subcritical flow accelerates to supercritical, which then undergoes a hydraulic jump to match the boundary condition at the right boundary. The numerical steady interface, as shown in Figure 8 (left), is in good agreement with the exact steady-state solution, and some discrepancy appears at the hydraulic jump location. At the end of the simulation time, the discharge curve in Figure 8 (right) shows a constant value, but the energy curve, especially E 2 , is not really constant. This indicates that our model is only an approximation of the rigid-lid model. Furthermore, the shear instability curve in Figure 8 (right) is less than one, suggesting that the problem remains hyperbolic.

Simulation of Conditionally Hyperbolic Cases
In this second subsection, we simulate steady interface in the case of channel with varying width. We conducted two test cases, both of which turned out to be conditionally hyperbolic case. Thus, in this case, we used a numerical scheme (27)-(30) that includes the artificial damping terms that is three-two orders of magnitude lower than the other terms.

Maximal Exchange Flows through a Contraction
This test case is known as the maximal exchange flow problem, which was introduced by [1]. On the computational domain x ∈ [−1.5, 1, 5], with spatial step size ∆x = 0.01, gravity acceleration g = 9.81, and damping parameter ν = 0.002. This calculation makes use of a flat bed channel with d(x) = 0 and varying width given by the following: with α = 1.073. The function (39) represents a normal channel width of b = 5 with a single contraction b = 1 at the origin x = 0 and b = 2 at x = ±0.5. Here, we use h 1 (x, 0) = 0.5, h 2 (x, 0) = 0.5, and u 1 (x, 0) = −u 2 (x, 0) = 0.2 as the initial conditions. In this simulation, the upper layer is flowing to the right, while the lower layer is flowing to the left. The boundary conditions u 1 (−1.5, t) = 0.2 and u 2 (1.5, t) = −0.2 are used in the numerical simulation, along with absorbing boundaries for the other two, u 1 (1.5, t) and u 2 (−1.5, t). Note that, for the absorbing boundary, we used the sponge layer technique here by [20]; for the computational interval −1.5 < x < 1.5, we place at the far right x = 1.5 a sponge layer for u 1 on 1.5 < x < 2 and at the far left x = −1.5 a sponge layer for u 2 on −2 < x < −1.5. When the steady state is reached, the interface forms the maximal exchange flow, as illustrated in Figure 9 (left).  The numerical steady interface shows a good agreement with the analytical steady interface of the rigid-lid model. We observe that the free surface remains flat during the simulation, indicating that our two-layer numerical scheme can be used to simulate the rigid-lid case. As shown in Figure 9 (right) both the discharges Q 1 and Q 2 and the energies E 1 and E 2 are constant, indicating that stability has been achieved. The G 2 curve shows that the maximal exchange flow is supercritical on both the upstream and downstream sides of the channel contraction. The control location with G 2 = 1 in this case is the channel contraction x = 0. In addition, this simulation requires the inclusion of a damping term; otherwise, instability will occur. Indeed, hyperbolicity is not guaranteed in this problem because the shear parameter curve is nearly one throughout the domain, as shown in Figure 9 (right).
In the following, we will recognize the steady interface on the Froude number plane. Following [6] and using (5) and (6), we first rewrite the energy difference formula for the zero channel bed d = 0 as follows.
The later expression is obtained after using (8). The given boundary conditions result in the upper layer flowing to the right, whereas the lower layer flowing to the left at the same rate yields Q r = | Q 1 Q 2 | = 1. In this case, the rigid-lid Equation (15) and the Bernoulli equation are read as follows.
Our interest is on the maximal two-way exchange flow between two infinite reservoirs; with density ρ 1 on the left and ρ 2 on the right. In this case, the inflowing layer spreads out to the point where the fluid thickness approaches zero. Since Q i is finite, the vanishing h i implies a Froude number F i that tends to infinity for i = 1, 2. The point A in Figure 10 (left) represents a position near reservoir ρ 1 , with F 2 → ∞, whereas the point C represents a position near reservoir ρ 2 , with F 1 → ∞. Consider Figure 10 (right), among several Bernoulli curves, there is only one level curve, i.e., with ∆E = 0.5 that satisfies these boundary conditions; that curve is plotted as a thick solid curve. This Bernoulli curve intersects the criticality condition (9) at a point (F 2 1 = 1 2 , F 2 2 = 1 2 ), i.e., the point B, located at the narrowest section, which becomes a control with G 2 = 1. Furthermore, the rigid-lid curve (40) with contourQ 2 b = 1 4 passes through the point B. Thus, the solid curve A − B − C in the Froude number plane represents the steady interface of the maximal exchange flow. Moreover, this steady flow is supercritical on both sides of the control.

Exchange Flows through Sill and Contraction
The final test case is taken from [1]. The steady interface in this case represents the exchange flow through a channel with contraction and a bottom sill. The computational domain is x ∈ [−1, 2] with the spatial step size ∆x = 0.01 and gravitational acceleration of g = 9.81, and the damping parameter is ν = 0.01. The channel beds and widths are given by the following: where α = 0.673 for x ≤ 1, α = 1.273 for x ≥ 1, and β = 3.75.
For this simulation, the initial conditions are h 1 (x, 0) = 1, h 2 (x, 0) = 1, u 1 (x, 0) = 0.2, and u 2 (x, 0) = −0.18. To achieve an exchange flow with the upper layer flowing to the right and the lower layer flowing to the left, we use boundary conditions u 1 (−1, t) = 0.2 and u 2 (2, t) = −0.18, with absorbing boundary conditions for the upper layer downstream and the lower layer upstream. In this case, the absorbing boundary is built using sponge layer technique similar to those described in Section 5.1. Figure 11 (left) depicts the steady interface and the composite Froude number G 2 . The numerical steady interface shows a good agreement with the analytical steady interface of the rigid-lid model. We observe that the free surface remains flat during the simulation, indicating that our two-layer numerical scheme can be used to simulate the rigid-lid case. Figure 11 (right) shows that both the discharges Q 1 and Q 2 , as well as the energies E 1 and E 2 , are constant, indicating that steady state has been achieved. The G 2 curve shows that the maximal exchange flow is supercritical on both the upstream and downstream sides of the channel, with a subcritical area trapped between the two controls, i.e., at the bottom sill x = 0 and at the channel contraction x = 1. The shear parameter curve shown in Figure 11 (right) is greater than one on the upstream side of the sill, therefore, hyperbolicity is not guaranteed for this problem. To avoid numerical instability, the numerical damping term must be incorporated Figure 11. (Left) Steady interface in the exchange flow through channel with contraction and a bottom sill. (Right) Plots of discharge and energy for each layer; the shear parameter (10) is greater than one, indicating a loss of hyperbolicity in that region.

Conclusions
We have proposed an extension of the momentum conserving staggered-grid scheme (MCS) for the two-layer shallow water flows in channels with varying bed and width. Formulation of the scheme is a direct analog of the MCS scheme for the one-layer shallow water equations. The scheme is explicit and well-balanced, and it has been tested with the steady-atrest case. The scheme was validated using various test cases of steady interface in a channel with bottom sill, including subcritical flow, smooth transcritical flow, and transcritical flow with hydraulic jump. The obtained results agree well with the exact steady interface of the two-layer rigid-lid model. In the simulation of the maximal exchange flow, the numerical scheme succeeded in producing a steady interface. In this simulation, an artificial numerical damping was used to avoid numerical instability. This was necessary because the problem is classified as conditionally hyperbolic. We concluded that the MCS scheme performs well in simulating various steady flows in channels with varying bed and width. However, the scheme developed here is limited to rectangular cross-section channels.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript or in the decision to publish the results.