https://orcid.org/0000-0002-8659-3753
Abstract
This paper proposes a novel energy-balanced numerical scheme for the two-layer shallow water equations (2LSWEs) that accurately captures internal hydraulic jumps without introducing spurious oscillations. The proposed scheme overcomes the problem of post-shock oscillations in the 2LSWE, a phenomenon commonly observed in numerical solutions of non-linear hyperbolic systems when shock-capturing schemes are used. The approach involves reconstructing the internal momentum equation of 2LSWEs using the correct Hugoniot curve via a set of shock wave fixes originally developed for single-layer shallow water equations. The scheme successfully preserves all stationary solutions, making it highly suitable for simulations of real-life scenarios involving small perturbations of these conditions.
HIGHLIGHTS
The paper proposes a remedy to the problem of post-shock oscillations in the two-layer shallow water equations.
It is an extension of a family of well-balanced and shock-resolving schemes that were originally developed for single-layer shallow water equations.
The resulting scheme, called HLL-MSF, successfully preserves all stationary solutions, including those with shock discontinuities.
INTRODUCTION
The two-layer shallow water equations (2LSWEs) are a set of mathematical equations in fluid dynamics that are widely used to model the behavior of fluids in a two-layer system. These equations are particularly useful for modeling geophysical flows in stratified conditions, which are frequently observed in coastal areas such as estuaries and sea straits (Schijf & Schönfled 1953; Ungarish 2013; Majd & Sanders 2014; Fernández-Nieto et al. 2016). An analytical solution for these equations cannot be found due to their non-linear and complex nature. As a result, numerical methods such as finite difference or finite volume are often used to solve these equations.
The numerical simulation of 2LSWEs faces several challenges, with the loss of hyperbolicity being the primary difficulty. The loss of hyperbolicity occurs when the internal eigenvalues of the system become complex. This has significant implications, as numerical methods that are designed to solve hyperbolic problems may not be applicable to elliptic or mixed-type systems. Therefore, specialized numerical techniques are required to accurately solve the 2LSWEs in the presence of hyperbolicity loss. To overcome this difficulty, various approaches have been proposed in the literature. These include ignoring the imaginary parts of complex eigenvalues, locally adding additional shear stress to restore the hyperbolic character (Castro-Díaz et al. 2010), introducing a third layer to avoid the loss of hyperbolicity and restore the stability (Bryson et al. 2010; Chertock et al. 2013) and using a splitting technique by solving each layer individually as a fully de-coupled system and then controlling the stability by an entropy condition (Bouchut & Morales de Luna 2008) or using a semi-coupled form by excluding some of the coupling terms (due to momentum exchange between the layer) from the coupled flux terms that may contribute to the loss of hyperbolicity (Swartenbroekx et al. 2010).
The 2LSWES contains non-conservative products arising from the layer coupling and geometric bottom. The existence of non-conservative products also presents a significant difficulty as these terms are not well-defined within the framework of distribution. The theory introduced by Dal Maso et al. (1995) tackles this problem by interpreting the non-conservative products as a bounded measure based on a family of paths. This interpretation enables the weak formulation of the equations and leads to the development of a new category of numerical schemes called ‘path-conservative’. The success of such schemes depends on the choice of the correct path, which should be based on the underlying physics of the problem. The presence of non-conservative products also allows for some non-trivial steady-state solutions. Preserving these steady states through numerical methods is essential since many real-life simulations of 2LSWEs are small perturbations of these conditions. Several techniques have been developed to maintain the steady state or the so-called well-balanced property, including pre-balanced formulations by modifying the governing equations to ensure the preservation of quiescent equilibrium (Rogers et al. 2003; Li et al. 2014), a path-conservative scheme which involves adding an extra equation related to temporal changes in bed elevation, hydrostatic relaxation approach (Audusse et al. 2004; Berthon & Foucher 2012; Gunawan & Lhébrard 2015), and the upwind discretization of the source terms by balancing the source terms according to the upwind fluxes (Bermudez & Vázquez 1994).
Some examples of the unwind approach include the lateralized HLL scheme (Soares-Frazão & Zech 2011; Meurice & Soares-Frazão 2020) developed for single-layer shallow water equations with fixed and movable bed and the augmented Riemann solvers such as ARoe or HLLS derived for single- and two-layer shallow water equations (George 2008; Murillo & García-Navarro 2010; Murillo et al. 2020). However, numerical schemes like Roe that rely on the explicit Jackson matrix can be computationally expensive for complex systems such as the two-layer shallow water equations where the analytical expressions of their eigenstructure are not readily available. The pre-balanced and relaxation approaches are only suitable for preserving the steady state at rest or ‘lake-at-rest’ situation, whereas the augmented solvers have proven effective in preserving not only the steady state at rest but also those containing contact waves with non-zero velocity (moving water equilibrium). Since the system is hyperbolic, stationary discontinuities in the form of shock waves may develop. Despite recent progress in using upwind approaches to resolve the stationary shock wave (the hydraulic jump) in single-layer shallow water equations (Navas-Montilla & Murillo 2019; Akbari & Pirzadeh 2022), accurately capturing the internal hydraulic jump in a two-layer shallow water system remains a challenge (Murillo et al. 2020). The numerical solution of the hydraulic jump is often associated with smears and oscillations, even for high-resolution schemes like Roe (Castro-Díaz et al. 2010).
The spurious oscillations near the jump are due to the numerical slowly moving shock anomaly, which can occur when shock moves slowly with respect to the mesh. This phenomenon was first observed by Woodward & Colella (1984) in the Euler equations when using the finite volume Godunov scheme. The consensus among researchers is that the presence of intermediate states due to the shock-capturing approach and the non-linearity of the jump conditions are the underlying causes of the numerical slowly moving shock anomaly (Arora & Roe 1997; Zaide & Roe 2012). The shock-capturing numerical schemes aim to capture the shock with a finite width, but this width can often be larger than the physical width of the shock. This results in the introduction of intermediate states that have no direct physical interpretation (Zaide & Roe 2012). These intermediate states can lead to spurious oscillations and errors in the numerical solution, which can compromise the accuracy of the numerical results. There are several techniques that can be used to mitigate or eliminate the numerical slowly moving shock anomaly, including adding artificial viscosity or numerical stabilization, and introducing shock-fitting or grid-aligning techniques (Jin & Liu 1996; Paciorri & Bonfiglioli 2009).
Earlier studies focused mainly on conservative systems like the Euler equations of gas dynamics. However, in recent years, researchers have developed treatments for the non-conservative systems as well, specifically the shallow water equations. Navas-Montilla & Murillo (2017) successfully extended the Zaide & Roe (2012) flux fixes, initially developed for the conservative Euler equations, to the single-layer shallow water equations, by taking into account the effect of the non-conservative products arising from the geometric bottom. These flux fixes involve replacing the original non-linear shock in the intermediate cell with a linear version derived from neighboring cells. The resulting scheme ensures a uniform mass-flow rate in smooth bottom cases without a ‘spike’ across the hydraulic jump. The technique also significantly reduces the post-shock oscillations caused by the slowly moving shock anomaly. However, this approach does not ensure the exact solution of the water depth inside the cell containing the stationary shock (hydraulic jump), and its ability is restricted to the numerical flux of the Roe scheme. Later, Akbari & Pirzadeh (2022) introduced a set of slowly moving and stationary shockwave fixes which are based on reconstructing momentum flux and its non-conservative product using the correct Hugoniot curve that connects the left and right limits of the jump. By reproducing the correct strength and speed of the jump inside the momentum equation, they are able to remove the non-physical components that cause the numerical artifacts. This novel technique offers several advantages over existing methods. It effectively reduces the post-shock oscillations caused by the slowly moving shock anomaly while also ensuring accurate solutions for water depth and mass-flow rates across the stationary shock wave (hydraulic jump). Additionally, this approach is compatible with less computationally demanding numerical schemes, such as Lax-Friedrich and HLL.
Although there has been some research on the behavior of numerical shock anomalies in the single-layer shallow water equations, to the best of our knowledge, there exists no published research on the numerical slowly moving shock anomaly in the context of the two-layer shallow water equations. In this paper, we propose a solution to this problem by extending the slowly moving and steady shock wave fixes introduced by Akbari & Pirzadeh (2022) to the two-layer shallow water equations. This involves deriving an energy-balanced HLL-type scheme for the 2LSWES, which is then treated with appropriate shock wave fixes to overcome the numerical shock wave anomaly. The resulting scheme, named HLL-MSF, is shown to capture internal hydraulic jumps, ideally by two cells, without producing smears or oscillations in the numerical solution. The numerical results also confirm convergence to the correct weak solution of the stationary shock. The proposed numerical scheme is suitable for simulation of nearly stationary flows offering an efficient alternative to the more expensive Roe scheme, which requires treatment for the loss of hyperbolicity and non-negativity of water layers.
The structure of the paper is as follows: In Section 2, the mathematical model of the two-layer shallow water equations in one space dimension is derived. The finite volume methods used to solve the conservation laws with source terms are introduced in Section 3. In Section 4, the problem of numerical slowly moving shock anomaly in the two-layer shallow water system is analyzed and then the proposed shock wave fixes are presented. The paper concludes with some numerical experiments to evaluate the performance of the proposed scheme.
GOVERNING EQUATIONS
, the fluxes 
, and the sources 
are defined by: 
and 
denote the non-negative layer depth, and the mass-flow rate of j-th layer, respectively, with 
represents the upper layer and 
represents the lower layer, 
denotes the density ratio, x refers to the spatial variable, t is the time, z is the variation due to bottom, and 
is the gravitational constant. The flow is also inviscid and no friction is assumed between the layers and the bottom. We can rewrite system (1) in the following quasi-linear form: 
is now the vector of source terms only due to the geometric bottom and 
is the Jacobian matrix of the flux vector coupled with the non-conservative products describing the momentum exchange between the layers, which can be defined by: 
and 
also are the horizontal velocity and the celebrity, respectively. The eigenvalues of matrix 
, denoted by 
, are the roots of the characteristic equation, namely: 
is specified as 
. The flow condition is said to be critical if one of the eigenvalues in (1) vanishes. The critical condition is equivalent to the following expression for the composite Froud number G (Abgrall & Karni 2009): 
is the internal Froude number and 
is the reduced gravity. The regime of the flow is, then, critical if 
, subcritical if 
and supercritical if 
. The system admits a family of steady states that satisfies the following ODEs: 
is the specific energy for the layers. System (1) also admits non-smooth ‘steady-state’ solutions that require a definition for the weak solution. For hyperbolic conservation laws, one can use the classical notion of weak solution in the distribution framework and derive the Rankine–Hugoniot jump conditions. However, such a notion cannot be used for hyperbolic systems in non-conservative form. The jump conditions for the two-layer SWEs are often derived under the rigid-lid assumption (Abgrall & Karni 2009). This assumption is considered reasonable when the density difference between the layers is small. Even after defining a weak solution for (1) under a certain assumption, the admissibility of such a weak solution should be checked with an entropy-like condition, as the weak solutions of conservation laws may not be unique.NUMERICAL METHODS
into control volumes 
whose lengths are 
. The finite volume method applied to (1) over this control volume can be written as: 
is the cell averaged value at a certain time step 
defined by: 
, one requires a suitable numerical flux, often a first-order Riemann solver, and an approximation to the integral of source terms. Several techniques have been developed for the treatment of non-conservative products (LeFloch 1989; Vázquez-Cendón 1999; LeFloch & Thanh 2011). In general, to preserve all steady states, the following discrete equality must hold: 
, one can rewrite (13) to: 
is defined by 
, and 
contains both the flux differences and the integral of non-conservative source terms. This coupled formulation gives rise to the so-called augmented numerical fluxes of the form: 
The numerical flux
and fastest 
wave speeds. So, let us define the two decomposed jumps by 
and 
which moves with the speeds 
and 
, respectively. By applying the Rankine–Hugoniot conditions across these two jumps, we have: 
and 
are defined, respectively, by: 
in (17) using the HLL intermediate state 
(Harten et al. 1983). The intermediate state and the jump conditions (18) described by the original scheme of the HLL are derived under the notion of weak solution, which only applies to the hyperbolic conservation laws in conservative form. To generalize the HLL scheme for hyperbolic systems with source terms such as (1), a similar technique offered by Akbari & Pirzadeh (2022) is adopted. This method can be considered as a modified version of the lateralized HLL scheme first developed by Soares-Frazão & Zech (2011) for the single-layer shallow water equations with source terms. Based on the approach used in Akbari & Pirzadeh (2022), to construct well-balanced numerical fluxes for the non-conservative system (1), the strength of the jumps (18) must be adjusted by: 
is now the intermediate state of the HLL for the non-conservative system which is computed by: 
is a vector that controls the jump strengths, such as: 
can be defined as the ratio between the determinant of uncoupled and coupled flux Jacobian. The second and fourth components of 
are not associated with any non-conservative term, so, their values are: 
in the numerator is the determinant of the coupled Jacobian associated with the jump strength 
which can be defined by: 
is the determinant of the uncoupled flux Jacobian given by: 
denotes the source terms corresponding to the momentum equations of the j-th layer. Notice that at critical points, when the determinant of the flux Jacobian is close to zero, the correct entropy may not be produced using (23). So, the following fixes are prescribed to treat the numerical viscosity terms in case of a transonic rarefaction: 
Other alternatives are available; see, for instance, Murillo et al. (2020). For most cases, (29) is enough to ensure the stability of the numerical solutions under the CFL condition (see Abgrall & Karni (2009) and Kurganov & Petrova (2009) for further discussions on the choice of wave bounds for the 2LSWEs).
Construction of an energy-balanced source term
is the integral of the source term that verifies the properties introduced in (10) and 
denotes as: 
and 
are the cell averaged source terms computed by the upper and lower bounds of the cell interval 
, respectively, defined by: 
and 
(31), one can use a weighted combination, namely: 
for each layer as: 
to 
at the hydraulic jump (the linear approximations usually provide correct dissipation across the discontinuity (Navas-Montilla & Murillo 2019)), and assuming 
, the weight coefficient (32) can now be defined by: 
Similarly, 
can be defined for negative flow
. It is worth noting that the procedures used here for the construction of the well-balanced numerical fluxes and source terms can be easily extended to other numerical methods, such as FORCE (Stecca et al. 2010; Martínez-Aranda et al. 2020) and central-upwind (Kurganov et al. 2001; Kurganov & Petrova 2009) schemes which are frequently used in the literature for the numerical approximation of the two-layer SWEs. In what follows, we first examine the ability of the proposed numerical scheme to preserve the smooth steady states of still and moving water, and then, we assess the failure of such a regular finite volume method to accurately capture the internal hydraulic jump. In all test cases we present here, the CFL number is set to 0.4.
Test 1) Quiescent equilibrium (lake-at-rest):
The simulation is conducted over smooth (case A) and discontinuous (case B) bottom configurations, which are respectively defined by:
- Case A) Smooth bed:
- Case B) Discontinuous bed:
- and 
-norm as follows: 
is the exact solution of a conserved variable and 
is the corresponding numerical solution for N compactional cells at the final time of simulation. As evident by 
and 
errors in Table 1, the proposed numerical scheme is exactly well-balanced for the steady state at rest, verifying the conservation property introduced in (9). Test 2) Trans-critical steady state without internal bore (shock):
Test 1 case A: The numerical solution of the water layers and the corresponding mass-flow rates at rest over a bump with 50 mesh cells.
Test 1 case A: The numerical solution of the water layers and the corresponding mass-flow rates at rest over a bump with 50 mesh cells.
Test 1 Case B: The numerical solution of the water layers and the corresponding mass-flow rates at rest over a discontinuous bed with 50 mesh cells.
Test 1 Case B: The numerical solution of the water layers and the corresponding mass-flow rates at rest over a discontinuous bed with 50 mesh cells.
and 
errors for Test 1
| Case |  |  |  |  | ||||
|---|---|---|---|---|---|---|---|---|
 |  |  |  |  |  |  |  | |
| A | 2.2204 × 10−16 | 1.1102 × 10−16 | 2.2653 × 10−12 | 4.6133 × 10−14 | 2.7756 × 10−15 | 1.1102 × 10−16 | 4.3155 × 10−12 | 4.6133 × 10−14 |
| B | 2.2204 × 10−16 | 1.1102 × 10−16 | 9.5406 × 10−12 | 1.9134 × 10−13 | 1.1102 × 10−15 | 1.1102 × 10−16 | 2.5442 × 10−11 | 5.0966 × 10−13 |
| Case | | | | | ||||
|---|---|---|---|---|---|---|---|---|
| | | | | | | | |
| A | 2.2204 × 10−16 | 1.1102 × 10−16 | 2.2653 × 10−12 | 4.6133 × 10−14 | 2.7756 × 10−15 | 1.1102 × 10−16 | 4.3155 × 10−12 | 4.6133 × 10−14 |
| B | 2.2204 × 10−16 | 1.1102 × 10−16 | 9.5406 × 10−12 | 1.9134 × 10−13 | 1.1102 × 10−15 | 1.1102 × 10−16 | 2.5442 × 10−11 | 5.0966 × 10−13 |
is given by: 
and the density ratio is taken as 0.98. To create a smooth trans-critical flow, a discharge rate 
is imposed on the upstream while at the downstream, the interface is set to 0.1 until 
when a free boundary condition is applied. The simulation is ended at 
. We compute the numerical solution using five uniform meshes 50, 100, 200, 400, and 800. Since no analytic solution can be defined, a reference solution is computed by a very fine mesh of 3,200 cells to validate numerical results. For a better comparison, the obtained reference solution is only used for the layer depths while a constant 
is considered as the exact solution of the mass-flow rates. The results of this comparison are reported in Table 2. The numerical solution computed with 100 meshes is also depicted in Figure 3. In Figure 4, we also provide the numerical results of the specific energy computed by the different mesh grids. The proposed numerical scheme, as observed in Figures 3 and 4, ensures a uniform solution for the specific mass-flow rates and energy in the entire computational domain, thus verifying the condition in (10). Test 3) Moving water steady state with internal bore (jump):
Test 2: The numerical solution of layer depths (left) and the mass-flow rates (right) using 100 meshes.
Test 2: The numerical solution of layer depths (left) and the mass-flow rates (right) using 100 meshes.
Test 2: The numerical solution of the specific energy for the lower (right) and upper (left) layers using different meshes.
Test 2: The numerical solution of the specific energy for the lower (right) and upper (left) layers using different meshes.
and 
errors for Test 2 using different meshes
| Cell |  |  |  |  | ||||
|---|---|---|---|---|---|---|---|---|
 |  |  |  |  |  |  |  | |
| 50 | 0.0015 | 5.8621 × 10−05 | 6.2794 × 10−14 | 1.8700 × 10−15 | 0.0368 | 7.6416 × 10−04 | 5.9022 × 10−14 | 2.1927 × 10−15 |
| 100 | 0.0063 | 1.1781 × 10−04 | 7.9395 × 10−14 | 2.1094 × 10−15 | 0.1467 | 0.0015 | 2.5327 × 10−13 | 4.8399 × 10−15 |
| 200 | 0.0135 | 2.9607 × 10−04 | 4.6779 × 10−13 | 5.7107 × 10−15 | 0.3100 | 0.0018 | 1.2363 × 10−12 | 1.3010 × 10−14 |
| 400 | 0.0050 | 1.1833 × 10−04 | 4.7766 × 10−12 | 2.6108 × 10−14 | 0.1146 | 3.9240 × 10−04 | 2.3205 × 10−11 | 1.1305 × 10−13 |
| 800 | 0.0057 | 6.0391 × 10−05 | 1.0448 × 10−11 | 2.8668 × 10−14 | 0.1331 | 1.7281 × 10−04 | 6.3716 × 10−12 | 3.3834 × 10−14 |
| Cell | | | | | ||||
|---|---|---|---|---|---|---|---|---|
| | | | | | | | |
| 50 | 0.0015 | 5.8621 × 10−05 | 6.2794 × 10−14 | 1.8700 × 10−15 | 0.0368 | 7.6416 × 10−04 | 5.9022 × 10−14 | 2.1927 × 10−15 |
| 100 | 0.0063 | 1.1781 × 10−04 | 7.9395 × 10−14 | 2.1094 × 10−15 | 0.1467 | 0.0015 | 2.5327 × 10−13 | 4.8399 × 10−15 |
| 200 | 0.0135 | 2.9607 × 10−04 | 4.6779 × 10−13 | 5.7107 × 10−15 | 0.3100 | 0.0018 | 1.2363 × 10−12 | 1.3010 × 10−14 |
| 400 | 0.0050 | 1.1833 × 10−04 | 4.7766 × 10−12 | 2.6108 × 10−14 | 0.1146 | 3.9240 × 10−04 | 2.3205 × 10−11 | 1.1305 × 10−13 |
| 800 | 0.0057 | 6.0391 × 10−05 | 1.0448 × 10−11 | 2.8668 × 10−14 | 0.1331 | 1.7281 × 10−04 | 6.3716 × 10−12 | 3.3834 × 10−14 |
Test 3: The numerical solution of the interface (left) and the mass-flow rate of the lower layer (right) using different meshes.
Test 3: The numerical solution of the interface (left) and the mass-flow rate of the lower layer (right) using different meshes.
Test 3: The numerical solution of the specific energy for the lower (left) and upper (right) layers using different meshes.
Test 3: The numerical solution of the specific energy for the lower (left) and upper (right) layers using different meshes.
INTERNAL BORE
In this section, we study the jump structure for the non-linear system of balanced laws (1). The Rankine–Hugoniot jump conditions for some limiting cases are derived, and the admissibility of the jump conditions is also discussed. We analyze the numerical slowly moving shock anomaly for these limiting cases. Finally, treatments are introduced to cure the post-shock oscillations caused by the numerical shock wave anomaly.
Rankine–Hugoniot conditions
connected by a set of left 
and right 
states as (see Figure 7): 
Schematic diagram of the internal bore moving with the speed 
.
Since the system is non-conservative, the standard notation of weak solution cannot be used to describe the jump conditions. To drive Rankine–Hugoniot jump conditions for the two-layer shallow water equations, we may use the rigid-lid assumption offered by Abgrall & Karni (2009). As there are no unique jump conditions, two limiting cases are studied, 
and 
, in which the jump conditions can either be described in the classical sense (or the general one due to the presence of the non-conservative bottom, see LeFloch & Thanh (2007)) or can be constructed under the rigid-lid assumption suggested by Abgrall & Karni (2009).
RH for
:
as an additional condition, where 
denotes a jump of any quantity. Assuming 
, the generalized version of the Rankine–Hugoniot conditions applied across the jump discontinuity can be described by: 
, the speed of jump vanishes and (38) gives: 
Equations (39), (40), and (42) show the non-linear relation between the two conserved variables. The non-linearity of jump conditions is often cited as the source cause of numerical shockwave anomaly observed in the solution of slowly moving shocks.
RH for
(rigid-lid):
, the Rankine–Hugoniot conditions applied across the jump can be defined by: 
: 
Equation (44) is similar to the one given by (42), except the effect of gravity pressure is reduced due to stratification.
Entropy satisfying jump
(and 
as soon as an explicit expression for the characteristic values is available). The jump discontinuity with the speed 
is entropy satisfying in the sense of Lax if it satisfies: 
) becomes analogous to: 
Similar conditions can be defined for 
. The system is also endowed with an entropy pair associated with the total energy, requiring the total energy to dissipate across the shock (see Abgrall & Karni (2009) for more detail on this).
Discrete jump and ideal internal bore
, 
, and 
. From (39)–(44), depending on the absence or the presence of the non-conservative terms, the numerical solution of the jump may lie on the following curves, respectively: 
and 
denote the right and left limits of the jump. The discrete jump is captured by the right and left limits and a number of intermediate states. A more dissipative scheme typically uses a larger number of intermediate states to resolve shocks, as the dissipation spreads out the shock wave over a greater number of cells. The exact position of the shock wave within these cells is commonly determined by the equal area rule. This rule states that the area above the pre-shock solution must be equal to the area below the post-shock, meaning that the total mass or momentum contained in the region above the pre-shock solution and below the post-shock solution must be conserved. From the mass conservation in the lower layer, the position of shock can be determined by: 
to avoid ambiguity in the exact shock position and also 
and 
to ensure the entropy stability. For a non-linear hyperbolic system, this is not, however, the case. In fact, the hydraulic jumps that satisfy such a condition must be ideal, that is, the following property hold for them: 
So, an ideal internal hydraulic jump must be connected with two states, the one that also satisfies the entropy inequality (either the Lax entropy condition or the entropy pair associated with the total energy). To identify the underlying problem of the post-shock oscillations, we examine the slowly moving and steady shock for the case of a zero-density ratio using the following numerical experiments.
Test 4) Slowly moving jumps for
over a flat-bed:
The test is taken from Akbari & Pirzadeh (2022) and consists of two cases of slowly right- and left-moving shocks for the single-layer SWEs. The current test aims to analyze the behavior of slowly moving shocks for the limiting case 
. For this purpose, the following initial conditions are taken.
- Case A) A right-moving shock:
and,
- Case B) A left-moving shock:
and 
are imposed. In Figure 8, the Hugoniot locus, together with the numerical solutions of the left- and right-moving shocks at 
, are presented. As we see from Figure 8, the numerical solutions of the slowly moving shocks lie on the Hugoniot curve shown with the solid line verifying the mass and momentum conservation. However, as observed in Figure 8, the intermediate states generate a set of shock curves that are not the ones connecting the right and left limits of the jump. In Figure 9, we plot the space-time evolution of the mass-flow rates. As evident in Figures 8 and 9, the numerical spike is responsible for the spurious waves generated downstream of the jump. 
is also used to denote the position of the intermediate state (middle state) containing the shock where the characteristic changes sign. The straight line that connects the right and left limits of the jump, depicted with the dashed line in Figure 8, is the slope (41) and can be computed by: 
Test 4: The Hugoniot locus compared with the numerical solution of slowly left- (left) and right-moving (right) shocks at 
.
Test 4: The Hugoniot locus compared with the numerical solution of slowly left- (left) and right-moving (right) shocks at .
Test 4: The space-time evolution of the mass-flow rates for the right - (right) and left-moving (left) jumps.
Test 4: The space-time evolution of the mass-flow rates for the right - (right) and left-moving (left) jumps.
Test 5: The Hugoniot locus compared with the numerical solution of a single-layer hydraulic jump.
Test 5: The Hugoniot locus compared with the numerical solution of a single-layer hydraulic jump.
The numerical solution is monotonicity preserving if the intermediate states lie on this slope. We see in Figure 8 that the intermediate states are above this line for both left- and right-moving shocks resulting in a numerical overshoot which is always ahead of the shock.
Test 5) Steady jump for
over a bed topography:
, where the exact solution can be explicitly defined. The test has the initial conditions defined by: 
denotes the cell that contains the shock. However, as expected, the discrete jump is smeared into more than one intermediate state. The closest cells to the exact left and right states in the Hugoniot locus are 
and 
, respectively. By specifying 
and 
as the left and right states, the following can be defined for the discrete jump condition: 
is the right state. To avoid imaginary solution, 
is only defined for 
. The discrete and exact shock curves overlap at the limit 
. The numerical solution of the jump lies on the discrete shock curve satisfying the discrete mass and momentum conservation. The numerical solution, however, is not monotonicity preserving. Similar to cases in Test 4, the mass and momentum fluxes must be constant 
across the shock in order to avoid the numerical artifacts.Numerical treatment of post-shock oscillations
Slowly moving shock wave fixes
, 
, and 
in the internal momentum equation by 
, 
, and 
when the shock is right moving and 
, 
, and 
when the shock is left moving, respectively. The new physical fluxes using moving shock wave fixes are given by: 
is reconstructed using 
, 
, and 
instead of 
, 
, and 
in (33), respectively, where, 
denotes as 
, and 
is also defined by:
is used to locate the cell (or cells) that contains the internal shock with 
is calculated using 
, and 
computes the direction of the internal shock containing in the cell i, namely: 
and 
is defined by 
and 
, respectively. Equation (56) is analogous to the sign 
(51) for the limiting case 
. After constructing the new physical fluxes and the non-conservative source terms in the lower layer with the proposed moving shock wave fixes, we can now define the fourth component of 
as: 
Stationary shock wave fixes
is the source terms (33) reconstructed using 
, 
, and 
instead of 
, 
, and 
, where, the notations 
and 
represent 
and 
, respectively, the superscript SF denotes the stationary shock wave fixes, and 
and 
are the one-sided fixes given by,
The moving and stationary shock wave fixes
is the mass-flow rate at the shock front defined by: 
is the mass-flow rate for 1-shock which can be approximated by: 
. To reconstruct the numerical flux differences (28), the flux limiters should be modified according to the new momentum equation defined for the lower layer; the first, second, and third components of 
are therefore redefined as: 
, and 
and 
are given by: 
, and 
are computed by: 
and 
can be calculated as: 
The same can be defined for 
and 
. Notice that the flux limiters 
should be bounded between 0 and 1. The proposed weighted technique allows the numerical scheme to single out the physically relevant shock from all shock curves involved in the Hugoniot locus. The HLL scheme using the moving and stationary shock wave fixes is hereafter called the HLL-MSF scheme.
Test 6) Revisiting Test 3: Trans-critical flow with internal hydraulic jump:
are depicted in Figures 11 and 12, respectively. As evident by the figures, the numerical solutions obtained by the proposed HLL-MSF scheme exhibit no smears or spurious oscillations near the discontinuity. The internal hydraulic jump is also captured ideally by two cells in all meshes. The proposed HLL-MSF ensures a uniform solution for the mass-flow rate in the entire computational domain. In Figure 13, it can be seen that the numerical solution of the specific energy is also constant away from the internal hydraulic jump, where energy is dissipated. Table 3 provides the positions of the left and right states connecting the internal hydraulic jump, as computed by the HLL-MSF scheme using different mesh sizes and density ratios. As indicated by Table 3 and Figure 12, in all cases, the numerical scheme accurately captures the internal hydraulic jump location. Compared to the reference solutions, the numerical scheme underestimates both the interface and mass-flow rates. However, the error due to these differences seems to remain within the numerical width of the shock (as the shock positions are accurately captured in all cases regardless). As seen in Figure 11, the HLL-MF scheme also gives similar results as HLL-MSF but remains oscillatory at the shock positions, indicating that the moving shock wave fixes, unlike the HLL-MSF scheme, are unable to construct the internal stationary shock wave accurately. To assess the scheme's performance with respect to the CFL number, simulations were also conducted using CFL values of 0.55 and 0.75. Irrespective of the CFL number chosen, the HLL-MSF method accurately captures the internal hydraulic jump without any smearing or oscillations, further highlighting the scheme's stability and robustness. Test 7) Perturbations of steady state with internal bore:
Test 6: The numerical solution of the interface (left) and the mass-flow rates in the lower layer (right) computed by the HLL-MF scheme for r = 0.98 using different meshes.
Test 6: The numerical solution of the interface (left) and the mass-flow rates in the lower layer (right) computed by the HLL-MF scheme for r = 0.98 using different meshes.
Test 6: The numerical solution of the interface (left) and the mass-flow rates in the lower layer (right) computed by the HLL-MSF scheme for r = 0.98 using different meshes.
Test 6: The numerical solution of the interface (left) and the mass-flow rates in the lower layer (right) computed by the HLL-MSF scheme for r = 0.98 using different meshes.
Test 6: The numerical solution of the specific energy for the upper (left) and lower (right) layers computed by the HLL-MSF scheme for r = 0.98 using different meshes.
Test 6: The numerical solution of the specific energy for the upper (left) and lower (right) layers computed by the HLL-MSF scheme for r = 0.98 using different meshes.
The position of internal jump for the reference solution and the proposed HLL-MSF scheme using different meshes and density ratios
| Position of internal bore | ||||||||
|---|---|---|---|---|---|---|---|---|
 |  |  |  | |||||
| xl | xr | xl | xr | xl | xr | xl | xr | |
| Ref. | 13.24375 | 13.25 | 13.29375 | 13.3 | 13.45 | 13.4563 | 13.6 | 13.6062 |
| 50 | 13.2 | 13.6 | 13.2 | 13.6 | 13.2 | 13.6 | 13.6 | 14 |
| 100 | 13.2 | 13.4 | 13.2 | 13.4 | 13.4 | 13.6 | 13.6 | 13.8 |
| 200 | 13.2 | 13.3 | 13.2 | 13.3 | 13.4 | 13.5 | 13.6 | 13.7 |
| 400 | 13.2 | 13.25 | 13.25 | 13.3 | 13.45 | 13.5 | 13.6 | 13.65 |
| 800 | 13.225 | 13.25 | 13.275 | 13.3 | 13.45 | 13.475 | 13.6 | 13.625 |
| Position of internal bore | ||||||||
|---|---|---|---|---|---|---|---|---|
| | | | |||||
| xl | xr | xl | xr | xl | xr | xl | xr | |
| Ref. | 13.24375 | 13.25 | 13.29375 | 13.3 | 13.45 | 13.4563 | 13.6 | 13.6062 |
| 50 | 13.2 | 13.6 | 13.2 | 13.6 | 13.2 | 13.6 | 13.6 | 14 |
| 100 | 13.2 | 13.4 | 13.2 | 13.4 | 13.4 | 13.6 | 13.6 | 13.8 |
| 200 | 13.2 | 13.3 | 13.2 | 13.3 | 13.4 | 13.5 | 13.6 | 13.7 |
| 400 | 13.2 | 13.25 | 13.25 | 13.3 | 13.45 | 13.5 | 13.6 | 13.65 |
| 800 | 13.225 | 13.25 | 13.275 | 13.3 | 13.45 | 13.475 | 13.6 | 13.625 |
in order to create an internal hydraulic jump. Once the solution is steady an internal perturbation is generated downstream of the jump as follows: 
is the internal perturbation which is set to 
. The perturbation is adequately strong with respect to the grid size (numerical width of the shock) in order to disturb the position of the internal hydraulic jump. At the downstream, free boundary conditions are imposed. We simulate the results using two uniform meshes, 200 and 800. The evolution of the internal wave at different times is depicted in Figure 14. The perturbation downstream caused the internal hydraulic jump to travel upstream. The new position of the stationary internal jump for both coarser and finer meshes is the same limit verifying the conservation property of the scheme. Based on the numerical results, we can conclude that the proposed numerical schemes can handle small perturbations of the steady state containing internal shock. 
Test 7: The numerical solution of the interface at time t = 0 (upper left), t = 5 (upper right), t = 50 (lower left), and t = 500 (lower right).
Test 7: The numerical solution of the interface at time t = 0 (upper left), t = 5 (upper right), t = 50 (lower left), and t = 500 (lower right).
CONCLUSION
This paper extends a family of exactly well-balanced numerical schemes, recently developed for the single-layer shallow water equations with variable bottom topography, to the two-layer model. The objective of this extension is to preserve all stationary solutions of 2LSWEs, including those with shock discontinuities (specifically the internal hydraulic jump). The paper starts by introducing an HLL-type energy-balanced scheme for the two-layer shallow water equations. However, it is shown that while this energy-balanced scheme can preserve smooth steady states, it fails to accurately capture the internal hydraulic jump, resulting in smears and oscillations that cannot be remedied by refining the mesh. The anomalous behavior of numerical schemes near the shock is caused by the non-linearity of jump conditions and the presence of intermediate states to resolve the shock at the mesh grid. To overcome this problem, a set of slowly moving and stationary shock wave fixes, recently developed for the single-layer SWEs, are extended to 2LSWEs. These proposed shock wave fixes aim to eliminate the spurious waves that are not associated with the physical shocks. In particular, the stationary shockwave fixes are introduced to accurately maintain the stationary internal hydraulic jump, while the slowly moving shockwave fixes are designed to dampen the unphysical oscillations produced by the left and right slowly moving shocks. To cure the numerical anomalies pertaining to both slowly moving and stationary shock waves, the internal momentum equation is reconstructed using a weighted combination of these fixes, where the weight coefficient is determined based on the jump conditions. The proposed approach is evaluated through several numerical experiments involving bottom topography, various density ratios, and flow conditions. The proposed scheme has been shown to capture the internal hydraulic jump ideally without introducing smears or spurious oscillations. The numerical results evidenced that the proposed scheme effectively preserves both smooth and non-smooth internal steady-state solutions of the two-layer shallow water equations. Future research directions include extending the approach to other balance law systems, particularly focusing on multilayer shallow water equations with the inclusion of Manning friction and varying cross-sections, adapting the method to higher-order numerical schemes, and exploring its applicability to two-dimensional systems.
DATA AVAILABILITY STATEMENT
All relevant data are available from an online repository or repositories: https://github.com/Shallowflows/Two-layer-SWEs.
CONFLICT OF INTEREST
The authors declare there is no conflict.
%20cropped.png?versionId=7142)