# A Two-Layer Hydrostatic-Reconstruction Method for High-Resolution Solving of the Two-Layer Shallow-Water Equations over Uneven Bed Topography

Authors: Jiang Lin; Bing Mao; Xinhua Lu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In this study, attention is focused on numerically solving the two-dimensional, two-layer nonlinear shallow-water equations (2LSWEs) over uneven bed topography. A two-layer hydrostatic-reconstruction method (2LHRM) is proposed for face value reconstructions. A numerical model is then developed based on the 2LHRM in the framework of finite-volume methods using the slope-limited centered (SLIC) approximate Riemann solver for flux evaluations. The validations against various benchmark tests show that the 2LHRM by working with the SLIC scheme predicts robust solutions around large gradients and is able to simulate lake-at-rest solutions without using any ad-hoc techniques.


---

## Body

## 1. Introduction

The last two decades have seen extensive use of upwind Godunov-type and centered finite-volume methods to solve the one-layer nonlinear shallow-water equations (abbreviated as 1LSWEs hereafter), of which the former is dominantly used to date. Numerical fluxes associated with these two types of numerical methods rely on reconstruction of cell interface values (Riemann states). A face value reconstruction algorithm may affect a numerical scheme to correctly predict stationary flow at rest (e.g., C-property [1]), as well as to determine quantitative accuracy and model robustness in capturing spatially large gradients. It is found that a single-valued definition of bed elevation at the cell-interfaces is convenient and of vital importance in devising a C-property satisfying numerical scheme [2–5]. In solving the 1LSWEs, the hydrostatic-reconstruction method (HRM) proposed by Audusse et al. [3] is widely adopted; in this method, single-valued inter-cell bed elevations are dynamically calculated based on the reconstructed inter-cell water surface elevation and water depth values. By working with appropriate approximate Riemann solvers and source term discretization method, the HRM helps construct a number of robust and well-balanced numerical schemes for solving the 1LSWEs (e.g., Liang and Borthwick [6]; Hou et al. [7]; Lu et al. [8]).As in solving the 1LSWEs, one encounters similar problems in numerically solving the two-layer shallow-water equations (2LSWEs) to predict lake-at-rest solutions, as well as to robustly capture spatially large gradients with high-resolution. Lee et al. [9] employed the surface-gradient splitting method, originally proposed for solving the 1LSWEs [6], to obtain a well-balanced scheme for solving the 2LSWEs. This scheme do maintain steady stationary flows at rest but experiences high spurious oscillations in the vicinity of discontinuities, as shown by Lee et al. [9]. Spinewine et al. [10] proposed an ad-hoc technique to switch from solving the water-level based continuity equation under a quiescent flow condition (QFC), to solving the water-depth based continuity equation not under a QFC. This method is able to numerically maintain quiescent flow at rest, yet it violates the principle of mass conservation because the temporal variation of the interface is unconsidered. An alternative ad-hoc technique was proposed by Lu et al. [11], by limiting the diffusion part of numerical fluxes so that a zero mass flux is predicted under a QFC. To achieve high-resolution numerical solutions, Abgrall and Karni [12] added two auxiliary equations to the conventional 2LSWEs with a relaxation parameter; their numerical results showed that this method generally works well, but spurious oscillation may occur when an inappropriate relaxation parameter is chosen and the predictions may be grid-dependent. Izem et al. [13] developed a finite-element discontinuous Galerkin method for solving the 2LSWEs with the local Lax–Friedrichs scheme for flux evaluations; their numerical tests show that spurious oscillations may be generated near the discontinuities even for a flat bed situation. Chertock et al. [14] and Frings [15] used a three-layer approximation for the 2LSWEs to enhance numerical stability; as studied by Chertock et al. [14], their method do help stabilize the model but cannot completely eliminate the numerical instability. Castro et al. [16] presented an approach via introducing artificial friction to suppress the interfacial oscillations. This method is simple to implement, yet further work needs to be done to evaluate if the introduced friction force contaminates the physical dissipations. Dudzinski and Lukác˘ová-Medvid’ová [17] developed a finite-volume evolution Galerkin scheme for solving the 2LSWEs, based on the theory of bicharacteristics [18]; a Newton–Raphson iteration is performed at each time step to compute the eigenvalues.As the HRM gains popularity and success in solving the 1LSWEs, while we found none of existing 2LSWEs solvers is developed based on HRM, we are curious that can HRM be extended for solving the 2LSWEs, and if so, how does it perform in terms of capturing discontinues and maintaining steady stationary flow at rest? Besides, lots of existing 2LSWEs solvers are developed based on upwind Godunov-type method which relies on evaluation of eigenvalues of the partial-differential equation system. Unlike the 1LSWEs, the eigenvalues of the 2LSWEs cannot be theoretically computed and thus in some sense nontrivial to implement. In viewing that the centered approximated Riemann solvers, e.g., the first-order centered (FORCE) and its second-order version, the slope-limited centered (SLIC) scheme [19] avoid using the eigenstructure information in evaluating the numerical fluxes and are capable of capturing discontinuities without oscillation, it is desirable to use them for developing 2LSWEs solvers. In this work, we extend the HRM originally proposed for solving the 1LSWEs to solve the 2LSWEs. For brevity, the extended HRM is called two-layer hydrostatic-reconstruction method (2LHRM) hereafter. A numerical model is developed based on the 2LHRM, and the SLIC scheme is invoked for flux evaluations in the framework of collocated-grid based finite-volume method. The performance of the 2LHRM and the developed model is then validated in terms of quantitative accuracy and C-property satisfying and shock-capturing properties.The remainder of this paper is organized as follows. Section2 presents the governing equations and 2LHRM. In Section 3, numerical tests are utilized to assess the performances of the developed 2LHRM and numerical model. Finally, discussions and conclusions are given in Section 4.


## 2. Governing Equations and 2LHRM

The 2LSWEs are useful for describing two-layer shallow-water flow motions (see Figure1). To purify the problem, as was done in previous relevant studies (e.g., Lee et al. [9]; Swartenbroekx et al. [20]; Hu et al. [21]), the flows in the layers are assumed to be immiscible and have constant densities in each layer. Following Lu et al. [11], we consider solving the following 2LSWEs system: (1)∂U∂t+∂F∂x+∂G∂y=S,with vectors defined by (2)
U
=
h
2
q
x
2
q
y
2
h
w
q
x
w
q
y
w
,
F
=
q
x
2
u
2
q
x
2
u
2
q
y
2
q
x
w
u
2
q
x
2
+
1
χ
u
1
q
x
1
u
2
q
y
2
+
1
χ
u
1
q
y
1
,
G
=
q
y
2
v
2
q
x
2
v
2
q
y
2
q
y
w
v
2
q
x
2
+
1
χ
v
1
q
x
1
v
2
q
y
2
+
1
χ
v
1
q
y
1
,
S
=
0
-
g
h
2
∂
z
2
∂
x
-
g
h
2
∂
z
2
∂
y
0
-
g
H
∂
z
2
∂
x
-
1
-
χ
χ
g
h
1
∂
z
1
∂
x
-
g
H
∂
z
2
∂
y
-
1
-
χ
χ
g
h
1
∂
z
1
∂
y
,where t denotes time; xl(l=1,2) are the Cartesian coordinates with x1=x and x2=y; hk(k=1,2) are the layer depths with k=1 and k=2 denoting the values of the lower- and upper-layers, respectively; qxk and qyk(k=1,2) are the discharges per unit width in the x- and y-directions, respectively; uk=qxk/hk,vk=qyk/hk (k=1,2) denote the layer-averaged velocity components; zb is the bed elevation above a reference level zr (see Figure 1); zk(k=1,2) denote the layer surface elevations with respect to zr and have the following relations: (3)z1=zb+h1,z2=z1+h2;χ=ρ2/ρ1 is the density ratio between upper and lower layers; and (4)H=h1+h2,hw=h2+1χh1,qxw=qx2+1χqx1,qyw=qy2+1χqy1.Note that the external shear stress (e.g., wind stress) at the upper-layer surface, the interfacial shear stress between the layers, and the bed friction force are omitted in this work to highlight the primary issues encountered in solving the 2LSWEs.Figure 1
Definition sketch of a two-layer shallow-water system (not to scale).To solve the 2LSWEs defined in (1) and (2), the finite-volume method with approximate Riemann solvers for flux evaluations is considered, which rely on reconstruction of face values. We here firstly summarize some fundamental requirements for a face value reconstruction algorithm to be satisfied, based on which, an approach is proposed for the reconstruction of face values in solving the 2LSWEs. Note that we do not consider wetting and drying processes in this work. Also, we only consider the QFC with χ≠1 in this section and thus both the upper-layer and lower-layer surfaces should be horizontally flat under a QFC. As for the QFC with χ=1, a discussion will be made in Appendix B.When solving the 1LSWEs (obtained by settingh2(x,y,t)≡0 in the 2LSWEs), it is required according to (3) that under general conditions (not only QFCs), e.g., in the x-direction, (5)z1i-1/2,jL,R=zbi-1/2,jL,R+h1i-1/2,jL,R;and under a QFC with z1(x,y,t)≡ω, in order to achieve a lake-at-rest solution, a reconstruction scheme should satisfy, e.g., in the x-direction, (6)z1i-1/2,jL=z1i-1/2,jR=ω,(7)h1i-1/2,jL=h1i-1/2,jR.Eqs. (5) and (7) lead to (8)zbi-1/2,jL=zbi-1/2,jR,meaning that the inter-cell bed elevation is single-valued. One can easily prove that (5)–(7), together, ensure lake-at-rest solutions in solving the 1LSWEs over an uneven bed topography, when the SLIC scheme is used for flux evaluations.When solving the 2LSWEs, except the requirements expressed in Eqs. (5)–(7), which are constraints on the lower-layer variables, the reconstructions in relation to the upper layer should satisfy similar requirements as, e.g., in the x-direction, under general conditions, (9)z2i-1/2,jL,R=zfi-1/2,jL,R+h2i-1/2,jL,R,where zf denotes the reconstructed inter-cell “bed elevation” for the upper layer, and under a QFC with a constant upper-layer surface elevation ζ, (10)z2i-1/2,jL=z2i-1/2,jR=ζ,(11)h2i-1/2,jL=h2i-1/2,jR.In addition to the constraints given in (5)–(11), physically, it is also required that under general conditions, (12)zfi-1/2,jL,R=z1i-1/2,jL,R,because the “bed elevation” of the upper layer is also the surface elevation of the lower layer.To realize the relations expressed in (5)–(7) in solving the 1LSWEs, Audusse et al. [3] proposed the HRM. In this method, z1 and h1, which are the water surface elevation and water depth in the 1LSWEs, respectively, are firstly reconstructed (e.g., by a first-order piecewise constant or piecewise linear method). Then, based on the reconstructed values, a single-valued inter-cell bed elevation zb is defined to satisfy (8), and after that, h1 and z1 are modified so that (5) is fulfilled. We here extend the idea of HRM from solving the 1LSWEs to 2LSWEs so that the relations given in (5)–(12) are satisfied. The extended method, i.e., the two-layer hydrostatic-reconstruction method (2LHRM), reconstructs the face values through the following successive procedures.Step 1.
Reconstruction of the upper-layer variables: in this step, the face values ofz2, h2, qx2, and qy2 are interpolated from those at the cell centers, e.g., in the x-direction, by using the monotonic upstream-centered scheme for conservation laws (MUSCL) [22] with the minmod slope limiter [23] for numerical stability(13a)fki+1/2,jL=fkni,j+12Ψri+1/2,jLfkni,j-fkni-1,j,(13b)fki+1/2,jR=fkni+1,j-12Ψri+1/2,jRfkni+2,j-fkni+1,j,where f can be any of the variables to be reconstructed; Ψ(r)=max0,min(1,r), with ri+1/2,jL=fi+1,jn-fi,jn/fi,jn-fi-1,jn, and ri+1/2,jR=fi,jn-fi+1,jn/fi+1,jn-fi+2,jn. Note that when setting Ψ(r)≡0, the first-order piecewise constant reconstruction scheme is recovered.Step 2.
Reconstruction of the lower-layer variables: in this step, the face values ofh1, qx1, and qy1 are reconstructed from (13a) and (13b) (with k=1). Then, the face value of z1 is computed from the relation z1=z2-h2 (see (3)). After that, we set zfL,R=z1L,R to satisfy (12).Step 3.
Calculation of cell-interface velocities by(14)ukL,R=qxkL,RhkL,R,vkL,R=qykL,RhkL,Rk=1,2.Step 4.
Calculation of the single-valued inter-cellzb as, e.g., in the x-direction, (15)zbi-1/2,j=Λz1i-1/2,jL-h1i-1/2,jL,z1i-1/2,jR-h1i-1/2,jR.Here, Λ(a,b)=a+b/2 is used. Once zb is computed, h1L,R is redefined as (16)h1i-1/2,jL,R=max⁡0,z1i-1/2,jL,R-zbi-1/2,j.Then, z1L,R and zfL,R are successively updated from (5) and (12), respectively.Step 5.
Becausezf may be modified in Step 4, we need to redefine h2L,R as (17)h2i-1/2,jL,R=max⁡0,z2i-1/2,jL,R-zfi-1/2,j.Then, z2 is recalculated from (9).Step 6.
Computation of layer discharges at the cell interface by(18)qxkL,R=ukL,RhkL,R,qykL,R=vkL,RhkL,Rk=1,2, where ukL,R and vkL,R are computed in Step 3. qxw and qyw are calculated according to (4) as, e.g.,(19a)qxwL,R=qx2L,R+1χqx1L,R,qywL,R=qy2L,R+1χqy1L,R.It is obvious that the reconstructed face values by using the 2LHRM satisfy all the relations in (5)–(12). As a result, lake-at-rest solutions are ensured because under a QFC, the numerical fluxes are zero since in (A.3)–(A.5) F(U~)≡F(U)≡0 and U~L≡U~R, and the source terms in (A.6a), (A.6b), and (A.6c) are also zero because under a QFC the upper- and lower-layer surfaces are horizontally flat and thus Θk=0(k=1,2) (see (A.7)). It is also meaningful to point out that the 2LHRM reduces to the HRM when the upper-layer does not exist (i.e., h2≡0 and z2≡z1).To gain a better understanding of the two-layer reconstruction approach, we illustrate the major reconstruction procedures in Figure2. Without loss of generality, we assume that after Step 3, (zb)L=z1L-h1L⩾(zb)R=z1R-h1R. According to the magnitudes of z1R and z2R obtained after Step 3 and (zb)i-1/2 (abbreviated as zbs hereafter) obtained in Step 4, three typical cases exist, i.e., z1R<z2R⩽zbs, z1R⩽zbs⩽z2R, and zbs⩽z1R<z2R (see Figure 2). We stress here again that we focus on simulations with all wet cells in this work and cases with vanishing layer depth(s) at cell centers are unconsidered; when dry cells are involved, more cases will exist. From Figure 2, it is clearly seen that after reconstruction, in comparison to the values obtained after Step 3, z2L and h2L do not change in all three cases; z2R does not change except when z2R<zbs; h2R does not change except when z1R<zbs; h1L and h1R change in all cases due to the single-valued definition of zbs. Note that for the tests studied later in Section 3, only the first case will be encountered.Figure 2
Sketch of the procedures of the two-layer reconstruction approach: (a), (b), (c) show the values after reconstruction Steps3, 4, and 5, respectively, for the case zbs⩽z1R<z2R. (d), (e), (f) show the results for the case z1R⩽zbs⩽z2R, and (g), (h), (i) for the case z1R<z2R⩽zbs. The pink filled circle show the elevation of zbf.
(a) (b) (c) (d) (e) (f) (g) (h) (i)Because an explicit scheme is used, the time step is constrained by the Courant–Friedrichs–Lewy (CFL) criterion:(20)Δt=Cr×minΔxmax⁡ξi-1/2,jL,ξi-1/2,jR,Δymax⁡ξi,j-1/2B,ξi,j-1/2T,where ξi-1/2,jL,r and ξi,j-1/2B,T are approximately local bounds of wave speeds [10, 11], and Cr is a predefined parameter.


## 3. Test of the 2LHRM

In this section, we test the performance of the 2LHRM. A numerical model is developed based on the 2LHRM using the finite-volume method. In the model, the slope-limited centered (SLIC) approximate Riemann solver is employed for numerical flux evaluations. The boundary conditions are treated using the ghost cell method [11, 24]. The details of the numerical methods for the numerical model are presented in Appendix A. Several numerical tests are used in the following to verify the various properties of the proposed method. The tests include 1D transitional and transcritical steady flows over a hump, 1D internal dam-break waves propagating over a flat bottom, and 2D internal dam-break waves propagating over a flat and an uneven bed topography. The first test is chosen to check whether the proposed method predicts lake-at-rest solutions and the capability of the model in simulating two-layer shallow-water flow motions under different flow regimes, with flows moving in the same or different directions in the layers. The second and third tests are employed to assess the robustness of the model in simulating flows with discontinuities. The third test is also used to assess the order-of-accuracy of the developed numerical model. The coefficient Cr in (20) is set to be 0.5 in this work. Note that the approximate eigenvalues of the 2LSWEs system are not used in flux evaluations and are only adopted for time step controlling. For 1D simulations in the following, the grid number in the y-direction is set to be unity.


### 3.1. 1D Steady Flows over a Hump

1D two-layer shallow-water flows over an uneven bed topography are simulated here. Depending on the boundary conditions, the layer surface levels of the steady flows can be smooth, or nonsmooth with a hydraulic jump [12]. Three test cases are considered, including steady transitional flows with a hydraulic jump, steady transcritical exchange flows, and transcritical parallel flows with smooth layer surface levels. For all the test cases, approximate analytic solutions are available, which were derived by Abgrall and Karni [12] under a rigid-lid assumption. The computational channel is 6 m long (x∈[-3m,3m]) with the bed topography defined as (21)zbx=18cos⁡πx2+12,x⩽2m,0,else.Different boundary and initial conditions are applied. For the steady transitional flows with a hydraulic jump, (22)qx1x=-3m,t=0.0928m2/s,h1x=3m,t=0.9205m,qx2x=-3m,t=0.0928m2/s,h2x=3m,t=0.5795m,z1x,t=0=0.9205m,z2x,t=0=1.5m,and for the steady transcritical exchange flows, (23)qx1x=-3m,t=0.0928m2/s,h1x=3m,t=0.1617m,qx2x=3m,t=-0.0928m2/s,h2x=-3m,t=0.4311m,z1x,t=0=1.0m,z2x,t=0=1.5m,and for the steady transcritical parallel flows, (24)qx1x=-3m,t=0.0928m2/s,h1x=3m,t=0.1617m,qx2x=-3m,t=0.0928m2/s,h2x=3m,t=1.3338m,z1x,t=0=1.0m,z2x,t=0=1.5m.In all the test cases, χ=0.98 and g=10m2/s. For each simulation, when the globally relative change indictor R, defined below, drops smaller than 1×10-7, a steady solution is thought to be achieved. (25)R=∑k=12∑i=1nx∑j=1nyhki,jn+1-hki,jnhki,jn2.Here, nx and ny denote, respectively, the mesh cell numbers in the x- and y-directions, which are chosen to ensure numerical accuracy according to a grid-independence study. For instance, for the steady transitional flows with a hydraulic jump test case, a series of predictions at steady states on successively finer grids (i.e., nx=200,400,…,3200) are obtained. As seen from Figures 3(a) and 3(b) (note Figure 3(b) is a zoom-in view of Figure 3(a) around the hydraulic jump location), the simulated layer surface profiles with nx=1600 and nx=3200 are rather small and nx=1600 is chosen for this test case to achieve a balance between numerical accuracy and computational efficiency.Figure 3
Computed layer surface levels at steady states for the transitional flow with a hydraulic jump over a hump test: (a) whole simulation domain; (b) a zoom-in view in the regionx∈[0,1m]. (c) and (d) show the computed versus analytic profiles of interface elevation and layer discharges per unit width at steady states, respectively. The bed topography is plotted on the right y-axis.
(a) (b) (c) (d)In Figure3(c), the computed z1 at steady states is shown for the first test case. Also shown is the analytic solution of Abgrall and Karni [12], which is piecewise smooth with a hydraulic jump around x=0.48m. It is seen that the simulated results match the analytic solution fairly well, with slight deviations around the predicted hydraulic jump location. These deviations are thought to be ascribed to the limitation of the rigid-lid assumption adopted in deriving the analytic solution [12], and also the numerical error in predicting constant discharges in the layers [11]. In Figure 3(d), the computed layer discharges per unit width along with the analytic solution are displayed. A zoom-in view around the hydraulic jump location is given in Figure 3(d), from which it is seen that the predicted layer discharges agree with the analytic solution well, but deviations between the computed and analytic values exist near the jump location, with a maximum relative difference about 16%. This type of deviation in the predicted discharge has been commonly seen when an approximate Riemann solver is used in simulating steady one-layer or two-layer shallow-water flows over a hump involving a hydraulic jump (e.g., Zhou et al. [2]; Valiani and Begnudelli [25]; Lu et al. [11]). To test the stability of the numerical model and to check whether (20) is appropriate to dynamically determine the time step, a stability test is implemented here. For the first test case, the predicted interface profiles with nx=1600 using different values of Cr are shown in Figure 4. Note that all results in Figure 4 are plotted at steady flow states, except for the fact that in Figure 4(d) the results at t=1.0s are shown. It is found that when Cr<0.89 (see Figures 4(a)–4(c)), the numerical model is stable. When Cr⩾0.89, the numerical solution exhibits numerical instability as can be seen clearly in Figure 4(d). The stability study here indicates that the time step defined in (20) is suitable for dynamically determining the time step.Figure 4
Simulated interface profiles for the transitional flow with a hydraulic jump test withnx=1600 using Cr=0.3 (a), Cr=0.5 (b), Cr=0.8 (c), and Cr=0.89 (d). The results are shown at steady states, except in (d) where the result is illustrated at t=1.0s. A zoom-in view in the region x∈[-2.5,-1.5m] is shown in (d). The thin solid line at the bottom of each figure represents the bed profile and is plotted on the right y-axis.
(a) (b) (c) (d)Figures5(a) and 5(b) illustrate the computed layer surface levels and the discharges per unit width at steady states, along with the analytic solution of Abgrall and Karni [12] for the exchange flow test case. The corresponding results for the parallel flow test case are given in Figures 5(c) and 5(d). In the simulations, nx=200 is chosen according to grid-independence studies. From Figure 5, it is seen that the numerical results match the analytic solution reasonably well. According to Abgrall and Karni [12], the layer surface levels in these two test cases should be the same, which is well predicted by the developed model (c.f. Figures 5(a) and 5(c)).Figure 5
Computed (symbols) versus analytic (thick lines) layer surface levels (a) and discharges per unit width (b) at steady states, for the transcritical exchange flow over a hump test. The results for the parallel flow test case are shown in (c) and (d). The bed topography is shown at the bottom of each figure by the thin solid line and is plotted on the righty-axis.
(a) (b) (c) (d)At last in this test, we verify the performance of the developed model in predicting lake-at-rest solutions. The simulation setup is the same as the first test case except for the fact that at the western boundary the discharges in the layers are fixed to be zero. Figures6(a) and 6(b) present the computed versus analytic profiles of layer surface levels and discharges per unit width at t=10s, respectively. From Figure 6, it is seen that lake-at-rest solution is achieved.Figure 6
Simulated versus analytic layer surface levels (a) and discharges per unit width (b) att=10 s for the lake-at-rest test with χ≠1. The bed topography is plotted on the right y-axis.
(a) (b)


### 3.2. 1D Internal Dam-Break Wave Propagations

The 1D internal dam-break flow tests are widely used to investigate the robustness of numerical schemes in simulating flows involving large gradients (e.g., Bouchut and de Luna [26]; Dudzinski and Lukác˘ová-Medvid’ová [17]; Lee et al. [27]). The considered channel is 10 m long (x∈[0,10m]) and is horizontally flat (zb(x)≡0). The initial conditions are prescribed as (26)h1x,t=0=0.8sgn⁡x-5+1.0m,h2x,t=0=-h1x,t=0+2.0m,qx1x,t=0=qx2x,t=0=0.In the simulations, χ=0.7 and g=9.81m2/s. The simulation domain is discretized with nx=1600, chosen via a grid-independence study.Figures7(a) and 7(b) illustrate the time evolutions of the simulated interface and upper-layer surface, respectively. As seen, the present model is robust to capture sharp gradients. In Figure 7(c), the computed interface profile at t=1s is compared with those obtained by Bouchut and de Luna [26] via decoupled solving two 1LSWEs, and by Dudzinski and Lukác˘ová-Medvid’ová [17] solving the 2LSWEs with a bicharacteristics-theory-based scheme; their models are referred to as BD08 and DL13 models hereafter, respectively. A zoom-in view of Figure 7(c) in the region x∈[5m,6.5m] is provided in Figure 7(d). It is seen that the BD08 model predicts slight oscillations around the initial discontinuity location (marked out by the vertically dash–dotted line in Figure 7(c)). Note that the BD08 model result shown here is computed by using a rather small time step with Cr=0.01. Bouchut and de Luna [26] found that their model suffers high oscillations if Cr=0.5 (used in the present simulation) is adopted. The DL13 model also predicts small numerical wiggles, which do not appear in the present model results.Figure 7
Time evolution of the simulated profiles ofz1 (a) and z2 (b) for the 1D internal dam-break test. (c) shows the comparison of simulated interface profile with the numerical results of BD08 and DL13 models at t=1s. (d) is a zoom-in view of (c) in the region x∈[5m,6.5m]. The reference solution in (c) and (d) is obtained by the present model with nx=6400.
(a) (b) (c) (d)The order-of-accuracy of the proposed numerical scheme is studied here based on theL1 and L2 errors defined as(27)Lp=1n∑ifs-feip1/p∑ifeip1/p,where p=1,2; fs and fe denote the simulated and exact (or reference) solutions of variable f, respectively; n denotes the grid number. Table 1 presents the L1 and L2 errors and the orders of accuracy for the computed interface profile results at t=1s. It should be noted that the exact solution for this test is unavailable and we use the reference solution simulated by the present model with nx=6400 as the approximately exact solution. It is seen from Table 1 that generally the order of accuracy of the proposed scheme is about but slightly less than second-order. This is as expected because near the discontinuities the order-of-accuracy of the second-order MUSCL scheme is locally reduced owing to the use of a slope limiter (see (13a) and (13b)).Table 1
L 1 and L2 errors and the orders of accuracy for the computed interface profile results at t=1s for the 1D internal dam-break test.
nx L 1 error L 1 Order L 2 error L 2 Order n x = 200 1.00 E - 04  1.98 E - 04 n x = 400 2.81 E - 05 1.84 7.16 E - 05 1.46 n x = 800 7.12 E - 06 1.98 2.43 E - 05 1.56 n x = 1600 1.66 E - 06 2.10 7.92 E - 06 1.62 n x = 3200 2.97 E - 07 2.48 1.93 E - 06 2.04


### 3.3. 2D Internal Dam-Break Wave Propagations

Extending the assessment to 2D space, we consider here 2D internal dam-break problems in a domain(x,y)∈[-1m,1m]2. The tests are previously studied by Kurganov and Petrova [28] and Dudzinski and Lukác˘ová-Medvid’ová [17]. The initial conditions are prescribed as (28)z1x,y,t=0=0.50m,ifx,y∈Ω,0.45m,else,z2x,y,t=0=1.0m,u1x,y,t=0=u2x,y,t=0=2.5m/s,where Ω={(x+0.5)2+(y+0.5)2<0.5}∪{x<-0.5,y<0}∪{x<0,y<-0.5}. At all the simulation boundaries, an open boundary condition is applied.Two test cases are considered, with the first case occurring on a horizontally flat bed (zb(x,y)≡0) and the second one on a Gaussian-shaped bed topography defined as (29)zbx,y=0.05exp⁡-100x2+y2-1m.For both test cases, χ=0.98 and g=9.81m2/s.In Figures8(a) and 8(b), the simulated interface profiles at t=0.1s with nx=ny=400 are illustrated for the flat and uneven bed cases, respectively. As seen in Figure 8, the present numerical model captures the discontinuity sharply and the simulated results match well with those predicted by Kurganov and Petrova [28] and Dudzinski and Lukác˘ová-Medvid’ová [17], but generally with less oscillations. Because the boundary and initial conditions and also the bed topography are symmetrical about the diagonal line x=y, the numerical results should be symmetrical about this line, which are well predicted by the present model. In Figure 9, x=y slices of the results predicted with different grid resolutions are shown. It is clearly seen that as the mesh refines, the predicted interface sharpens and the multiwave structure is more obvious.Figure 8
Simulated interface profile att=0.1s with nx=ny=400 for the flat bed case (a), and the uneven bed case (b). The dashed lines represent the initial discontinuity location in the interface.
(a) (b)Figure 9
x = y slices of the computed interface profiles at t=0.1s on different meshes for the 2D internal dam-break flow over a flat bed (a), and an uneven bed (b). xD denotes the distance from the point (-1m,-1m) along the line x=y.
(a) (b)Table2 presents the L1 and L2 errors and the orders of accuracy for the computed interface profile results of flat and uneven bed cases shown in Figure 9 (note f=z1 and n=nx=ny in (27)). Note that the reference solutions simulated by the present model with nx=ny=3200 are used as the approximately exact solutions. It is seen from Table 2 that owing to the use of a slope limiter to enhance numerical stability, the order of accuracy of the proposed scheme is about but slightly less than second-order in this test (see also the results in the 1D internal dam-break test in Section 3.2).Table 2
L 1 and L2 errors and the orders of accuracy for the computed interface profile results at t=0.1s (c.f. Figure 9).
Case Numbers of cells L 1 error L 1 Order L 2 error L 2 Order Flat bed n x = n y = 100 6.00 E - 06  2.66 E - 05 n x = n y = 200 1.86 E - 06 1.69 9.88 E - 06 1.43 n x = n y = 400 7.34 E - 07 1.34 4.43 E - 06 1.16 n x = n y = 800 1.92 E - 07 1.94 1.32 E - 06 1.74 n x = n y = 1600 3.59 E - 08 2.42 2.82 E - 07 2.23  Uneven bed n x = n y = 100 1.67 E - 05  4.51 E - 05 n x = n y = 200 5.06 E - 06 1.72 1.70 E - 05 1.41 n x = n y = 400 1.31 E - 06 1.95 4.98 E - 06 1.77 n x = n y = 800 4.26 E - 07 1.63 1.52 E - 06 1.71 n x = n y = 1600 1.26 E - 07 1.75 4.51 E - 07 1.76


## 3.1. 1D Steady Flows over a Hump

1D two-layer shallow-water flows over an uneven bed topography are simulated here. Depending on the boundary conditions, the layer surface levels of the steady flows can be smooth, or nonsmooth with a hydraulic jump [12]. Three test cases are considered, including steady transitional flows with a hydraulic jump, steady transcritical exchange flows, and transcritical parallel flows with smooth layer surface levels. For all the test cases, approximate analytic solutions are available, which were derived by Abgrall and Karni [12] under a rigid-lid assumption. The computational channel is 6 m long (x∈[-3m,3m]) with the bed topography defined as (21)zbx=18cos⁡πx2+12,x⩽2m,0,else.Different boundary and initial conditions are applied. For the steady transitional flows with a hydraulic jump, (22)qx1x=-3m,t=0.0928m2/s,h1x=3m,t=0.9205m,qx2x=-3m,t=0.0928m2/s,h2x=3m,t=0.5795m,z1x,t=0=0.9205m,z2x,t=0=1.5m,and for the steady transcritical exchange flows, (23)qx1x=-3m,t=0.0928m2/s,h1x=3m,t=0.1617m,qx2x=3m,t=-0.0928m2/s,h2x=-3m,t=0.4311m,z1x,t=0=1.0m,z2x,t=0=1.5m,and for the steady transcritical parallel flows, (24)qx1x=-3m,t=0.0928m2/s,h1x=3m,t=0.1617m,qx2x=-3m,t=0.0928m2/s,h2x=3m,t=1.3338m,z1x,t=0=1.0m,z2x,t=0=1.5m.In all the test cases, χ=0.98 and g=10m2/s. For each simulation, when the globally relative change indictor R, defined below, drops smaller than 1×10-7, a steady solution is thought to be achieved. (25)R=∑k=12∑i=1nx∑j=1nyhki,jn+1-hki,jnhki,jn2.Here, nx and ny denote, respectively, the mesh cell numbers in the x- and y-directions, which are chosen to ensure numerical accuracy according to a grid-independence study. For instance, for the steady transitional flows with a hydraulic jump test case, a series of predictions at steady states on successively finer grids (i.e., nx=200,400,…,3200) are obtained. As seen from Figures 3(a) and 3(b) (note Figure 3(b) is a zoom-in view of Figure 3(a) around the hydraulic jump location), the simulated layer surface profiles with nx=1600 and nx=3200 are rather small and nx=1600 is chosen for this test case to achieve a balance between numerical accuracy and computational efficiency.Figure 3
Computed layer surface levels at steady states for the transitional flow with a hydraulic jump over a hump test: (a) whole simulation domain; (b) a zoom-in view in the regionx∈[0,1m]. (c) and (d) show the computed versus analytic profiles of interface elevation and layer discharges per unit width at steady states, respectively. The bed topography is plotted on the right y-axis.
(a) (b) (c) (d)In Figure3(c), the computed z1 at steady states is shown for the first test case. Also shown is the analytic solution of Abgrall and Karni [12], which is piecewise smooth with a hydraulic jump around x=0.48m. It is seen that the simulated results match the analytic solution fairly well, with slight deviations around the predicted hydraulic jump location. These deviations are thought to be ascribed to the limitation of the rigid-lid assumption adopted in deriving the analytic solution [12], and also the numerical error in predicting constant discharges in the layers [11]. In Figure 3(d), the computed layer discharges per unit width along with the analytic solution are displayed. A zoom-in view around the hydraulic jump location is given in Figure 3(d), from which it is seen that the predicted layer discharges agree with the analytic solution well, but deviations between the computed and analytic values exist near the jump location, with a maximum relative difference about 16%. This type of deviation in the predicted discharge has been commonly seen when an approximate Riemann solver is used in simulating steady one-layer or two-layer shallow-water flows over a hump involving a hydraulic jump (e.g., Zhou et al. [2]; Valiani and Begnudelli [25]; Lu et al. [11]). To test the stability of the numerical model and to check whether (20) is appropriate to dynamically determine the time step, a stability test is implemented here. For the first test case, the predicted interface profiles with nx=1600 using different values of Cr are shown in Figure 4. Note that all results in Figure 4 are plotted at steady flow states, except for the fact that in Figure 4(d) the results at t=1.0s are shown. It is found that when Cr<0.89 (see Figures 4(a)–4(c)), the numerical model is stable. When Cr⩾0.89, the numerical solution exhibits numerical instability as can be seen clearly in Figure 4(d). The stability study here indicates that the time step defined in (20) is suitable for dynamically determining the time step.Figure 4
Simulated interface profiles for the transitional flow with a hydraulic jump test withnx=1600 using Cr=0.3 (a), Cr=0.5 (b), Cr=0.8 (c), and Cr=0.89 (d). The results are shown at steady states, except in (d) where the result is illustrated at t=1.0s. A zoom-in view in the region x∈[-2.5,-1.5m] is shown in (d). The thin solid line at the bottom of each figure represents the bed profile and is plotted on the right y-axis.
(a) (b) (c) (d)Figures5(a) and 5(b) illustrate the computed layer surface levels and the discharges per unit width at steady states, along with the analytic solution of Abgrall and Karni [12] for the exchange flow test case. The corresponding results for the parallel flow test case are given in Figures 5(c) and 5(d). In the simulations, nx=200 is chosen according to grid-independence studies. From Figure 5, it is seen that the numerical results match the analytic solution reasonably well. According to Abgrall and Karni [12], the layer surface levels in these two test cases should be the same, which is well predicted by the developed model (c.f. Figures 5(a) and 5(c)).Figure 5
Computed (symbols) versus analytic (thick lines) layer surface levels (a) and discharges per unit width (b) at steady states, for the transcritical exchange flow over a hump test. The results for the parallel flow test case are shown in (c) and (d). The bed topography is shown at the bottom of each figure by the thin solid line and is plotted on the righty-axis.
(a) (b) (c) (d)At last in this test, we verify the performance of the developed model in predicting lake-at-rest solutions. The simulation setup is the same as the first test case except for the fact that at the western boundary the discharges in the layers are fixed to be zero. Figures6(a) and 6(b) present the computed versus analytic profiles of layer surface levels and discharges per unit width at t=10s, respectively. From Figure 6, it is seen that lake-at-rest solution is achieved.Figure 6
Simulated versus analytic layer surface levels (a) and discharges per unit width (b) att=10 s for the lake-at-rest test with χ≠1. The bed topography is plotted on the right y-axis.
(a) (b)


## 3.2. 1D Internal Dam-Break Wave Propagations

The 1D internal dam-break flow tests are widely used to investigate the robustness of numerical schemes in simulating flows involving large gradients (e.g., Bouchut and de Luna [26]; Dudzinski and Lukác˘ová-Medvid’ová [17]; Lee et al. [27]). The considered channel is 10 m long (x∈[0,10m]) and is horizontally flat (zb(x)≡0). The initial conditions are prescribed as (26)h1x,t=0=0.8sgn⁡x-5+1.0m,h2x,t=0=-h1x,t=0+2.0m,qx1x,t=0=qx2x,t=0=0.In the simulations, χ=0.7 and g=9.81m2/s. The simulation domain is discretized with nx=1600, chosen via a grid-independence study.Figures7(a) and 7(b) illustrate the time evolutions of the simulated interface and upper-layer surface, respectively. As seen, the present model is robust to capture sharp gradients. In Figure 7(c), the computed interface profile at t=1s is compared with those obtained by Bouchut and de Luna [26] via decoupled solving two 1LSWEs, and by Dudzinski and Lukác˘ová-Medvid’ová [17] solving the 2LSWEs with a bicharacteristics-theory-based scheme; their models are referred to as BD08 and DL13 models hereafter, respectively. A zoom-in view of Figure 7(c) in the region x∈[5m,6.5m] is provided in Figure 7(d). It is seen that the BD08 model predicts slight oscillations around the initial discontinuity location (marked out by the vertically dash–dotted line in Figure 7(c)). Note that the BD08 model result shown here is computed by using a rather small time step with Cr=0.01. Bouchut and de Luna [26] found that their model suffers high oscillations if Cr=0.5 (used in the present simulation) is adopted. The DL13 model also predicts small numerical wiggles, which do not appear in the present model results.Figure 7
Time evolution of the simulated profiles ofz1 (a) and z2 (b) for the 1D internal dam-break test. (c) shows the comparison of simulated interface profile with the numerical results of BD08 and DL13 models at t=1s. (d) is a zoom-in view of (c) in the region x∈[5m,6.5m]. The reference solution in (c) and (d) is obtained by the present model with nx=6400.
(a) (b) (c) (d)The order-of-accuracy of the proposed numerical scheme is studied here based on theL1 and L2 errors defined as(27)Lp=1n∑ifs-feip1/p∑ifeip1/p,where p=1,2; fs and fe denote the simulated and exact (or reference) solutions of variable f, respectively; n denotes the grid number. Table 1 presents the L1 and L2 errors and the orders of accuracy for the computed interface profile results at t=1s. It should be noted that the exact solution for this test is unavailable and we use the reference solution simulated by the present model with nx=6400 as the approximately exact solution. It is seen from Table 1 that generally the order of accuracy of the proposed scheme is about but slightly less than second-order. This is as expected because near the discontinuities the order-of-accuracy of the second-order MUSCL scheme is locally reduced owing to the use of a slope limiter (see (13a) and (13b)).Table 1
L 1 and L2 errors and the orders of accuracy for the computed interface profile results at t=1s for the 1D internal dam-break test.
nx L 1 error L 1 Order L 2 error L 2 Order n x = 200 1.00 E - 04  1.98 E - 04 n x = 400 2.81 E - 05 1.84 7.16 E - 05 1.46 n x = 800 7.12 E - 06 1.98 2.43 E - 05 1.56 n x = 1600 1.66 E - 06 2.10 7.92 E - 06 1.62 n x = 3200 2.97 E - 07 2.48 1.93 E - 06 2.04


## 3.3. 2D Internal Dam-Break Wave Propagations

Extending the assessment to 2D space, we consider here 2D internal dam-break problems in a domain(x,y)∈[-1m,1m]2. The tests are previously studied by Kurganov and Petrova [28] and Dudzinski and Lukác˘ová-Medvid’ová [17]. The initial conditions are prescribed as (28)z1x,y,t=0=0.50m,ifx,y∈Ω,0.45m,else,z2x,y,t=0=1.0m,u1x,y,t=0=u2x,y,t=0=2.5m/s,where Ω={(x+0.5)2+(y+0.5)2<0.5}∪{x<-0.5,y<0}∪{x<0,y<-0.5}. At all the simulation boundaries, an open boundary condition is applied.Two test cases are considered, with the first case occurring on a horizontally flat bed (zb(x,y)≡0) and the second one on a Gaussian-shaped bed topography defined as (29)zbx,y=0.05exp⁡-100x2+y2-1m.For both test cases, χ=0.98 and g=9.81m2/s.In Figures8(a) and 8(b), the simulated interface profiles at t=0.1s with nx=ny=400 are illustrated for the flat and uneven bed cases, respectively. As seen in Figure 8, the present numerical model captures the discontinuity sharply and the simulated results match well with those predicted by Kurganov and Petrova [28] and Dudzinski and Lukác˘ová-Medvid’ová [17], but generally with less oscillations. Because the boundary and initial conditions and also the bed topography are symmetrical about the diagonal line x=y, the numerical results should be symmetrical about this line, which are well predicted by the present model. In Figure 9, x=y slices of the results predicted with different grid resolutions are shown. It is clearly seen that as the mesh refines, the predicted interface sharpens and the multiwave structure is more obvious.Figure 8
Simulated interface profile att=0.1s with nx=ny=400 for the flat bed case (a), and the uneven bed case (b). The dashed lines represent the initial discontinuity location in the interface.
(a) (b)Figure 9
x = y slices of the computed interface profiles at t=0.1s on different meshes for the 2D internal dam-break flow over a flat bed (a), and an uneven bed (b). xD denotes the distance from the point (-1m,-1m) along the line x=y.
(a) (b)Table2 presents the L1 and L2 errors and the orders of accuracy for the computed interface profile results of flat and uneven bed cases shown in Figure 9 (note f=z1 and n=nx=ny in (27)). Note that the reference solutions simulated by the present model with nx=ny=3200 are used as the approximately exact solutions. It is seen from Table 2 that owing to the use of a slope limiter to enhance numerical stability, the order of accuracy of the proposed scheme is about but slightly less than second-order in this test (see also the results in the 1D internal dam-break test in Section 3.2).Table 2
L 1 and L2 errors and the orders of accuracy for the computed interface profile results at t=0.1s (c.f. Figure 9).
Case Numbers of cells L 1 error L 1 Order L 2 error L 2 Order Flat bed n x = n y = 100 6.00 E - 06  2.66 E - 05 n x = n y = 200 1.86 E - 06 1.69 9.88 E - 06 1.43 n x = n y = 400 7.34 E - 07 1.34 4.43 E - 06 1.16 n x = n y = 800 1.92 E - 07 1.94 1.32 E - 06 1.74 n x = n y = 1600 3.59 E - 08 2.42 2.82 E - 07 2.23  Uneven bed n x = n y = 100 1.67 E - 05  4.51 E - 05 n x = n y = 200 5.06 E - 06 1.72 1.70 E - 05 1.41 n x = n y = 400 1.31 E - 06 1.95 4.98 E - 06 1.77 n x = n y = 800 4.26 E - 07 1.63 1.52 E - 06 1.71 n x = n y = 1600 1.26 E - 07 1.75 4.51 E - 07 1.76


## 4. Discussions and Conclusions

In this work, a two-layer hydrostatic-reconstruction method (2LHRM) was proposed for solving the 2LSWEs. The 2LHMR are designed to satisfy various fundamental requirements, for example, the reconstructed “bed elevation” of the upper layer should be exactly the same with the surface level of the lower layer at the cell interfaces. A numerical model was developed in the framework of finite-volume methods to test the performance of the 2LHRM. A centered scheme, in particular, the SLIC scheme is chosen for numerical flux evaluations, making the proposed algorithms flexible for numerical implementations because no eigenstructure information of the 2LSWEs system is required in evaluating the numerical fluxes. The performances of the 2LHRM and the developed numerical model are assessed using various numerical tests, showing that the proposed method is capable of maintaining steady stationary flow at rest and captures spatially large gradients sharply. An extension of the proposed 2LHRM to handle wetting and drying processes and simulations of practical stratified flows would be interesting and are left for future researches.We should remark that our attention in this study is focused on stratified flows with different densities, i.e.,χ≠1. The relevant discussions in numerically predicting lake-at-rest solutions in Sections 2 and 3 are also for flows with χ≠1. In Appendix B, we tested the proposed method in predicting lake-at-rest solutions with χ=1, which shows that the proposed method is capable of predicting lake-at-rest solutions with χ=1, but with a smearing interface at locations where the reconstructed interface value of h2L is unequal to h2R. To improve the present two-layer reconstruction scheme to be able to predict lake-at-rest solutions while keeping the interface shape unchanged with χ=1 would be meaningful and needs future investigations.It is worth remarking that physically, Kelvin–Helmholtz instabilities may occur at the interface under some situations and the present model should also predict interfacial instabilities even if the eigenvalues are not explicitly computed. However, this is out of the scope of using the 2LSWEs and if one focuses on the interfacial instabilities, more complicated models (usually computationally more expensive), e.g., a multiphase Navier–Stokes equations based model (e.g., Hu et al. [21]; Yang et al. [29]) should be used.



---

*Source: 5064171-2019-07-07.xml*