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.

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

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.

The system of two-layer shallow water equations in one space dimension can be written in the following differential form:
(1)
where the vectors of the conserved variables, , the fluxes , and the sources are defined by:
(2)
where 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:
(3)
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:
(4)
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:
(5)
The quartic Equation (5) does not provide an easy explicit expression for the eigenvalues. In the absence of a simple analytic solution, the eigenvalues of (1) can be computed approximately (see Krvavica 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:
(6)
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):
(7)
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:
(8)
The simplest of them is the steady state at rest described by the zero mass-flow rate as follows:
(9)
In general, the smooth solutions of (8) are expressed by the conservation of mass and energy in each layer, namely:
(10)
where 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.
In this section, we briefly introduce the finite volume method regularly used for the numerical solution of the balanced laws. We begin by partitioning the space-time into control volumes whose lengths are . The finite volume method applied to (1) over this control volume can be written as:
(11)
where is the cell averaged value at a certain time step defined by:
(12)
To update the numerical solution , 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:
(13)
which is the direct integral to (8). By setting the integral of source terms as , one can rewrite (13) to:
(14)
where the notation 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:
(15)
Plug (15) into (11) results in:
(16)

The numerical flux

In the present work, an HLL-type upwind scheme is derived for the two-layer SWEs. The method defined here does not need an explicit computation of eigenvalues; thus, no additional treatment is required for the loss of hyperbolicity. The scheme of Harten 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:
(17)
where and are defined, respectively, by:
(18)
One can compute the numerical fluxes 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:
(19)
where is now the intermediate state of the HLL for the non-conservative system which is computed by:
(20)
and is a vector that controls the jump strengths, such as:
(21)
The vector 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:
(22)
To be more precise, the ratio between the determinant of the uncoupled and coupled flux Jacobian for the second and fourth components is equal to one. Therefore, no change for the jump strengths corresponding to these two is required. For the first and third components, one can use the following expressions:
(23)
The term in the numerator is the determinant of the coupled Jacobian associated with the jump strength which can be defined by:
(24)
(25)
and the denominator is the determinant of the uncoupled flux Jacobian given by:
(26)
where 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:
(27)
The numerical fluxes of the HLL scheme can now be defined from (20) as:
(28)
Following the suggestion by Toro (2009), the minimum and maximum wave bounds can be given also as:
(29)

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

To guarantee the well-balanced property, the integral of source terms must exactly balance the flux gradients, meaning that the mass and momentum conservation must hold at the discrete level. To verify (10) for all smooth steady states, the numerical solutions must be also energy conservative. In a smooth steady state, the discrete momentum and energy equations can be written by:
(30a)
(30b)
where is the integral of the source term that verifies the properties introduced in (10) and denotes as:
(31)
where and are the cell averaged source terms computed by the upper and lower bounds of the cell interval , respectively, defined by:
(32)
Following Murillo & García-Navarro (2013) to formulate an energy-balanced scheme, instead of a linear combination between and (31), one can use a weighted combination, namely:
(33)
and then from system (1), obtain the weight coefficients for each layer as:
(34)
The energy-balanced source terms (33) ensure energy conservation for smooth solutions. However, across the jump, energy must be dissipated. By setting 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:
(35)
where A is the determinant of the Jacobian matrix given as:
(36)

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 first test case is the simulation of a steady state at rest or the so-called lake-at-rest. The example is taken from Abgrall & Karni (2009). The problem has the following lake-at-rest initial conditions in a domain with the interval [−1, 1]:

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:

The density ratio is set to 0.98. The domain is discretized into 50 uniform computational cells. The numerical solutions are generated at the time 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):

Figure 1

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.

Figure 1

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.

Close modal
Figure 2

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.

Figure 2

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.

Close modal
Table 1

and errors for Test 1

Case



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



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 
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 
The current test case is intended to verify the ability of the proposed numerical scheme to preserve the moving water's steady states without an internal hydraulic jump. The bed function is given by:
The initial conditions are defined by:
The channel interval is 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):

Figure 3

Test 2: The numerical solution of layer depths (left) and the mass-flow rates (right) using 100 meshes.

Figure 3

Test 2: The numerical solution of layer depths (left) and the mass-flow rates (right) using 100 meshes.

Close modal
Figure 4

Test 2: The numerical solution of the specific energy for the lower (right) and upper (left) layers using different meshes.

Figure 4

Test 2: The numerical solution of the specific energy for the lower (right) and upper (left) layers using different meshes.

Close modal
Table 2

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 
In this test case, we assess the failure of the traditional finite volume method to accurately capture the internal stationary jump. The bed function is the same as Test 2 over the same interval [0 20]. The initial condition is an internal dam defined by:
We consider free boundary conditions at both ends. The density ratio is set to 0.98. The numerical simulation is run until the steady state is reached. The steady state should contain an internal hydraulic jump at the downstream of the bed. The numerical solution is computed using four uniform meshes 50, 10, 200, and 400. The computed results are shown in Figures 5 and 6. No exact or reference solution can be considered here. Nevertheless, a constant ‘uniform’ mass-flow rate and an energy loss are expected across the jump discontinuity, whereas in the absence of the jump, a smooth solution must be ensured. However, as observed in Figures 5 and 6, the numerical method produces a smeared shock that spans over a number of intermediate cells far away from the jump. The spurious oscillations are due to the numerical slowly moving shock anomaly. The numerical anomaly observed across the jump is a common phenomenon of the hyperbolic systems of conservation laws that have non-linear jump conditions. It is described by a start-up error (a numerical spike) followed by spurious waves that ruin the downstream condition. Since the jump conditions are non-linear unless the intermediate state overlaps with the left or right states, the numerical spike cannot be resolved using traditional finite volume methods. It is worth noting that the smeared shock is due to the excessive amount of numerical diffusion produced by the HLL scheme, which is more than the correct entropy required. The scale of these shedding waves is proportional to the size of mesh. For instance, in the case of the coarser mesh with 50 or 100 cells, no internal hydraulic jump is formed due to this excessive diffusion (the internal jump moves outside of the boundary). On the other hand, the smears seem to be less visible on the fine meshes. Nevertheless, the numerical spike appears even for a very fine mesh.
Figure 5

Test 3: The numerical solution of the interface (left) and the mass-flow rate of the lower layer (right) using different meshes.

Figure 5

Test 3: The numerical solution of the interface (left) and the mass-flow rate of the lower layer (right) using different meshes.

Close modal
Figure 6

Test 3: The numerical solution of the specific energy for the lower (left) and upper (right) layers using different meshes.

Figure 6

Test 3: The numerical solution of the specific energy for the lower (left) and upper (right) layers using different meshes.

Close modal

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

Let us define a jump discontinuity with the speed connected by a set of left and right states as (see Figure 7):
(37)
Figure 7

Schematic diagram of the internal bore moving with the speed .

Figure 7

Schematic diagram of the internal bore moving with the speed .

Close modal

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:

The lower layer will behave as a single-layer system when the density ratio is zero. By dropping the coupling terms, the only non-conservative terms that remain are due to the geometric bottom. Since the system is non-conservative, the classical notion of weak solutions does not apply. To give a sense to the notion of weak solutions for the single-layer shallow system with a non-flat bottom, LeFloch & Thanh (2007) suggested using 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:
(38a)
(38b)
(38c)
Depending on the third condition (38c), two possible scenarios can occur. If the bottom is discontinuous, that is, , the speed of jump vanishes and (38) gives:
(39)
Considering that the jump condition (39) contains an additional term arising from the singular bottom, one requires to define a correct path to connect the right and left states of the jump. The jump condition (39) is, therefore, path-dependent. In the case where the bottom is flat, the system (38) reduces to the standard RH relations, namely:
(40a)
(40b)
From the first condition (41a), the speed of jump is:
(41)
Substituting (41) into (40b) gives:
(42)

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 rigid-lid assumption is considered reasonable when the density ratio between the layers is close to one. Following Abgrall & Karni (2009), the rigid-lid assumption can be used as an extra condition to describe the internal hydraulic jump. Assuming , the Rankine–Hugoniot conditions applied across the jump can be defined by:
(43a)
(43b)
(43c)
In (43), only the conditions for the lower layer are provided with the third condition due to the rigid-lid assumption. From the above system, the following can be computed for :
(44)

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

Entropy satisfying jump

To check the admissibility of the jump conditions, we follow Lax (1973) and define the entropy satisfying jump at least for the limiting case (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:
(45)
where the superscript p denotes p-characteristic. For a steady shock where the characteristic changes sign across the discontinuity, the Lax entropy condition (assuming ) becomes analogous to:
(46)

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

Let us assume a discrete 1-jump located inside the cell m using , , 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:
(47)
or
(48)
where 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:
(49)
The jump location given by (49) should be consistent with the other conserved variables namely 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:
(50)

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:

For both test cases, 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:
(51)
Figure 8

Test 4: The Hugoniot locus compared with the numerical solution of slowly left- (left) and right-moving (right) shocks at .

Figure 8

Test 4: The Hugoniot locus compared with the numerical solution of slowly left- (left) and right-moving (right) shocks at .

Close modal
Figure 9

Test 4: The space-time evolution of the mass-flow rates for the right - (right) and left-moving (left) jumps.

Figure 9

Test 4: The space-time evolution of the mass-flow rates for the right - (right) and left-moving (left) jumps.

Close modal
Figure 10

Test 5: The Hugoniot locus compared with the numerical solution of a single-layer hydraulic jump.

Figure 10

Test 5: The Hugoniot locus compared with the numerical solution of a single-layer hydraulic jump.

Close modal

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:

The test aims to analyze the numerical solution of a steady jump generated over a varying topography for the limiting case , where the exact solution can be explicitly defined. The test has the initial conditions defined by:
The bed function is the same as in Test 2. A mass-flow rate of 0.6 m2/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

In general, when the shock is mobile, we follow the technique proposed by Akbari & Pirzadeh (2022) for the single-layer SWEs to correct the internal momentum flux of two-layer SWEs and its corresponding non-conservative products. Since the spurious waves are ahead of the shock, the strength of upwind running waves must vanish to avoid the post-shock oscillations. So according to the direction of the shock wave, we must replace , , 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:
(52)
and the energy-balanced source terms in the lower layer is reconstructed using , , and instead of , , and in (33), respectively, where, denotes as , and is also defined by:
If
(53)
else
(54)
where the first condition, i.e., 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:
(55)
By expanding (55) using a similar projection of the source term as (25), where the linear approximation of the source term is used for the lower layer and the energy-balanced formulation for the upper layer, we can have:
(56)
where the notations 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:
(57)
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

A constant momentum flux must be ensured across the steady shock. Following Akbari & Pirzadeh (2022), the fix for the stationary shock is given as follows:
(58)
where 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,
If
(59)
else
(60)

The moving and stationary shock wave fixes

Finally, to have a robust algorithm that handles both slowly moving and steady internal jumps, a weighted technique based on the discrete jump conditions defined (48) is used. The proposed weighted technique is as follows:
(61)
where D is a linear weight given by:
(62)
and M is the discrete jump condition as follows:
(63)
where is the mass-flow rate at the shock front defined by:
(64)
and is the mass-flow rate for 1-shock which can be approximated by:
(65)
We can similarly compute . 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:
(66)
where , and and are given by:
(67)
and
(68)
where D is the weight coefficient (62), , and are computed by:
(69)
and and can be calculated as:
(70)
and
(71)

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:

In this test case, we evaluate the performance of the proposed HLL-MF and HLL-MSF schemes to simulate the internal hydraulic jump over a non-flat bottom. For this purpose, we revisit Test 3. This time however, we use four different density ratios 0.99, 0.98, 092, and 0.8. We compute numerical solutions with the proposed HLL-MF and HLL-MSF scheme using also five different meshes 50, 100, 200, and 800. The corresponding reference solutions are also computed by the proposed HLL-MSF scheme using a very fine mesh with 3,200 cells. The numerical results compared with the reference solution using the HLL-MF and HLL-MSF schemes for 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:

Figure 11

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.

Figure 11

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.

Close modal
Figure 12

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.

Figure 12

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.

Close modal
Figure 13

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.

Figure 13

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.

Close modal
Table 3

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




xlxrxlxrxlxrxlxr
Ref.13.2437513.2513.2937513.313.4513.456313.613.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




xlxrxlxrxlxrxlxr
Ref.13.2437513.2513.2937513.313.4513.456313.613.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 
The current test is devoted to the simulation of a moving internal bore and is intended to verify the ability of the proposed HLL-MSF scheme to capture small perturbations of the steady states with the internal hydraulic jump. To introduce a perturbation to the steady state, we revisit Test 2, and the same conditions are imposed except this time the interface at the downstream boundary is set to in order to create an internal hydraulic jump. Once the solution is steady an internal perturbation is generated downstream of the jump as follows:
where 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.
Figure 14

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

Figure 14

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

Close modal

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.

All relevant data are available from an online repository or repositories: https://github.com/Shallowflows/Two-layer-SWEs.

The authors declare there is no conflict.

Abgrall
R.
&
Karni
S.
2009
Two-layer shallow water system: A relaxation approach
.
SIAM Journal on Scientific Computing
31
(
3
),
1603
1627
.
https://doi.org/10.1137/06067167x
.
Akbari
M.
&
Pirzadeh
B.
2022
Implementation of exactly well-balanced numerical schemes in the event of shockwaves: A 1D approach for the shallow water equations
.
International Journal for Numerical Methods in Fluids
94
(
7
),
849
895
.
https://doi.org/10.1002/fld.5076
.
Arora
M.
&
Roe
P. L.
1997
On postshock oscillations due to shock capturing schemes in unsteady flows
.
Journal of Computational Physics
130
(
1
),
25
40
.
Audusse
E.
,
Bouchut
F.
,
Bristeau
M.
,
Klein
R.
&
Perthame
B.
2004
A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows
.
SIAM Journal on Scientific Computing
25
(
6
),
2050
2065
.
https://doi.org/10.1137/S1064827503431090
.
Bermudez
A.
&
Vázquez
M. E.
1994
Upwind methods for hyperbolic conservation laws with source terms
.
Computers & Fluids
23
(
8
),
1049
1071
.
https://doi.org/10.1016/0045-7930(94)90004-3
.
Berthon
C.
&
Foucher
F.
2012
Efficient well-balanced hydrostatic upwind schemes for shallow-water equations
.
Journal of Computational Physics
231
(
15
),
4993
5015
.
https://doi.org/10.1016/j.jcp.2012.02.031
.
Bouchut
F.
&
Morales de Luna
T.
2008
An entropy satisfying scheme for two-layer shallow water equations with uncoupled treatment
.
ESAIM: Mathematical Modelling and Numerical Analysis
42
(
4
),
683
698
.
https://doi.org/10.1051/m2an:2008019
.
Bryson
S.
,
Epshteyn
Y.
,
Kurganov
A.
&
Petrova
G.
2010
Well-balanced positivity preserving central-upwind scheme on triangular grids for the Saint-Venant system
.
ESAIM: Mathematical Modelling and Numerical Analysis
45
(
3
),
423
446
.
https://doi.org/10.1051/m2an/2010060
.
Castro-Díaz
M. J.
,
Fernández-Nieto
E. D.
,
González-Vida
J. M.
&
Parés-Madroñal
C.
2010
Numerical treatment of the loss of hyperbolicity of the two-layer shallow-water system
.
Journal of Scientific Computing
48
(
1–3
),
16
40
.
https://doi.org/10.1007/s10915-010-9427-5
.
Chertock
A.
,
Kurganov
A.
,
Kurganov
A.
,
Qu
Z.
&
Wu
T.
2013
Three-layer approximation of two-layer shallow water equations
.
Mathematical Modelling and Analysis
18
(
5
),
675
693
.
https://doi.org/10.3846/13926292.2013.869269
.
Dal Maso
G.
,
LeFloch
P.
&
Murat
F.
1995
Definition and weak stability of nonconservative products
.
Journal de Mathematiques Pures et Appliquees
74
,
483
548
.
Fernández-Nieto
E. D.
,
Garres-Díaz
J.
,
Mangeney
A.
&
Narbona-Reina
G.
2016
A multilayer shallow model for dry granular flows with the μ(I)-rheology: Application to granular collapse on erodible beds
.
Journal of Fluid Mechanics
798
,
643
681
.
https://doi.org/10.1017/jfm.2016.333
.
George
D. L.
2008
Augmented Riemann solvers for the shallow water equations over variable topography with steady states and inundation
.
Journal of Computational Physics
227
(
6
),
3089
3113
.
https://doi.org/10.1016/j.jcp.2007.10.027
.
Gunawan
P. H.
&
Lhébrard
X.
2015
Hydrostatic relaxation scheme for the 1D shallow water – Exner equations in bedload transport
.
Computers & Fluids
121
,
44
50
.
https://doi.org/10.1016/j.compfluid.2015.08.001
.
Harten
A.
,
Lax
P. D.
&
van Leer
B.
1983
On upstream differencing and Godunov-type schemes for hyperbolic conservation laws
.
SIAM Review
25
(
1
),
35
61
.
https://doi.org/10.1137/1025002
.
Jin
S.
&
Liu
J.-G.
1996
The effects of numerical viscosities
.
Journal of Computational Physics
126
(
2
),
373
389
.
https://doi.org/10.1006/jcph.1996.0144
.
Krvavica
N.
,
Tuhtan
M.
&
Jelenić
G.
2018
Analytical implementation of Roe solver for two-layer shallow water equations with accurate treatment for loss of hyperbolicity
.
Advances in Water Resources
122
,
187
205
.
https://doi.org/10.1016/j.advwatres.2018.10.017
.
Kurganov
A.
&
Petrova
G.
2009
Central-upwind schemes for two-layer shallow water equations
.
SIAM Journal on Scientific Computing
31
(
3
),
1742
1773
.
https://doi.org/10.1137/080719091
.
Kurganov
A.
,
Noelle
S.
&
Petrova
G.
2001
Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton–Jacobi equations
.
SIAM Journal on Scientific Computing
23
(
3
),
707
740
.
https://doi.org/10.1137/s1064827500373413
.
Lax
P. D.
1973
Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock
.
Courant Institute of Mathematical Science, Philadelphia
.
LeFloch
P.
1989
Shock Waves for Nonlinear Hyperbolic Systems in Nonconservative Form
.
Institute for Mathematics and its Applications, University of Minnesota
,
Minneapolis
.
LeFloch
P. G.
&
Thanh
M. D.
2007
The Riemann problem for the shallow water equations with discontinuous topography
. Communications in Mathematical Sciences 5 (4),
865
885
. https://doi.org/10.4310/CMS.2007.v5.n4.a7.
LeFloch
P. G.
&
Thanh
M. D.
2011
A Godunov-type method for the shallow water equations with discontinuous topography in the resonant regime
.
Journal of Computational Physics
230
(
20
),
7631
7660
.
http://dx.doi.org/10.1016/j.jcp.2011.06.017
.
Li
G.
,
Caleffi
V.
&
Gao
J.
2014
High-order well-balanced central WENO scheme for pre-balanced shallow water equations
.
Computers & Fluids
99
,
182
189
.
http://dx.doi.org/10.1016/j.compfluid.2014.04.022
.
Majd
M. S.
&
Sanders
B. F.
2014
The LHLLC scheme for two-layer and two-phase transcritical flows over a mobile bed with avalanching, wetting and drying
.
Advances in Water Resources
67
,
16
31
.
https://doi.org/10.1016/j.advwatres.2014.02.002
.
Martínez-Aranda
S.
,
Ramos-Pérez
A.
&
García-Navarro
P.
2020
A 1D shallow-flow model for two-layer flows based on FORCE scheme with wet–dry treatment
.
Journal of Hydroinformatics
22
(
5
),
1015
1037
.
https://doi.org/10.2166/hydro.2020.002
.
Meurice
R.
&
Soares-Frazão
S.
2020
A 2D HLL-based weakly coupled model for transient flows on mobile beds
.
Journal of Hydroinformatics
22
(
5
),
1351
1369
.
https://doi.org/10.2166/hydro.2020.033
.
Murillo
J.
&
García-Navarro
P.
2010
Weak solutions for partial differential equations with source terms: Application to the shallow water equations
.
Journal of Computational Physics
229
(
11
),
4327
4368
.
https://doi.org/10.1016/j.jcp.2010.02.016
.
Murillo
J.
&
García-Navarro
P.
2013
Energy balance numerical schemes for shallow water equations with discontinuous topography
.
Journal of Computational Physics
236
,
119
142
.
https://doi.org/10.1016/j.jcp.2012.11.003
.
Murillo
J.
,
Martinez-Aranda
S.
,
Navas-Montilla
A.
&
García-Navarro
P.
2020
Adaptation of flux-based solvers to 2D two-layer shallow flows with variable density including numerical treatment of the loss of hyperbolicity and drying/wetting fronts
.
Journal of Hydroinformatics
22
(
5
),
972
1014
.
https://doi.org/10.2166/hydro.2020.207
.
Navas-Montilla
A.
&
Murillo
J.
2017
Overcoming numerical shockwave anomalies using energy balanced numerical schemes. Application to the shallow water equations with discontinuous topography
.
Journal of Computational Physics
340
,
575
616
.
https://doi.org/10.1016/j.jcp.2017.03.057
.
Navas-Montilla
A.
&
Murillo
J.
2019
Improved Riemann solvers for an accurate resolution of 1D and 2D shock profiles with application to hydraulic jumps
.
Journal of Computational Physics
378
,
445
476
.
https://doi.org/10.1016/j.jcp.2018.11.023
.
Paciorri
R.
&
Bonfiglioli
A.
2009
A shock-fitting technique for 2D unstructured grids
.
Computers & Fluids
38
(
3
),
715
726
.
https://doi.org/10.1016/j.compfluid.2008.07.007
.
Rogers
B. D.
,
Borthwick
A. G. L.
&
Taylor
P. H.
2003
Mathematical balancing of flux gradient and source terms prior to using Roe's approximate Riemann solver
.
Journal of Computational Physics
192
(
2
),
422
451
.
https://doi.org/10.1016/j.jcp.2003.07.020
.
Schijf
J. B.
&
Schönfled
J. C.
1953
Theoretical considerations on the motion of salt and fresh water
. In:
Proceedings Minnesota International Hydraulic Convention
.
http://resolver.tudelft.nl/uuid:5d1c2eb0-d51c-4b3c-ad77-a77513941c6c
.
Soares-Frazão
S.
&
Zech
Y.
2011
HLLC scheme with novel wave-speed estimators appropriate for two-dimensional shallow-water flow on erodible bed
.
International Journal for Numerical Methods in Fluids
66
(
8
),
1019
1036
.
https://doi.org/10.1002/fld.2300
.
Stecca
G.
,
Siviglia
A.
&
Toro
E. F.
2010
Upwind-biased FORCE schemes with applications to free-surface shallow flows
.
Journal of Computational Physics
229
(
18
),
6362
6380
.
https://doi.org/10.1016/j.jcp.2010.05.001
.
Swartenbroekx, C., Soares-Frazão, S., Spinewine, B., Guinot, V. & Zech, Y. 2010. Hyperbolicity preserving HLL solver for two-layer shallow-water equations applied to dam-break flows, Proceedings of the River Flow 2010 Conference. Braunschweig, Germany, pp. 1379–1387
.
Toro
E. F.
2009
Riemann Solvers and Numerical Methods for Fluid Dynamics
.
Springer-Verlag
,
Berlin
.
Ungarish
M.
2013
Two-layer shallow-water dam-break solutions for gravity currents in non-rectangular cross-area channels
.
Journal of Fluid Mechanics
732
,
537
570
.
https://doi.org/10.1017/jfm.2013.417
.
Vázquez-Cendón
M. E.
1999
Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry
.
Journal of Computational Physics
148
(
2
),
497
526
.
https://doi.org/10.1006/jcph.1998.6127
.
Woodward
P.
&
Colella
P.
1984
The numerical simulation of two-dimensional fluid flow with strong shocks
.
Journal of Computational Physics
54
(
1
),
115
173
.
https://doi.org/10.1016/0021-9991(84)90142-6
.
Zaide
D. W.
&
Roe
P. L.
2012
Flux functions for reducing numerical shockwave anomalies
. In:
ICCFD7
,
Big Island, Hawaii
, pp.
9
13
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).