A 1D shallow-flow model for two-layer flows based on FORCE scheme with wet–dry treatment

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.

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

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

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

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.

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.

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

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.

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.

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.

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

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.

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.

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

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

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

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

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.

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.

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.

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.

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.

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.

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.