Abstract

The two-layer problem is defined as the coexistence of two immiscible fluids, separated by an interface surface. Under the shallow-flow hypothesis, 1D models are based on a four equations system accounting for the mass and momentum conservation in each fluid layer. Mathematically, the system of conservation laws modelling 1D two-layer flows has the important drawback of loss of hyperbolicity, causing that numerical schemes based on the eigenvalues of the Jacobian become unstable. In this work, well-balanced FORCE scheme is proposed for 1D two-layer shallow flows. The FORCE scheme combines the first-order Lax–Friedrichs flux and the second-order Lax–Wendroff flux. The scheme is supplemented with a hydrostatic reconstruction procedure in order to ensure the well-balanced behaviour of the model for steady flows even under wet–dry conditions. Additionally, a method to obtain high-accuracy numerical solutions for two-layer steady flows including friction dissipation is proposed to design reference benchmark tests for model validation. The enhanced FORCE scheme is faced to lake-at-rest benchmarking tests and steady flow cases including friction, demonstrating its well-balanced character. Furthermore, the numerical results obtained for highly unsteady two-layer dambreaks are used to analyse the robustness and accuracy of the model under a wide range of flow conditions.

HIGHLIGHTS

  • In this work, a new 1D well-balanced FORCE scheme is proposed able to deal with steady and highly unsteady two-layer shallow flows, including wet–dry fronts.

  • The enhanced FORCE scheme is supplemented with a new hydrostatic reconstruction procedure, based on the definition of a virtual free surface and a virtual pressure surface for each layer, in order to ensure the well-balanced behaviour at the intercell edges for steady flows.

  • The proper treatment of the wet–dry fronts in both the upper and the lower layer is ensured by the correct balance between the convective fluxes and the integrated source terms at the wet–dry intercell edges.

  • Additionally, a method to obtain high-accuracy numerical solutions for two-layer steady flows including friction dissipation is also proposed in order to design reference benchmark tests for two-layer models validation.

  • The numerical results obtained with the well-balanced FORCE scheme for steady and highly unsteady two-layer benchmark tests are used to analyse the robustness and accuracy of the proposed model under a wide range of flow conditions.

INTRODUCTION

Some natural geophysical flows can be modelled as a layered system consisting of two shallow fluids of distinct densities. The density ratio between the lighter upper layer and the heavier lower layer can be close to unity, as in stratified marine currents with different salinity, or extremely large as in the case of mine tailing dams. Under the hypothesis of two immiscible fluids where the vertical characteristic dimension of the flow is much smaller than their horizontal extension, these kinds of geophysical flows have been often modelled using the depth-averaged shallow-flow equations in both layers and involving a free surface (Kurganov & Petrova 2009; Castro-Díaz et al. 2011; Spinewine et al. 2011).

When considering flow per unit width, the 1D two-layer shallow-flow problem, from the mathematical point of view, consists of a system of four conservation laws involving mass and momentum conservation for each layer. It has been studied by a large number of researchers (Ardron 1980; Armi 1986; Chen & Peng 2006; Bouchut & de Luna 2008), not only because of its applications but also because of the mathematical challenge itself. The two-layer system is conditionally hyperbolic and stable. Both layers interact by means of the hydrostatic pressure terms caused by the slope of both the free surface and the internal interface separating the two fluids, as well as by the momentum exchange at the interface due to friction.

One of the main features of the two-layer shallow-flow system is the loss of hyperbolicity which arises under a wide range of flow conditions (Ardron 1980; Armi 1986). The appearance of this complex behaviour is related to the density ratio between fluids and the velocity difference between layers (Spinewine et al. 2011). This conditionally non-hyperbolic character of the two-layer shallow-flow system has given rise to an extensive list of proposed solutions in order to avoid the appearance of numerical instabilities (Chen & Peng 2006; Bouchut & de Luna 2008; Kurganov & Petrova 2009; Abgrall & Karni 2010; Castro-Díaz et al. 2011; Mandli 2013; Chiapolino & Saurel 2018). Traditionally, numerical schemes based on the computation of the two-layer system eigenvalues, such as the HLLC or the A-Roe schemes, have failed in the treatment of the loss of hyperbolicity and avoiding the appearance of numerical instabilities. Furthermore, many of the reported numerical schemes for two-layer or multi-layer systems are complex to implement as the numerical scheme usually includes a specific method to deal with those instabilities. For instance, Chiapolino & Saurel (2018) solves the problem by introducing a term associated with the fluid compressibility which manages the instabilities changing the density of the layers; in addition, an interface artificial turbulence is added to increase diffusivity. The use of an entropy satisfying scheme (Bouchut & de Luna 2008) or an optimal amount of friction in a corrector step (Castro-Díaz et al. 2011) which keeps the flow inside the hyperbolic region has also been proposed.

The list of authors who have addressed the problem from the analytical point of view is much shorter (Schijf & Schönfled 1954; Farmer & Armi 1986). In the case of Farmer & Armi (1986), the solutions are restricted to the rigid-lid assumption (a uniform free surface) to a similar density of both layers, smooth bed functions and frictionless conditions. Up to now, the validation of new two-layer numerical schemes relies on a few analytic solutions for specific cases which, in general, do not include friction. In Krvavica et al. (2017), a two-layer shallow-flow model based on Roe-scheme including friction terms is validated using a salt-wedge case with a steady reference solution obtained from the ODE system previously proposed by Balloffet & Borah (1985).

This article presents a new 1D FORCE scheme (Toro 1996) for the two-layer shallow-flow equations, including friction source terms. A hydrostatic reconstruction procedure is also proposed to ensure the well-balanced behaviour of the scheme under steady flow conditions involving source terms. Furthermore, in order to design reference benchmark tests for two-layer model validation, a method to obtain high-accuracy numerical solutions for two-layer steady flows including friction dissipation is also proposed, based on the single-layer method reported by MacDonald (1996). The proposed scheme is faced with lake-at-rest, steady and transient flow benchmarking tests, proving its well-balanced character, robustness and accuracy under a wide range of fluid characteristic and flow conditions.

The text is organised as follows: in section ‘Mathematical model and governing equations’, the two-layer governing equations are stated, including the coupling pressure terms between layers and the friction source terms; section ‘Force numerical method for two-layer systems with source terms’ is devoted to the proposed well-balanced FORCE scheme supplemented with the hydrostatic reconstruction; in section ‘Reference solutions for two-layer steady flow’, the new procedure to obtain high-accuracy numerical solutions for two-layer steady problems is presented. The numerical results obtained with the proposed scheme for different steady and transient benchmarking tests are presented in section ‘Numerical results’. Finally, the conclusion of this work has been drawn in the last section.

MATHEMATICAL MODEL AND GOVERNING EQUATIONS

The one-dimensional two-layer shallow water equations result from depth integration of the 3D incompressible Navier–Stokes mass and momentum equations, according to the scheme shown in Figure 1. The lower layer, referred to as layer 2, is integrated between the bed surface and the interphase surface , whereas the upper layer, referred to as layer 1, is integrated between the interphase surface and the free surface . Considering a density ratio , being the density of the lighter upper layer and the density of the heavier lower layer , the following system of partial differential equations can be obtained for the two-layer flow per unit width:
formula
(1)
where is the acceleration of gravity, are the layer depths, are the depth-averaged velocities of each layer in the x direction and is the fixed bottom level. The boundary shear stresses and are those acting at the interface surface between fluid layers and at the fixed bottom. The pressure coupling terms and between both layers are expressed as:
formula
(2)
and they appear as a result of the vertical averaging of the three-dimensional problem (Toro 1997), transforming the system in a non-homogeneous problem even in cases over frictionless surfaces and flat bottom. Those terms characterise the two-layer flow equations as they describe the pressure momentum exchange between layers, i.e. accounts for the hydrostatic pressure of the upper layer 1 on the lower layer 2 and considers the interface surface slope, which is equivalent to a bed slope term for the upper layer 1.
Figure 1

Scheme for 1D two-layer shallow-flow problem.

Figure 1

Scheme for 1D two-layer shallow-flow problem.

There is a formal analogy with the 1D single-layer shallow water flow (Bouchut & de Luna 2008). The main difference is that the coupling terms, and , are not fixed but are also time-dependent. These source terms make the double-layer shallow water a complex flow to simulate numerically. It can easily lead to strong instabilities associated with physical phenomena as layer mixing and shear flows (Abgrall & Karni 2010).

The boundary shear stress between the lower layer and the fixed bed , assuming a turbulent flow for the lower layer, is defined with the squared flow velocity (Naef et al. 2006):
formula
(3)
being the friction turbulent factor, according to the Manning's law (Manning 1890) it is defined as , where is the bed roughness coefficient. The interface friction term is Chézy-type (Swartenbroekx et al. 2013), depending on the square of the relative velocity between the upper and the lower layer:
formula
(4)
being and the interface friction coefficients for the lower layer and the upper layers, respectively.
The 1D system (1) can be reordered in vector form as:
formula
(5)
where is the conserved variables vector and is the convective flux vector:
formula
(6)
The source vector includes the bottom slope, the friction stress terms at the fixed bed surface, the boundary shear stress at the interface boundary, as well as the pressure coupling terms and between both layers:
formula
(7)
This kind of formulation for the two-layer problem is referred to as the uncoupled formulation, since all the coupling terms between layers, i.e. and , are included in the source vector (Spinewine et al. 2011). The Jacobian matrix of the convective fluxes is defined as leading to the following quasi-lineal partial derivative system:
formula
(8)
with:
formula
(9)
being and the infinitesimal wave celerity in the upper and lower layer, respectively. The Jacobian matrix has four different real eigenvalues:
formula
(10)

The eigenvalues, and in Equation (10) are related to the upper layer whereas and are associated with the bottom layer. They are considered to control the numerical stability of explicit schemes through the CFL condition. If the convective fluxes vector is defined including the coupling terms and , leading to the coupled formulation of the two-layer problem, the eigenvalues of the Jacobian are not explicit (Abgrall & Karni 2010; Castro-Díaz et al. 2011) and can only be approximated when the density ratio . Closed-form solutions for eigenvalues of the coupled Jacobian matrix have been recently presented in Krvavica et al. (2018) and tested in Krvavica (2020). Furthermore, this coupled formulation which includes the pressure terms and into the convective fluxes, in order to concentrate in the Jacobian of the system all the relevant information, might lead to the loss of the hyperbolic character of the system (5) in a wide range of flow regimes (Spinewine et al. 2011).

FORCE NUMERICAL METHOD FOR TWO-LAYER SYSTEMS WITH SOURCE TERMS

The FORCE scheme (Toro 1996; Toro & Billett 2000; Toro et al. 2009; Dumbser et al. 2010) is an explicit-centred finite volume method (FVM) with proven numerical stability and robustness. For the sake of simplicity, a uniform spatial and variable temporal discretisation (Figure 2) is applied to solve the two-layer shallow water flow equations (Equation (1)), which are integrated spatially and temporally as follows:
formula
(11)
where the convective fluxes and are defined as:
formula
(12)
Figure 2

FVM sketch for two-layer FORCE scheme.

Figure 2

FVM sketch for two-layer FORCE scheme.

Assuming a piecewise constant representation for the conserved variables at the i cell for the time level n as:
formula
(13)
it is possible to define the numerical flux at the intercell edge separating cells i and as:
formula
(14)
Therefore, the 1D scheme for the conserved variables updating at the i cell can be expressed as:
formula
(15)
where is the time step and is the cell size. The source contribution, including bed slope, pressure and shear stress components, is approximated using a spatially integrated term defined as follows:
formula
(16)
The FORCE numerical flux at the intercell edges is formulated in a single step as a combination of the Lax–Friedrichs (first order and highly diffusive) (Lax 1954) and the Lax–Wendroff (second order and conditionally oscillatory) (Lax & Wendroff 1960) fluxes:
formula
(17)
where is a weighting parameter, is the Lax–Friedrichs flux at the intercell edge , computed as:
formula
(18)
and is the Lax–Wendroff numerical fluxes at the intercell edge , computed as:
formula
(19)
being and the set of conserved variables provided by the hydrostatic reconstruction (see section ‘Well-balanced hydrostatic reconstruction’) at the left and right side of the intercell edge , respectively.
The spatial integration of the source terms vector at the cell i, , is approximated as:
formula
(20)
where and account for the integrated pressure forces at the upper and lower layers of cell i, respectively; and and are the integrated friction forces at the bed and interface surfaces along cell i, respectively. The terms and are computed as:
formula
(21)
defining the virtual pressure surface for the upper and lower layers, and respectively, as:
formula
(22)

The values for the virtual pressure surface variation and , and the averaged layer depth and at cell i are provided by the hydrostatic reconstruction for the upper and lower layers. This reconstruction has been designed to ensure a well-balanced equilibrium of the fluxes at the intercell edges under static equilibrium conditions, as well as the correct treatment of the wet–dry fronts for both layers (see section ‘Well-balanced hydrostatic reconstruction’).

Furthermore, the explicit FORCE scheme requires a time step limitation to ensure the numerical stability of the method:
formula
(23)
being the eigenvalue the maximum of all the eigenvalues over the domain at every time and the Courant–Friedrichs–Lewy (Kubrusly & de Moura 2013) condition.

Well-balanced hydrostatic reconstruction

The original FORCE scheme (Toro 1996) was proposed for homogeneous systems, but including source terms into an enhanced formulation requires a careful treatment in order to ensure the well-balanced property of the numerical scheme when dealing with flow steady states. The well-balanced FORCE version proposed here, which balances the source terms integral with the intercell numerical fluxes, uses the hydrostatic reconstruction method (Audusse et al. 2004; Mahdavi & Talebbeydokhti 2009) as well as a correct treatment of the wet–dry front for both layers.

The main idea is to formulate the discrete source term at cell i and the numerical fluxes at the intercell edges and to achieve the balance between the convective fluxes and the source terms during the time step:
formula
(24)

The condition (24) requires to reconstruct the discrete values of the set of conserved variables at the left and right sides of the intercell edge, and , respectively. Moreover, it is necessary to define the degrees of freedom on the integrated source term (Equation (20)), i.e. the virtual pressure surface variations (,) and the averaged depth ( and ).

There exist two different conditions which require the maintenance of the static equilibrium state. First, under null velocity and horizontal free surface level conditions, , the solution must remain in equilibrium if the density of both layers is equal, i.e. , regardless of the position of the interface surface (Figure 3(a)). Second, under null velocity and horizontal free and interface surfaces, and , the static equilibrium must remain unaltered for any value of the density ratio r (Figure 3(b)). The proposed hydrostatic reconstruction method is based on the proper definition of the virtual pressure surfaces and , as well as the virtual free surfaces and , for both the upper and the lower layers respectively.

Figure 3

Hydrostatic equilibrium conditions: (a) static equilibrium with density ratio and (b) static equilibrium with horizontal free and interface surfaces.

Figure 3

Hydrostatic equilibrium conditions: (a) static equilibrium with density ratio and (b) static equilibrium with horizontal free and interface surfaces.

For the upper layer, the virtual pressure surface at cells i and is defined as:
formula
(25)
and the virtual free surface at cells i and can be expressed as:
formula
(26)
which agree with the physical interface surface and the physical free surface , respectively (Figure 4). The average value of the virtual pressure surface for the upper layer at the edge can be computed as:
formula
(27)
Figure 4

Hydrostatic reconstruction of the upper layer under equilibrium condition.

Figure 4

Hydrostatic reconstruction of the upper layer under equilibrium condition.

The first component of the reconstructed conserved variables and , i.e. the upper layer depth, is defined as:
formula
(28)
From this, the second component of the reconstructed conserved variables at the intercell edge , associated with the upper layer discharge, is defined as:
formula
(29)
being and the upper layer velocities at cells i and , respectively. Under hydrostatic equilibrium, and there exists a constant free surface for both configurations depicted in Figure 3. Therefore, in this case . Replacing in Equation (24), the homogeneous numerical fluxes balance at cell i for the upper layer can be expressed as:
formula
(30)
Therefore, to ensure the static equilibrium, the averaged depth and the virtual pressure surface variation present in the pressure force term for the upper layer at cell i (21) should be defined as:
formula
(31)
and hence .
Analogously, for the lower layer, the virtual pressure surface at cells i and is:
formula
(32)
and the virtual free surface at cells i and is:
formula
(33)
The reconstruction scheme focusing on the lower layer is depicted in Figure 5. Note that the effect of the upper layer is included in the virtual pressure surface. The average value of the virtual pressure surface for the lower layer at the edge can be computed as:
formula
(34)
Figure 5

Hydrostatic reconstruction of the lower layer under equilibrium condition.

Figure 5

Hydrostatic reconstruction of the lower layer under equilibrium condition.

The third component of the reconstructed conserved variables and , the lower layer depth, is defined as:
formula
(35)
From this, the fourth component of the reconstructed conserved variables at the intercell edge , associated with the lower layer discharge, is:
formula
(36)
being and the lower layer velocities at cells i and , respectively. Under hydrostatic equilibrium, and the virtual free surface for the lower layer is constant for both equilibrium configurations depicted in Figure 3, hence also for the lower layer .
Following the same procedure as for the upper layer, the averaged depth and the virtual pressure surface variation present in the pressure force term for the lower layer at cell i Equation (21) should be defined as:
formula
(37)
and hence .

Wet–dry fronts treatment

In order to ensure a correct treatment of the wet–dry fronts in both layers, the pressure force source term and at the upper and lower layer Equation (21), respectively, can be decomposed as:
formula
(38)
and hence Equation (24) can be rewritten as:
formula
(39)
Therefore, considering a wet–dry configuration for the upper layer at the intercell edge defined as:
formula
(40)
in order to guarantee that the well-balanced property of the numerical scheme is maintained for the upper layer even at wet–dry conditions, it is enough to ensure that:
formula
(41)
independently of the values of and .
Similarly, considering a wet–dry configuration for the lower layer at the intercell edge defined as:
formula
(42)
ensuring the well-balanced property for the lower layer requires that:
formula
(43)
independently of the values of and .
The terms and in Equations (41) and (43) are the integrated source term at the intercell edge for the upper and lower layers, respectively, which can be expressed as:
formula
(44)
being and the integrated pressure terms Equation (38); and and the integrated friction term at the left side of the intercell edge for the upper and lower layers, respectively, which can be expressed as:
formula
(45)

REFERENCE SOLUTIONS FOR TWO-LAYER STEADY FLOW

A simple method, based on the procedure reported by MacDonald (1996) for steady single-layer flows, is presented in this section to obtain high-accuracy numerical solutions for 1D steady two-layer flow cases. In order to validate the proposed model, reference solutions obtained with this method for different steady flows are going to be compared with the numerical results obtained for the same flows with the well-balanced FORCE scheme.

Following MacDonald (1996), in steady states , therefore system (1) can be simplified to:
formula
(46)
being and the flow discharge per unit width of the channel at the upper and lower layers, respectively, and . Defining the dimensionless Froude number as for the upper layer and for the lower layer, system (46) can be simplified to:
formula
(47)
Energy equations for both layers can also be obtained with a spatial integration of (47):
formula
(48)
with:
formula
(49)

Starting from an arbitrary function of the upper layer depth , in order to obtain reference solutions for the lower layer depth and the bed elevation , system (47) is solved using an ODE-solver. In this work, the solver implemented in Matlab was used to obtain all the exact two-layer steady solutions. It is worth noting that the reference solutions obtained for and using this method represent high-accuracy numerical approximations to those exact solutions.

The necessary input flow features for each layer are the densities and , the unit flow discharges and , the inlet boundary values for the lower layer depth and bed elevation , the friction coefficient at the interface surface and the Manning roughness coefficient at the bed surface. The method is summarised on Figure 6.

Figure 6

Scheme for the two-layer exact steady solution method.

Figure 6

Scheme for the two-layer exact steady solution method.

Three steady two-layer flow cases are presented with different density ratios, flow regimes at upper and lower layers and friction coefficients with the purpose of showing cases of interest whose high-accuracy numerical solution has been obtained with the method proposed. All of them consider the spatial domain . The upper layer depth function used in all the cases was similar to the one proposed by MacDonald (1996) for the single-layer analytical case:
formula
(50)

The required input values for each two-layer steady flow case are shown in Table 1. Case 1 consists of two layers with nearly equal densities (), with the lower layer under supercritical conditions and the upper layer under subcritical flow regime. Case 2 consists of a two-layer steady flow with both layers under subcritical flow regime and density ratio . Both cases 1 and 2 are frictionless at the interface and the bed surfaces. Finally, Case 3 includes a density ratio and frictional momentum exchange at the interface and bed surface. The reference solutions obtained for each hypothetical case have been shown and compared with the numerical results in section ‘Two-layer steady state tests’.

Table 1

Input values for the two-layer steady cases considered

Case
1,150 1,000 1.0 1.0 1.0 0.4 0.0 0.0 
3,000 1,000 0.8 1.2 1.0 0.4 0.0 0.0 
3,000 1,000 0.4 1.0 1.5 0.4 0.01 0.04 
Case
1,150 1,000 1.0 1.0 1.0 0.4 0.0 0.0 
3,000 1,000 0.8 1.2 1.0 0.4 0.0 0.0 
3,000 1,000 0.4 1.0 1.5 0.4 0.01 0.04 

NUMERICAL RESULTS

In this section, the proposed numerical scheme for the two-layer system is used to solve a set of different problems. Three different conditions are considered in order to validate the numerical behaviour of the proposed scheme: the lake-at-rest benchmarking test, steady flows and highly transient dam-break flows.

Two-layer lake-at-rest tests

This benchmarking test aims to demonstrate the well-balanced character of the FORCE scheme supplemented with the proposed hydrostatic reconstruction, even including the presence of wet–dry fronts in both layers. Any numerical method is said to be well balanced if under lake-at-rest conditions, i.e. horizontal free surface and null velocity, no instabilities appear and the static equilibrium is preserved. The hydrostatic reconstruction (3.1) included in FORCE allows the method to fulfil this condition for any value of the parameter (17) or density ratio r.

A spatial domain is considered with . Both the bed level and the interface level are arbitrarily defined as shown in Figure 7, including different wet–dry configurations in both layers. The weighting parameter is set to , the CFL and the final simulation time are 0.5 and 2.0 s, respectively. Frictionless conditions were also set at the bed and interface surfaces.

Figure 7

Lake-at-rest initial conditions.

Figure 7

Lake-at-rest initial conditions.

Two cases are considered with different density ratio r between upper and lower layer. First, is set, i.e. two layers of equal density. Figure 8 shows the layer levels and discharges for . Under this condition, the hydrostatic pressure of the upper layer on the lower layer is exactly balanced by the term accounting for the variation of the upper layer depth. Therefore, the numerical scheme is able to hold the static equilibrium state without suffering from any perturbation and preserving the initial interface surface unaltered with a horizontal upper layer free surface, as well as null velocities in both layers.

Figure 8

Lake-at-rest with density ratio : (left) layer levels at and (right) layer discharges at .

Figure 8

Lake-at-rest with density ratio : (left) layer levels at and (right) layer discharges at .

Second, the density ratio between both layers is set to , i.e. an upper light fluid over a dense fluid. Under this condition, the hydrostatic pressure of the upper layer on the lower layer does not balance exactly the upper layer depth variation term , and hence the initial static equilibrium is lost at the first simulation stages. However, as time progresses, a new static equilibrium state is reached and maintained with horizontal interface levels. Figure 9 (top) shows the interface and free surface temporal evolution. For the final simulation time , the interface and free surfaces are horizontal, and the discharge is null in both layers (Figure 9 (bottom)).

Figure 9

Lake-at-rest with density ratio : (left) layer levels temporal evolution and (right) layer discharges temporal evolution.

Figure 9

Lake-at-rest with density ratio : (left) layer levels temporal evolution and (right) layer discharges temporal evolution.

Two-layer steady state tests

In this section, the two-layer steady state cases with reference solution described in section ‘Reference solutions for two-layer steady flow’ are simulated using the proposed well-balanced FORCE scheme. Case 1 and 2 are frictionless with density ratio and , respectively. Case 3 considers a density ratio and involves friction momentum exchange terms at both the bed surface and the interface surface. The set-up parameters for each two-layer steady state case have been reported in Table 1.

A spatial domain is considered with . The chosen CFL is for all the simulations and the weighted parameter is set to . At the left boundary cell a constant discharge is set at both layers, whereas at the right boundary cell , constant values for the layer depths and energies are also imposed. Each simulation is stopped once steady flow regimes are reached at both layers.

Figure 10 shows the initial condition for Case 1, together with the bed profile and the reference solution for the two-layer steady flow calculated in section ‘Reference solutions for two-layer steady flow’. An initial constant elevation is imposed for the interface, as well as for the free surface. Constant discharges and are also imposed at the initial time along the whole spatial domain. The interface friction coefficient is set null, as well as the Manning's roughness coefficient at the bed surface.

Figure 10

Case 1: Initial condition for the interface and free surfaces, bed profile and reference solution for the two-layer steady flow.

Figure 10

Case 1: Initial condition for the interface and free surfaces, bed profile and reference solution for the two-layer steady flow.

For the two-layer steady flow Case 1, numerical instabilities do not appear using the proposed well-balanced FORCE scheme. Figure 11 (left) depicts the steady numerical interface and free surfaces elevation compared to the reference solution. Although a reasonable agreement is observed for the interface and free surface elevation, some differences appear for the flow velocities (see Figure 11 (right)), especially in the lower layer.

Figure 11

Case 1: comparison of numerical and reference solution for (left) interface and free surface elevations and (right) upper and lower layer velocities.

Figure 11

Case 1: comparison of numerical and reference solution for (left) interface and free surface elevations and (right) upper and lower layer velocities.

In order to identify the influence of the mesh discretisation, a more refined mesh with is also tested. Table 2 shows the root mean square error (RMSE) for the numerical depth (,) and velocities (,) of both layers compared with the reference solution. The increase in the mesh refinement leads to slightly lower errors in all the considered variables, with exception of the lower layer velocity . Furthermore, the influence of the weighting parameter is also assessed. The reduction of , maintaining , reduces the diffusivity of the scheme and also leads to slightly lower errors in all the considered variables, with exception of the lower layer velocity . However, the mesh refinement is not able to avoid the differences in the lower layer velocity between the numerical and the reference solutions, nor the reduction of the weighting parameter .

Table 2

Root mean square errors (RMSE) for Case 1: influence of the mesh discretisation and the weighting parameter

Case 1
     
     
     
     
     
     
Case 1
     
     
     
     
     
     

Figure 12 shows the initial condition, bed profile and reference solution for the two-layer steady flow Case 2. Constant elevations are again imposed for both the interface and the free surfaces. Constant discharges and are also imposed at the initial time along the whole spatial domain. The interface friction coefficient and the Manning's roughness coefficient at the bed surface are null.

Figure 12

Case 2: Initial condition for the interface and free surfaces, bed profile and reference solution for the two-layer steady flow.

Figure 12

Case 2: Initial condition for the interface and free surfaces, bed profile and reference solution for the two-layer steady flow.

Figure 13 (left) depicts the steady numerical interface and free surfaces elevation compared to the reference solution. As in the previous case, although a good agreement is observed for the interface and free surface elevation (see Table 3) between the numerical results obtained with the well-balanced FORCE scheme and the high-accuracy solution obtained with the ODE-solver, marked differences appears again in the lower layer flow velocity (see Figure 13 (right)).

Table 3

Root mean square errors (RMSE) for Case 2 and Case 3 with and

Test
Case 2     
Case 3     
Test
Case 2     
Case 3     
Figure 13

Case 2: comparison of numerical and reference solution for (left) interface and free surface elevations and (right) upper and lower layer velocities.

Figure 13

Case 2: comparison of numerical and reference solution for (left) interface and free surface elevations and (right) upper and lower layer velocities.

Figure 14 shows the initial condition, bed profile and reference solution for the two-layer steady flow Case 3. As in the previous cases, constant elevations are imposed for both the interface and the free surfaces. Constant discharges and at the initial time are also imposed. The interface friction coefficient is set to , whereas the Manning's roughness coefficient at the bed surface is .

Figure 14

Case 3: Initial condition for the interface and free surfaces, bed profile and reference solution for the two-layer steady flow.

Figure 14

Case 3: Initial condition for the interface and free surfaces, bed profile and reference solution for the two-layer steady flow.

For this two-layer steady flow Case 3, numerical instabilities do not appear either using the proposed well-balanced FORCE scheme, despite friction terms being included. Figure 15 (left) shows the steady numerical interface and free surfaces elevation compared to the high-accuracy solution obtained with the ODE-solver, whereas Figure 15 (right) depicts the flow velocities at both layers. A good agreement is observed for both the interface and free surface elevation, as well as for the flow velocities (see Table 3).

Figure 15

Case 3: comparison of numerical and reference solution for (left) interface and free surface elevations and (right) upper and lower layer velocities.

Figure 15

Case 3: comparison of numerical and reference solution for (left) interface and free surface elevations and (right) upper and lower layer velocities.

Finally, Table 3 shows the RMSE in Case 2 and Case 3 with and for the main variables of both layers. The highest errors occurs for Case 2, specially for the layer averaged velocity . Case 3 shows lower errors despite it including the friction dissipation terms.

Two-layer dam-break test

In this section, the proposed well-balanced FORCE scheme for two-layer flows is faced to a theoretical dam-break case. Dambreak is a highly transient benchmarking test widely used for the validation of shallow water numerical methods. In order to be able to use the single-layer analytical solution for comparison, the density ratio r is set to , i.e. both fluids with approximately the same density. A spatial domain is considered with . The final simulation time is , the CFL condition is set to and the parameter is equal to . Friction terms are null at the bed and interface surfaces. Both layers are set initially at rest, with a sharp discontinuity at for the upper layer depth:
formula
(51)
Figure 16 shows the numerical results for the interface and free surface levels at , and . The exact solution for the single-layer free surface level has also been plotted. As the densities of the upper and lower fluids are approximately equal, the free surface level evolution agrees with the exact single-layer solution, despite a marked change appearing at the interface surface. The dam-break shock wave going to the right direction moves with the same velocity as the single-layer exact solution. Furthermore, also the rarefaction wave going to the left direction agrees in height and velocity with the exact solution.
Figure 16

Dam-break temporal evolution with and : (top) , (centre) and (bottom) .

Figure 16

Dam-break temporal evolution with and : (top) , (centre) and (bottom) .

In order to analyse the effects of the density ratio variations, the idealised dam-break test case is also simulated considering , i.e. both fluid with exactly the same density, and accounting for a light fluid over a dense layer. Figure 17 shows the numerical results for the interface and free surface levels at , and . On one hand, for , numerical instabilities appear at the initial discontinuity region in the interface between both layers and they propagate as the dambreak progresses. These instabilities do not disappear when reducing the CFL as they are a consequence of the loss of the hyperbolic character of the mathematical system (Castro-Díaz et al. 2011). However, as the density of both fluids is equal, these oscillations do not affect the free surface, which agrees with the exact single-layer solution. On the other hand, when the density ratio between layers is , the level changes in the interface are less marked as the pressure term of the upper layer over the lower layer is reduced. In this case, the numerical solution is more stable. The free surface evolution cannot be compared with the exact solution for the single-layer problem in this case.

Figure 17

Dam-break temporal evolution for and , with : (top) , (centre) and (bottom) .

Figure 17

Dam-break temporal evolution for and , with : (top) , (centre) and (bottom) .

In order to analyse the effects of the weighted parameter on the stability and accuracy of the proposed numerical scheme, the same dam-break test case is simulated considering and varying , and . As shown in Figure 18, values of near to lead to diffusive solutions, as the first-order Lax–Friedrich flux dominates the numerical flux computation. When is reduced to , the Lax–Wendroff formulation, which is second order in time and space, has a more marked influence on composed numerical FORCE flux and hence the solution becomes less diffusive. Nevertheless, when is exactly , a predictable oscillatory behaviour was observed in the solution, derived from the pure second-order Lax–Wendroff formulation.

Figure 18

Dam-break temporal evolution with and , , : (top) , (centre) and (bottom) .

Figure 18

Dam-break temporal evolution with and , , : (top) , (centre) and (bottom) .

Finally, to demonstrate the capability of the well-balanced FORCE scheme to deal with two-layer flow configuration which causes loss of hyperbolicity of the coupled system, the idealised dambreak with is analysed. This case shows complex eigenvalues of the coupled Jacobian matrix (Krvavica et al. 2018), which make difficult the computation of the problem based on the coupled Jacobian. With the well-balanced FORCE scheme, numerical oscillations appear at the interface since the beginning and remain afterward (see Figure 17), but the scheme remains stable for the free-surface computation without showing a dramatic crash.

Furthermore, although it is not possible to totally avoid these oscillations, their amplitude can be controlled by increasing the diffusivity of the scheme. The most effective way to increase the numerical diffusivity is to decrease the CFL. The undesirable effects of lower CFL numbers on the solution accuracy can be partially corrected by tuning the weighting parameter . Figure 19 shows the dam-break temporal evolution with , CFL = 0.1 and on the left. Also the eigenvalues of the coupled Jacobian matrix have been plotted on the right to show that this case corresponds to a situation in which the loss of hyperbolicity exits. Comparison with Figure 17 shows that the numerical oscillation amplitude is reduced at the interface and the free surface is accurately predicted.

Figure 19

Dam-break temporal evolution with , CFL = 0.1 and : (left column) free surface and interface evolution, (right column) eigenvalues of the coupled Jacobian matrix.

Figure 19

Dam-break temporal evolution with , CFL = 0.1 and : (left column) free surface and interface evolution, (right column) eigenvalues of the coupled Jacobian matrix.

Experimental dam-break over movable bed

The proposed scheme is tested against a laboratory case consisting of a water dambreak over a flat movable bed. This experimental test (Spinewine & Zech 2007) consists of an idealised dam-break flow over a sediment flat bed made of sand with grain diameter , density and porosity . The experiment was carried out in a long and wide flume. The thickness of the sand layer was over the flume floor. Breaking of the dam was reproduced by the downward movement of a pneumatically actuated thin gate placed at the middle of the flume. The initial water level was upstream the gate and nil downstream. The temporal evolution of the free surface and bed surface were reported experimentally until after the gate opening each .

The simulations have been performed with a cell size , and . The turbulent friction term at the interface between the upper layer (water) and the lower layer (sand bed) is approximated using a variable friction coefficient , being a Manning-type roughness coefficient. The friction term between the sand layer and the flume floor is computed using a Chézy-type formulation, , being . The upper layer density was , whereas the bulk density of the sand layer is estimated as .

Figure 20 (left) shows the numerical results for the temporal evolution of the water free surface and bed level. The experimental results have been also plotted for comparison. In general, a good agreement with measured data can be found for all the times reported. The free water surface and the bed level are well predicted and the propagation velocity of the dam-break wave is accurately captured. Furthermore, the temporal evolution of the velocity in both layers has been also plotted in Figure 20 (right). The upper layer velocity maximum increases progressively as the dam-break wave progressed downstream the flume, whereas the lower layer averaged velocity remains stable along the dam-break flow.

Figure 20

Experimental dambreak over movable bed: (left) free surface and bed levels and (right) upper and lower layer velocities. From top to bottom, , , , , and .

Figure 20

Experimental dambreak over movable bed: (left) free surface and bed levels and (right) upper and lower layer velocities. From top to bottom, , , , , and .

CONCLUSIONS

This work demonstrates the possibility of simulating two-layer shallow flows with well-balanced FORCE schemes. In order to ensure this property, a new hydrostatic reconstruction procedure has also been proposed which ensures a correct balance between the homogeneous fluxes and the source terms at the intercell edges even under wet–dry conditions. Furthermore, in order to obtain new benchmarking tests for steady two-layer flows, a pioneering method to obtain high-accuracy reference solutions has been reported. This method allows to consider in the reference solution the dissipation terms caused by the frictional effects at both the interface separating fluid layers and the bed surface.

The numerical scheme has been tested against steady and transient two-layer cases. First, the FORCE scheme supplemented with the hydrostatic reconstruction has been faced to the lake-at-rest test considering different density ratios r between fluids and wet–dry configurations in both layers. The model is able to reach equilibrium states regardless of the value of r imposed, proving the well-balanced property of the scheme. Then, the FORCE model has been tested against three different two-layer steady flow cases with reference solution, which have been obtained by means of the proposed high-accuracy numerical procedure and including different flow regimes, density ratios and friction terms. The numerical results show good agreement with the reference solutions for the three steady cases and avoiding the appearance of numerical instabilities. The numerical scheme has also been applied to a two-layer dam-break test, considering different density ratios. The sensitivity of the FORCE scheme to the value of r has been analysed. Generally, the model becomes more stable as the density ratio decreases, i.e. as the lower layer increases its density respect to that of the upper layer. When the density ratio between layers is , the mathematical two-layer coupled system loses its hyperbolic character. However, this method is able to deal with the instabilities associated by increasing the numerical viscosity. The effects of the weighting parameter , balancing the first-order Lax–Friedrichs flux and the second-order Lax–Wendroff flux at the cell edges, have also been studied. Results demonstrate that values of smaller than the typical value are enough to avoid the oscillatory character of the pure Lax–Wendroff formulation and to reduce the excessive numerical diffusion associated with the Lax–Friedrichs scheme. Finally, the proposed scheme has been tested against an experimental dam-break over movable flat bed case, involving the advance of a wet–dry front. The scheme is able to predict with reasonable accuracy the free surface evolution and the movable bed changes. The obtained results show the big potential of FORCE scheme, supplemented with the proposed hydrostatic reconstruction method, to be used as modelled tool on real two-layer shallow-flow applications.

ACKNOWLEDGEMENTS

This work was developed under a CSIC-JAE grant in the framework of the project PGC2018-094341-B-I00, which is funded by MINECO/FEDER and by Diputacion General de Aragon, DGA, through Fondo Europeo de Desarrollo Regional, FEDER.

REFERENCES

Abgrall
R.
Karni
S.
2010
Two-layer shallow water system: a relaxation approach
.
SIAM Journal of Scientific Computing
3
,
1603
.
Ardron
K. H.
1980
One-dimensional two-fluid equations for horizontal stratified two-phase flow
.
International Journal of Multiphase Flow
6
(
4
),
295
304
.
Audusse
E.
Bouchut
F.
Bristeau
M.-O.
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
,
2050
2065
.
Balloffet
A.
Borah
D. K.
1985
Lower Mississippi salinity analysis
.
Journal of Hydraulic Engineering
111
(
2
),
300
315
.
Bouchut
F.
de Luna
T. M.
2008
An entropy satisfying scheme for two-layer shallow water equations with uncoupled treatment
.
Mathematical Modelling and Numerical Analysis
4
,
683
.
Castro-Díaz
M. J.
Parés-Madroñal
C.
Fernández-Nieto
E. D.
González-Vida
J. M.
2011
Numerical treatment of the loss of hyperbolicity of the two-layer shallow-water system
.
Journal of Scientific Computing
48
(
1–3
),
16
40
.
Chiapolino
A.
Saurel
R.
2018
Models and methods for two-layer shallow water flows
.
Journal of Computational Physics
371
,
1043
1066
.
Dumbser
M.
Hidalgo
A.
Castro
M.
Parés
C.
Toro
E. F.
2010
Force schemes on unstructured meshes ii: Non-conservative hyperbolic systems
.
Computer Methods in Applied Mechanics and Engineering
199
(
9
),
625
647
.
Krvavica
N.
Kožar
I.
Travaš
V.
Ožanić
N.
2017
Numerical modelling of two-layer shallow water flow in microtidal salt-wedge estuaries: finite volume solver and field validation
.
Journal of Hydrology and Hydromechanics
65
(
1
),
49
59
.
Kubrusly
C. S.
de Moura
C. A.
2013
The Courant–Friedrichs–Lewy (CFL) Condition, 80 Years After Its Discovery
.
Extras Springer/Birkhäuser
,
Berlin
.
Kurganov
A.
Petrova
G.
2009
Central-upwind schemes for two-layer shallow water equations
.
SIAM Journal of Scientific Computing
31
,
1742
1773
.
Lax
P. D.
1954
Weak solutions of nonlinear hyperbolic equations and their numerical computation
.
Communications on Pure and Applied Mathematics
7
(
1
),
159
193
.
Lax
P.
Wendroff
B.
1960
Systems of conservation laws
.
Communications on Pure and Applied Mathematics
13
(
2
),
217
237
.
MacDonald
I.
1996
Analysis and computation of steady open channel flow. PhD thesis, University of Reading, Reading, UK
.
Manning
R.
1890
On the Flow of Water in Open Channels and Pipes
.
Institute of Civil Engineers of Ireland
,
Dublin
.
Naef
D.
Rickenmann
D.
Rutschmann
P.
McArdell
B. W.
2006
Comparison of flow resistance relations for debris flows using a one-dimensional finite element simulation model
.
Natural Hazards and Earth System Sciences
6
,
155
165
.
Schijf
J. B.
Schönfled
J. C.
1954
Theoretical considerations on the motion of salt and fresh water
.
Proceedings Minnesota International Hydraulic Convention
.
Spinewine
B.
Zech
Y.
2007
Small-scale laboratory dam-break waves on movable beds
.
Journal of Hydraulic Research
45
(
suppl. 1
),
73
86
.
Spinewine
B.
Guinot
V.
Soares-Frazão
S.
Zech
Y.
2011
Solution properties and approximate Riemann solvers for two-layer shallow flow models
.
Computers & Fluids
44
(
1
),
202
220
.
Swartenbroekx
C.
Zech
Y.
Soares-Frazão
S.
2013
Two-dimensional two-layer shallow water model for dam break flows with significant bed load transport
.
International Journal for Numerical Methods in Fluids
73
,
477
508
.
Toro
E. F.
1996
On Glimm-Related Schemes for Conservation Laws.
Technical Report mmu-9602
.
Department of Mathematics and Physics, Manchester Metropolitan University
,
UK
.
Toro
E. F.
1997
Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction
.
Springer
,
Berlin, Heidelberg
.
Toro
E. F.
Billett
S. J.
2000
Centred TVD schemes for hyperbolic conservation laws
.
IMA Journal of Numerical Analysis
20
(
1
),
47
79
.
Toro
E. F.
Hidalgo
A.
Dumbser
M.
2009
Force schemes on unstructured meshes i: conservative hyperbolic systems
.
Journal of Computational Physics
228
(
9
),
3368
3389
.