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

In gravity currents such as saline gravity currents in submarine environments, turbidity currents in reservoirs or hypersaline plume generated by desalination plants, density may change in processes of entrainment, transport, and deposition, generating layers of different density (Fernández-Torquemada et al. 2009; Zordan et al. 2018). As density may change, it is appropriate to define the bulk density of each layer as a mixture of different species, each one with different concentration and relative density. Most of the gravity currents happen on complex topographies and over varied bed compositions and the erosion, transport, and depositional processes depend on the particle characteristics and on the flow hydrodynamics. Even though the literature on sediment transport and gravity currents is rich, the hydrodynamics of density-driven flows and two-phase flows and the internal flow structure are not completely understood (Zordan et al. 2018), and support through experimental (Capart & Young 1998) and numerical approaches is required (Spinewine et al. 2011).

Successful simulation of shallow water flows has been the target of many researchers for more than one decade and the approaches in the literature have continuously improved the quality of the numerical results and the reliability of numerical predictions. Advances based on Godunov-type methods (Godunov 1959) involve different topics: from the original idea of the C-property (Bermúdez & Vázquez-Cendón 1994) or well-balanced schemes (Greenberg & Leroux 1996) enforcing equilibrium (Rogers et al. 2003; Chacon et al. 2007; Murillo et al. 2007a, 2007b, 2008a, 2008b, 2009; Murillo & Garcia-Navarro 2010a, 2010b) and their extension to energy or exactly balanced schemes (Chinnayya et al. 2004; Noelle et al. 2007; Fjordholm et al. 2011; Murillo & García-Navarro 2012a, 2012b, 2014; Castro Díaz et al. 2013; Xing 2014; Navas-Montilla & Murillo 2015, 2016; Murillo & Navas-Montilla 2016; Valiani & Caleffi 2017), to the definition of appropriate entropy fixes in critical points where one of the eigenvalues vanishes (Harten 1984; Castro et al. 2001, 2004, 2005, 2006; Abgrall & Karni 2009; Murillo & Navas-Montilla 2016), the extension to very high order of accuracy (Castro & Toro 2008; Vignoli et al. 2008; Castro et al. 2012; Montecinos et al. 2012; Navas-Montilla & Murillo 2015, 2016) or the use of dynamically adaptive grid generation (Liang et al. 2004; Rogers et al. 2004; Lee et al. 2011). There is still room for improvement in new areas, such as the problems related to numerical shockwave anomalies (Cameron 1966; Emery 1968; Zaide & Roe 2012; Navas-Montilla & Murillo 2017), that still claim our attention.

In general, for both one- and two-layer shallow water systems of equations, the preservation of steady state solutions for general moving water equilibrium (Noelle et al. 2007; Fjordholm et al. 2011; Castro Díaz et al. 2013; Xing 2014) and water height positivity are a challenge (Cozzolino et al. 2012; Xing & Shu 2014; Zhang et al. 2014; Haleem et al. 2015; Liang & Smith 2015; Barzgaran et al. 2019). Non-physical negative water height becomes traumatic as the eigenvalues do not determine the time step size. This is a result of the non-purely hyperbolic characteristic of the system of equations (Xing et al. 2010; Xing & Zhang 2013; Xing & Shu 2014). Imposing case-dependent threshold values for the water depth, the limitation of the computational domain in the computation of the flow advance over dry bed is the current tendency to avoid numerical difficulties (Delis et al. 2011a, 2011b). Other techniques involve velocity-based limiters (Delis et al. 2011a, 2011b; Vater et al. 2015) to control the stability at wet/dry fronts.

In this work we focus on the combination of two solvers, which are extensions of the well-known HLL solver developed by Harten et al. (1983) and the Roe solver (Roe 1981), maybe the most disseminated approximate Riemann solvers for the resolution of shallow flows. Both methods provide an approximate flux at cell interfaces using an intermediate or star region (Toro 2009), ⁠, between cells. In the HLL, the numerical flux arises directly from the definition of a control volume that is centered on a cell interface and integrated over space and time. For the Roe method, a wave decomposition of the Jacobian matrix has to be defined. In contrast to the HLL solver, where alternative sets of wave estimates can be provided (Davis 1988; Einfeldt 1988), the approximate Roe solver provides directly a set of wave celerities that will determine the structure of the numerical solution. Being more computationally expensive, the Roe scheme is usually more suitable than the HLL in terms of accuracy, as the latter only considers two waves, that result in an increased amount of numerical diffusion. The HLLC solver (Toro et al. 1994) incorporates a contact wave in the solution leading to a competitive solver, but is not easily extended to any set of equations. It is also worth saying that the Roe scheme is not either easily applicable to any set of conservation equations as appropriate averages have to be defined (Rosatti et al. 2008; Murillo & García-Navarro 2010a, 2010b; Murillo et al. 2012).

In practical applications of shallow flows, the evaluation of the source term is an aspect of great difficulty. It is possible, for instance, to define Riemann problems (RPs) with large gradients in the bed step that may lead to gross estimations when evaluating its variation in time. In this work, two augmented solvers, which will be referred to as ARoe (augmented Roe) (Murillo & Garcia-Navarro 2010a, 2010b; Murillo & Navas-Montilla 2016) and HLLS (Murillo & García-Navarro 2012a, 2012b) solvers, including the extra wave associated with the presence of bed discontinuity, are used to analyze the effect of the integral approximations made over the source term. They were designed to consider that, in the presence of source terms, the constant state in the approximate solution, ⁠, is no longer single-valued. In general problems considering a bottom step, where water elevation is discontinuous, the solutions are bi-valuated at the discontinuity. Both solvers have been successfully used in complex cases described by the one-layer shallow water equations (SWE), as they provide a comprehensive description of the inner states defined in the presence of wetting/drying fronts in steady/unsteady solutions and applied to complex rheology, where the well-balanced property must be redefined to provide accurate stop-and-go triggering mechanisms (Murillo & García-Navarro 2011; Juez et al. 2014). Also, accurate results are obtained for both solvers if using appropriate fixes in the vicinity of critical points (Murillo & Garcia-Navarro 2010a, 2010b; Murillo & Navas-Montilla 2016).

The two-layer SWE represents a complex problem observed in nature driven by internal and external waves, making the application of numerical methods that can assist in the formulation of all the waves of the system, as in the Roe method, desirable. In the two-layer SWE, approximate eigenvalues cannot be described explicitly and numerical treatment of critical cases is of importance. In real flows, intense mixing of the two layers produces shear layer Kelvin–Helmholtz instabilities that may lead to intense mixing of the two layers. Loss of the hyperbolic character of the equations (Abgrall & Karni 2009) is associated with this phenomenon. This problem can be overcome by performing a non-trivial uncoupled treatment of the layers (Bouchut & Morales de Luna 2008), by incorporating in the equations a quadratic friction law between layers that returns the system to the hyperbolic region (Castro et al. 2011) or by remaining close to the boundary of the hyperbolicity (Krvavica et al. 2018). It is worth mentioning the work in Spinewine et al. (2011), where Riemann solvers derived from the HLL formalism were presented including strict preservation of static configurations. The coupling between the two layers did not appear in the Jacobian matrix for the upper layer, allowing a set of explicit and real wave celerites to be obtained.

On the other hand, when the hyperbolic character of the equations is lost, HLLS is still applicable to the unmodified two-layer system. The apparent simplicity of the HLLS solver may hide its numerical capabilities, as it contains the extraordinary power of the approximate Roe's Jacobian matrix. It allows the HLL to include an extra wave accounting for the source term, regardless of the hyperbolic character of the system. Its simplicity makes the HLLS a good candidate to manage wet/dry fronts, allowing positivity over the water depths to be enforced in a simple way, without requiring any tuning parameter.

In the two-layer model equations, layers of different density are considered, resulting from the transport of material at different concentrations in each layer, but assuming a constant value in each one. In environmental problems, transport of mass including spatial variations is frequently analyzed, where the definition of suitable solver is a not minor task (Fernández-Nieto & Narbona-Reina 2008; Murillo et al. 2008a, 2008b; Bai & Jin 2009; Morales et al. 2009; Murillo et al. 2009, 2012; Juez et al. 2013, 2016; Franzini & Soares-Frazao 2018; Gordillo et al. 2020). It is important to remark that the two-layer shallow flow model cannot be expressed in conservation-law form, making unfeasible the use of flux-based solvers, that is, the set of eigenvalues and eigenvectors cannot be determined through a Jacobian matrix derived from a flux function, as done in the one-layer shallow flow model. The same difficulty appears when modeling flow in elastic vessels. In Murillo et al. (2019), the Roe and HLL schemes were adapted to deal with test cases involving extreme collapsing conditions in elastic vessels. To achieve this goal, the system is transformed to provide a conservation-law form, allowing the proper definition of Rankine–Hugoniot conditions. Here, the range of applicability of the presented combination between the ARoe and HLLS is applied to the two-layer SWE with local variable density. Density variation in each layer is the result of variations of a bulk density driven by the mixture of different constituents (Juez et al. 2016).

In order to present and justify the numerical techniques presented here, the structure of the manuscript starts with the mathematical model, where first the integral and next the differential formulations are presented for the -split two-dimensional two-layer SWE. The differential formulation is commonly presented to define the eigenstructure of the problem. Here, differential formulation will also be used as the starting point in the definition of integral approaches over the source terms. We will focus on defining an integral formulation by defining suitable discretizations in scenarios of interest, arriving finally at discretizations that will prove successful when demanding energy balanced solutions, with independence of the numerical solver selected. This will allow us to alternate between HLLS and ARoe solvers depending on the particular needs of the case. In the next section, the linear approximate RP together with a suitable flux function, that will allow description of the updating scheme in terms of a wave-propagation function or in terms of fluxes, will be defined. Once the approximate RP is defined, the ARoe and HLLS will be recalled, respectively, presenting the approximate Jacobian Roe matrix used in both solvers, presenting next the correction/definition of wave speed estimates for ARoe and HLLS solvers. Analysis of transcritical cases is of relevance as they may produce unstable results that have to be properly circumvented. These results are translated to the two-layer system. Next, positivity fixes required to ensure an accurate tracking of wetting/drying fronts in the context of the HLLS solver are presented. The choice between the HLLS or the ARoe schemes in each RP places emphasis on the evolution of the transported materials and is also discussed. Finally, the extension to two-dimensional problems is defined.

Numerical results involve test cases of great difficulty and are presented. A wide variety of test cases are presented, confirming the properties of the combination presented here. Such cases comprise exactly balanced scenarios in subcritical and subcritical-transcritical scenarios, dam-break problems over bed variations and wet/dry fronts, non-hyperbolic situations and transcritical exchange flows with loss of hyperbolicity. The effect of a perturbation in the bulk density and the resulting oscillatory effect is shown. Finally, a 2D dam break around an obstacle that illustrates the stability and robustness of the numerical schemes in the presence of non-uniform bed topography and wetting/drying fronts for both layers, is presented.

Instead of addressing directly the 2D two-layer shallow water model, in this section, the -split 2D two-layer shallow water equations will be explored first. In this way, the relevant information in the coordinate y is retained while reducing the complexity of the techniques presented here. The resulting analysis of this set of equations can be straightforwardly translated to 2D without loss of generality and no numerical approximations by exploiting the rotational invariance of the governing equations, avoiding redefinition of the numerical fluxes depending on the mesh characteristics (Toro 2009). The coordinate x will be the normal direction (normal to the intercell boundary in the computational domain) and coordinate y will be the tangential direction (parallel to the intercell boundary in the computational domain).

Figure 1 depicts two control volumes, one for each layer. The bed boundary is fixed in time and a unit width configuration is used. Index 1 refers to the upper layer and control volume 1. Index 2 refers to the lower layer and control volume 2. In each control volume the velocity is depth averaged. The relevant integral formulation of the interactions between models derives from the depth-averaged equations for the conservation of mass and momentum.

In lower control volume is the bathymetry function and ⁠, ⁠, are the water depths and averaged velocities in x and in y of layers = 1,2, respectively. Herein, function is the surface level of the lower layer 2 and is the surface level of the upper layer 1. Furthermore, the mathematical model adopted allows a variable density in space in each layer. As a result, a variable density ratio between layers can be defined.

Prior to defining or applying any approximate solver, the differential form of the two-layer equations in (38) must be integrated ensuring that equilibrium is recovered in any control volume defined in a computational domain. Approximate solvers provide solutions to general initial value problems and must also guarantee equilibrium in steady state conditions.

To that end, it is necessary to analyze the physical meaning of the relation of fluxes and source terms. At the bed step discontinuity, a contact wave appears and a choice regarding energy/momentum conservation across such wave must be made. In Valiani & Caleffi (2017), the authors claim that the conservation of energy across the bed step contact wave seems to be the only approach consistent with the classical SWE (LeFloch et al. 2011) and the momentum balance can also be satisfied.

Now that the energy functions have been defined, they can be integrated in the computational domain. Their integration will be based on keeping discrete equilibrium while ensuring that momentum balance is also satisfied. In order to clarify the integration steps, quiescent equilibrium is examined first. Next, exact equilibrium under non-zero velocity conditions is analyzed. The integral relations provided here will be used here to guarantee equilibrium in steady state condition with independence of the approximate solvers used.

is a piecewise constant approximation of the solution at time of the conserved variables in cell i, is the length of the computational cells considered constant, is the time step, and and are the numerical fluxes, yet to be defined, which are computed solving the RPs at the interfaces by means of a suitable Riemann solver.

With independence of the approximate solver used, the consistency condition must be satisfied, i.e., that the integral of the solution of the linearized RP over a suitable control volume must be equal to the integral of the exact solution over the same control volume. Figure 2 shows an RP, with initial values and a control volume given by the time interval and the space interval ⁠, where ⁠, being the minimum and maximum wave velocities, respectively, in the domain at ⁠. Numerical estimations for and will be defined when discussing the wave-speed estimates of the problem. The notation L and R for the waves at the discrete level is retained here, as in general and can be selected using a wide variety of combinations of the initial information stored in cells i and ⁠.

One result of Roe's linearization is that the resulting approximate Riemann solution consists of only discontinuities and is constructed as a sum of jumps or shocks. The solutions for are governed by the celerities in and each one consists of regions connected by waves, one of them steady, with celerity ⁠.

Unlike the ARoe solver, in the HLLS solver, which we will use in this work, only two inner states are defined, limiting the number of waves involved in the solution. The selection of the HLLC solver (Toro et al. 1994) or the HLLCS (Murillo & García-Navarro 2012a, 2012b) solver, can be considered a better choice, as both of them avoid excessive smearing of contact waves and vortices in 2D flows if compared with using the HLL or HLLS solvers. Despite these improvements, the HLLS solver has been selected, as it retains a simpler and more efficient analysis of the inner states in wetting/drying processes for both layers while introducing a sufficient level of numerical diffusion in cases of loss of hyperbolicity.

In order to generate the intercell fluxes it is necessary to compute the speeds and ⁠. Suitable approximations are discussed next.

Wave-speed estimates are analyzed in this section. When roots in are real, the ARoe scheme provides the set of approximate eigenvalues or waves of the system, internal and external ones, while the HLLS requires estimation of the external ones. In any case, the time step will be restricted by the largest approximate wave celerity selected when using one specific solver.

Roe's linearization provides an approximate solution constructed as a sm of jumps or shocks, leading to a good approximation for contact and shock waves. In rarefaction waves, a continuous change in flow variables appears, and approximations based on shocks may lead to inaccurate results, especially in transcritical cases. The accuracy of the numerical predictions for the one layer is certainly improved if using the Harten–Hyman entropy fix in (Harten 1984; Toro 2009) in transcritical dam-break problems as well as in steady solutions. This procedure also allows reproducing exactly the sonic point in the presence of source terms (Murillo & García-Navarro 2012a, 2012b).

Although the two-layer model presents a rich variety of waves, when using the HLLS, only two inner states are defined. The application of the HLLS does not require a decomposition in waves, only computation of appropriate wave bounds and ⁠. Direct wave speed estimate of the external waves, and ⁠, is the simplest method providing minimum and maximum signal velocities.

The relevant conclusion is that HLLS can be used to strongly simplify our numerical solution and in cases where the set of eigenvalues is not provided, that is, hyperbolicity is lost, generating suitable values for and ⁠. Loss of hyperbolicity is a challenge when numerically dealing with two-layer flows and is a limitation for the application of more complete solvers, such as the Roe solver. Whether appropriate wave estimates are provided, the HLLS can be directly applied as the Jacobian matrix is defined even when the roots of the characteristic polynomial in (93) become complex.

In order to overcome possible numerical anomalies, a critical flow-fix correction is proposed, based on the evaluation of the characteristic polynomial in cells i and in (39).

The numerical integration of the source terms can be performed exactly in steady state RPs, but when moving to transient variations, gross estimations of the source term may lead to unphysical results, as negative values of water depth, with independence of the time step selected. A positivity fix of the source terms can be enforced by means of the choice of a suitable control volume in the cell averaging step, or by analyzing each inner state of the approximate solution separately. This procedure has been applied successfully to the one-layer shallow water equations in Murillo & Navas-Montilla (2016) for the ARoe scheme, but cannot be easily extended to the two-layer shallow flows where a large number of waves is included and an explicit evaluation of the approximate eigenvalues cannot be provided. This issue is bypassed here by using the HLLS solver, not only in cases where hyperbolicity is lost, but also in those problems where the water depth of either the upper or lower layer vanishes, as will be seen in the applications.

Using the definition of the inner states in (120) and their relation in (119), it is possible to enforce positivity conditions over the two inner states that shape the approximate solution given by the HLLS solver. Positivity conditions enforce positive values of water depth for both layers at the left and right side of the cell interface, and ⁠, respectively. Both states depend on source terms and ⁠, and therefore source corrections cannot be directly done by modifications over each value separately. Nevertheless, it is still possible to analyze separately and ⁠, allowing a direct control over each layer. If ⁠, the following results can be derived for layers ⁠, and for with ⁠:

• Positive values of require the following limit over

  

  (139)
• In the case that becomes negative, can be replaced by ⁠, leading to the following left and right inner states

  

  (140)
• ensuring positive values of water depth on the right side of the plane solution. In the case that the left side of the initial solution is dry, if setting the discharge across becomes nil,

  

  (141)

• and the modification of the source term acts as the imposition of a reflecting wall.

• Positive values of lead to the following limit over

  

  (142)
• In the case that becomes negative, can be replaced by with the following consequences

  

  (143)

• ensuring positive values of water depth in both sides of the solution. In the particular case, with the left side of the initial solution dry, ⁠, if enforcing ⁠, intercell flux given by (121) becomes zero and the modification of the source strengths acts as the imposition of a reflecting wall in the wet/dry RP.

If integral approaches of the source terms are modified to ensure the condition in (144) positive values of the space solution are ensured and if one side of the initial solution is dry, a reflexion boundary condition is generated independently over each layer.

In the previous sections the HLLS and the ARoe scheme with application to the two-layer SWE have been detailed. The HLLS can be used even in cases where this system looses the hyperbolic character, by using appropriate wave estimates and ⁠. The exclusive use of the HLLS in two-layer flow solver provides inaccurate and excessively diffusive results, as in the hyperbolic region the problem is governed by two internal and two external waves. Only two bounding waves cannot accurately define the solution.

The use of the HLLS is, therefore, enforced here in cases of loss of hyperbolicity and in wetting-drying RPs when one of the water depths at each side of the RP vanishes, allowing a systematic treatment of wetting-drying fronts without imposing any tuning parameter. In any other case, the ARoe scheme is applied.

As we are considering the transport of species with variable concentration in each layer, the numerical update of their value over time must be carried out using a unified framework for the two numerical schemes used, HLLS and ARoe, for the resolution of the flow.

The 2D extension to the two-layer shallow water system is presented, and the updating numerical scheme written in terms of an unstructured mesh can be easily written for a Cartesian mesh if desired. The definition of the updating numerical scheme, based again in intercell flux approximations, arises from exploiting the rotational in variance of the governing equations (Toro 2009).

The suitability of the combination of the HLLS and ARoe in cases of loss of hyperbolicity or in wetting-drying fronts when one of the water depths of one or both layers vanishes is shown in the applications section. A wide variety of test cases is presented confirming the properties of this combination. They include exactly balanced scenarios in subcritical and subcritical-transcritical scenarios, dam-break problems over bed variations and wet/dry fronts for both layers, non-hyperbolic conditions, transcritical exchange flow with loss of hyperbolicity. Despite the complexity of the test cases presented here, accurate and stable simulations are guaranteed, ensuring positivity of the solution and not requiring a decrease of the time step. Most of the numerical tests presented in this work assume constant density in each layer and constant density ratio, as the test cases presented have been selected from the literature published for two-layer flows. But in environmental flows, density may change, making it necessary to consider for each layer a mixture of different species, each one with different concentration and relative density. In order to demonstrate the capability of the tools presented here in cases where initial conditions involve variable density in space, a new numerical test is presented in this work. This test case, that in conditions of constant density ratio remains in static equilibrium, is modified introducing a mixture of three constituents leading to horizontal variations of density and density ratio. The solution generates an oscillatory movement that is analyzed using cyclic boundary conditions.

Figure 7 shows the numerical solution using 100 cells. In all cases the unit discharge is constant. This result is independent of the integral approximation over the source term, as it is ensured by the properties of the solvers themselves. However, depending on the approximation done over the trust term, differences in the head energy are observed. Numerical integration in (82) (⁠⁠) deviates the total head energy with respect to the exact solution. The symmetry in the bed elevation is responsible for recovering the exact level downstream. When using the approach in (83) (⁠⁠) numerical integration conserves exactly at every cell the head energy, with independence of the grid refinement, as confirmed by results in Figure 8, for 1,000 cells.