## 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

*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: where 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: where and also are the horizontal velocity and the celebrity, respectively. The eigenvalues of matrix , denoted by , are the roots of the characteristic equation, namely:

*et al.*(2018) for the recent analytic solution of eigenvalues and Schijf & Schönfled (1953) for the approximate one). The Jacobian matrix (4) has four, two external and two internal, eigenvalues corresponding to barotropic and baroclinic components of the flow, respectively. The internal eigenvalues also appear to be complex under certain conditions (Schijf & Schönfled 1953; Abgrall & Karni 2009), indicating that the system is conditionally hyperbolic. To overcome the loss of hyperbolicity, Swartenbroekx

*et al.*(2010) suggest rewriting system (3) in a semi-coupled form where the non-conservative terms that contribute to the loss of hyperbolicity are accounted by the source terms rather than the coupled flux terms. Their proposed system admits four real and distinct eigenvalues, which can be defined by: where the subscripts upp and low denote the upper and lower layers, respectively, and 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): where 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:

## NUMERICAL METHODS

### The numerical flux

*et al.*(1983) is based on dividing the Reimann problem into three states separated by the slowest 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: where and are defined, respectively, by:

*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: where is now the intermediate state of the HLL for the non-conservative system which is computed by: and is a vector that controls the jump strengths, such as:

*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

*A*is the determinant of the Jacobian matrix given as:

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:

*t*= 5. The results for cases A and B are depicted in Figures 1 and 2, respectively, and the errors are reported in Table 1. These errors are computed using - and -norm as follows: where 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):

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} |

**Test 3)**Moving water steady state with internal bore (jump):

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} |

## 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

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**:

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):**

Equation (44) is similar to the one given by (42), except the effect of gravity pressure is reduced due to stratification.

### Entropy satisfying jump

*p*denotes

*p*-characteristic. For a steady shock where the characteristic changes sign across the discontinuity, the Lax entropy condition (assuming ) 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

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.

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:

^{2}/s is imposed on the upstream, and at the downstream, the interface is set to 0.63 m. In Figure 10, the curve shown with a solid line is the exact solution defined by (42). The curve shown with the dashed line is called the discrete jump condition (48), which is an imaginary curve that satisfies the mass and momentum conservation at the discrete level. To define the discrete jump conditions, we need to specify the position of the left and right states. In Figure 10, 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: where 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

*.*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: where the superscript

*MF*denotes the moving shock wave fixes. The numerical HLL scheme that uses moving shock wave fixes is hereafter called the HLL-MF scheme.

#### Stationary shock wave fixes

*SF*denotes the stationary shock wave fixes, and and are the one-sided fixes given by,

#### The moving and stationary shock wave fixes

*D*is a linear weight given by: and

*M*is the discrete jump condition as follows: where is the mass-flow rate at the shock front defined by: and is the mass-flow rate for 1-shock which can be approximated by:

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:

**Test 7)**Perturbations of steady state with internal bore:

. | 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 |

## 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.

## REFERENCES

*Communications in Mathematical Sciences*

**5**(4),