Research Article  Open Access
Jiang Lin, Bing Mao, Xinhua Lu, "A TwoLayer HydrostaticReconstruction Method for HighResolution Solving of the TwoLayer ShallowWater Equations over Uneven Bed Topography", Mathematical Problems in Engineering, vol. 2019, Article ID 5064171, 14 pages, 2019. https://doi.org/10.1155/2019/5064171
A TwoLayer HydrostaticReconstruction Method for HighResolution Solving of the TwoLayer ShallowWater Equations over Uneven Bed Topography
Abstract
In this study, attention is focused on numerically solving the twodimensional, twolayer nonlinear shallowwater equations (2LSWEs) over uneven bed topography. A twolayer hydrostaticreconstruction method (2LHRM) is proposed for face value reconstructions. A numerical model is then developed based on the 2LHRM in the framework of finitevolume methods using the slopelimited 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 lakeatrest solutions without using any adhoc techniques.
1. Introduction
The last two decades have seen extensive use of upwind Godunovtype and centered finitevolume methods to solve the onelayer nonlinear shallowwater 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., Cproperty [1]), as well as to determine quantitative accuracy and model robustness in capturing spatially large gradients. It is found that a singlevalued definition of bed elevation at the cellinterfaces is convenient and of vital importance in devising a Cproperty satisfying numerical scheme [2–5]. In solving the 1LSWEs, the hydrostaticreconstruction method (HRM) proposed by Audusse et al. [3] is widely adopted; in this method, singlevalued intercell bed elevations are dynamically calculated based on the reconstructed intercell 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 wellbalanced 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 twolayer shallowwater equations (2LSWEs) to predict lakeatrest solutions, as well as to robustly capture spatially large gradients with highresolution. Lee et al. [9] employed the surfacegradient splitting method, originally proposed for solving the 1LSWEs [6], to obtain a wellbalanced 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 adhoc technique to switch from solving the waterlevel based continuity equation under a quiescent flow condition (QFC), to solving the waterdepth 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 adhoc 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 highresolution 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 griddependent. Izem et al. [13] developed a finiteelement 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 threelayer 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áováMedvid’ová [17] developed a finitevolume 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 Godunovtype method which relies on evaluation of eigenvalues of the partialdifferential 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 firstorder centered (FORCE) and its secondorder version, the slopelimited 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 twolayer hydrostaticreconstruction 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 collocatedgrid based finitevolume method. The performance of the 2LHRM and the developed model is then validated in terms of quantitative accuracy and Cproperty satisfying and shockcapturing properties.
The remainder of this paper is organized as follows. Section 2 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 twolayer shallowwater flow motions (see Figure 1). 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: with vectors defined by where denotes time; ) are the Cartesian coordinates with and ; ) are the layer depths with and denoting the values of the lower and upperlayers, respectively; and ) are the discharges per unit width in the  and directions, respectively; () denote the layeraveraged velocity components; is the bed elevation above a reference level (see Figure 1); ) denote the layer surface elevations with respect to and have the following relations: is the density ratio between upper and lower layers; and Note that the external shear stress (e.g., wind stress) at the upperlayer 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.
To solve the 2LSWEs defined in (1) and (2), the finitevolume 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 in this section and thus both the upperlayer and lowerlayer surfaces should be horizontally flat under a QFC. As for the QFC with , a discussion will be made in Appendix B.
When solving the 1LSWEs (obtained by setting in the 2LSWEs), it is required according to (3) that under general conditions (not only QFCs), e.g., in the direction, and under a QFC with , in order to achieve a lakeatrest solution, a reconstruction scheme should satisfy, e.g., in the direction, Eqs. (5) and (7) lead to meaning that the intercell bed elevation is singlevalued. One can easily prove that (5)–(7), together, ensure lakeatrest 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 lowerlayer variables, the reconstructions in relation to the upper layer should satisfy similar requirements as, e.g., in the direction, under general conditions, where denotes the reconstructed intercell “bed elevation” for the upper layer, and under a QFC with a constant upperlayer surface elevation , In addition to the constraints given in (5)–(11), physically, it is also required that under general conditions, 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, and , which are the water surface elevation and water depth in the 1LSWEs, respectively, are firstly reconstructed (e.g., by a firstorder piecewise constant or piecewise linear method). Then, based on the reconstructed values, a singlevalued intercell bed elevation is defined to satisfy (8), and after that, and 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 twolayer hydrostaticreconstruction method (2LHRM), reconstructs the face values through the following successive procedures.
Step 1. Reconstruction of the upperlayer variables: in this step, the face values of , , , and are interpolated from those at the cell centers, e.g., in the direction, by using the monotonic upstreamcentered scheme for conservation laws (MUSCL) [22] with the minmod slope limiter [23] for numerical stabilitywhere can be any of the variables to be reconstructed; , with , and . Note that when setting , the firstorder piecewise constant reconstruction scheme is recovered.
Step 2. Reconstruction of the lowerlayer variables: in this step, the face values of , , and are reconstructed from (13a) and (13b) (with ). Then, the face value of is computed from the relation (see (3)). After that, we set to satisfy (12).
Step 3. Calculation of cellinterface velocities by
Step 4. Calculation of the singlevalued intercell as, e.g., in the direction, Here, is used. Once is computed, is redefined as Then, and are successively updated from (5) and (12), respectively.
Step 5. Because may be modified in Step 4, we need to redefine as Then, is recalculated from (9).
Step 6. Computation of layer discharges at the cell interface by where and are computed in Step 3. and are calculated according to (4) as, e.g.,
It is obvious that the reconstructed face values by using the 2LHRM satisfy all the relations in (5)–(12). As a result, lakeatrest solutions are ensured because under a QFC, the numerical fluxes are zero since in (A.3)–(A.5) and , and the source terms in (A.6a), (A.6b), and (A.6c) are also zero because under a QFC the upper and lowerlayer surfaces are horizontally flat and thus (see (A.7)). It is also meaningful to point out that the 2LHRM reduces to the HRM when the upperlayer does not exist (i.e., and ).
To gain a better understanding of the twolayer reconstruction approach, we illustrate the major reconstruction procedures in Figure 2. Without loss of generality, we assume that after Step 3, . According to the magnitudes of and obtained after Step 3 and (abbreviated as hereafter) obtained in Step 4, three typical cases exist, i.e., , , and (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, and do not change in all three cases; does not change except when ; does not change except when ; and change in all cases due to the singlevalued definition of . Note that for the tests studied later in Section 3, only the first case will be encountered.
(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: where and are approximately local bounds of wave speeds [10, 11], and 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 finitevolume method. In the model, the slopelimited 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 dambreak waves propagating over a flat bottom, and 2D internal dambreak waves propagating over a flat and an uneven bed topography. The first test is chosen to check whether the proposed method predicts lakeatrest solutions and the capability of the model in simulating twolayer shallowwater 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 orderofaccuracy of the developed numerical model. The coefficient 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 direction is set to be unity.
3.1. 1D Steady Flows over a Hump
1D twolayer shallowwater 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 rigidlid assumption. The computational channel is 6 m long () with the bed topography defined as Different boundary and initial conditions are applied. For the steady transitional flows with a hydraulic jump, and for the steady transcritical exchange flows, and for the steady transcritical parallel flows, In all the test cases, and . For each simulation, when the globally relative change indictor , defined below, drops smaller than , a steady solution is thought to be achieved. Here, and denote, respectively, the mesh cell numbers in the  and directions, which are chosen to ensure numerical accuracy according to a gridindependence 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., ) are obtained. As seen from Figures 3(a) and 3(b) (note Figure 3(b) is a zoomin view of Figure 3(a) around the hydraulic jump location), the simulated layer surface profiles with and are rather small and is chosen for this test case to achieve a balance between numerical accuracy and computational efficiency.
(a)
(b)
(c)
(d)
In Figure 3(c), the computed 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 . 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 rigidlid 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 zoomin 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 . This type of deviation in the predicted discharge has been commonly seen when an approximate Riemann solver is used in simulating steady onelayer or twolayer shallowwater 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 using different values of 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 are shown. It is found that when (see Figures 4(a)–4(c)), the numerical model is stable. When , 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.
(a)
(b)
(c)
(d)
Figures 5(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, is chosen according to gridindependence 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)).
(a)
(b)
(c)
(d)
At last in this test, we verify the performance of the developed model in predicting lakeatrest 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. Figures 6(a) and 6(b) present the computed versus analytic profiles of layer surface levels and discharges per unit width at , respectively. From Figure 6, it is seen that lakeatrest solution is achieved.
(a)
(b)
3.2. 1D Internal DamBreak Wave Propagations
The 1D internal dambreak 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áováMedvid’ová [17]; Lee et al. [27]). The considered channel is 10 m long () and is horizontally flat (). The initial conditions are prescribed as In the simulations, and . The simulation domain is discretized with , chosen via a gridindependence study.
Figures 7(a) and 7(b) illustrate the time evolutions of the simulated interface and upperlayer surface, respectively. As seen, the present model is robust to capture sharp gradients. In Figure 7(c), the computed interface profile at is compared with those obtained by Bouchut and de Luna [26] via decoupled solving two 1LSWEs, and by Dudzinski and LukáováMedvid’ová [17] solving the 2LSWEs with a bicharacteristicstheorybased scheme; their models are referred to as BD08 and DL13 models hereafter, respectively. A zoomin view of Figure 7(c) in the region 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 . Bouchut and de Luna [26] found that their model suffers high oscillations if (used in the present simulation) is adopted. The DL13 model also predicts small numerical wiggles, which do not appear in the present model results.
(a)
(b)
(c)
(d)
The orderofaccuracy of the proposed numerical scheme is studied here based on the and errors defined aswhere ; and denote the simulated and exact (or reference) solutions of variable , respectively; denotes the grid number. Table 1 presents the and errors and the orders of accuracy for the computed interface profile results at . 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 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 secondorder. This is as expected because near the discontinuities the orderofaccuracy of the secondorder MUSCL scheme is locally reduced owing to the use of a slope limiter (see (13a) and (13b)).

3.3. 2D Internal DamBreak Wave Propagations
Extending the assessment to 2D space, we consider here 2D internal dambreak problems in a domain . The tests are previously studied by Kurganov and Petrova [28] and Dudzinski and LukáováMedvid’ová [17]. The initial conditions are prescribed as where . 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 () and the second one on a Gaussianshaped bed topography defined as For both test cases, and .
In Figures 8(a) and 8(b), the simulated interface profiles at with 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á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 , the numerical results should be symmetrical about this line, which are well predicted by the present model. In Figure 9, 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.
(a)
(b)
(a)
(b)
Table 2 presents the and errors and the orders of accuracy for the computed interface profile results of flat and uneven bed cases shown in Figure 9 (note and in (27)). Note that the reference solutions simulated by the present model with 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 secondorder in this test (see also the results in the 1D internal dambreak test in Section 3.2).

4. Discussions and Conclusions
In this work, a twolayer hydrostaticreconstruction 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 finitevolume 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., . The relevant discussions in numerically predicting lakeatrest solutions in Sections 2 and 3 are also for flows with . In Appendix B, we tested the proposed method in predicting lakeatrest solutions with , which shows that the proposed method is capable of predicting lakeatrest solutions with , but with a smearing interface at locations where the reconstructed interface value of is unequal to . To improve the present twolayer reconstruction scheme to be able to predict lakeatrest solutions while keeping the interface shape unchanged with 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.
Appendix
A. Numerical Methods
In the framework of FVMs, (1) may be explicitly discretized on a uniform Cartesian grid as Here, and denote cell indexes; is the time step; and denote the space intervals in the  and directions, respectively; superscripts “” and “” denote the values at the old and new time levels, respectively; and are the numerical flux vectors; the tilde on a term denotes that this term is evaluated by the temporally evolved values described below in Section Appendix A.1.
A.1. Numerical Flux Evaluation
To calculate the numerical fluxes in (A.1), the SLIC scheme is employed. This scheme involves three steps to compute a numerical flux. First, the face values are reconstructed from those at the cell centers by using (13a) and (13b). Then, the face values are evolved over a half timestep as where the superscripts “L”/“R” and “B”/“T” denote the reconstructed values in the first step at the left/right side of the cell interface in the direction, and bottom/top side of the cell interface in the direction, respectively. The vector is discretized by using the method described later in Section Appendix A.2. Finally, the numerical flux vectors in (A.1) are computed by arithmetically averaging the twostep Lax–Wendroff flux and the Lax–Friedrichs flux as [19], e.g., in the direction, where
A.2. Discretizations of
The discretization of is involved in calculating in (A.1), and in (A.2c). Taking in (A.2c) as an example, we discretize it as follows:where Here, denotes the th component of . and are discretized similarly as and , respectively. Note that the method of using (A.7) to discretize the gradient of layer surface elevation was originally proposed by Ying and Wang [30] in solving the 1LSWEs.
The vector in (A.1) is discretized in the same manner as that in (A.2c), but calculated using the values being temporally evolved (see (A.2a), (A.2b), (A.2c)).
B. Property in Predicting LakeatRest Solutions with = 1
For the twolayer shallowwater flows, a QFC is different whether is equal to 1 or not. When , the initially quiescent twolayer flows can be maintained at rest only when both the upper and lowerlayer surface elevations are constants (for wet regions). However, when , it is only required that the upperlayer surface elevation is constant, while the interface profile can be of arbitrary shape. Few previous studies have considered a QFC with , but numerical schemes are expected to maintain stationary flow at rest when . Here, we consider a 1D twolayer lakeatrest test with . The simulation domain is 5 m long () with the bed topography defined as Initially, and is defined as The flows are initially at rest. The initial layer surface levels along with the bed terrain are shown in Figure 10(a). As defined in (B.1) and (B.2), there is a discontinuity in the bed terrain and the initial interface profile at . There is another discontinuity in the bed profile at , where the initial layer surface elevations are constants. The simulation domain is discretized uniformly with . Open boundary conditions are prescribed at both ends of the simulation domain. It is obvious that the initially flows should be maintained at rest for this test.
(a)
(b)
The simulated time evolution of the profiles of the layer surfaces and discharges per unit width are illustrated in Figures 10(a) and 10(b), respectively. As can be seen in Figure 10(b), no spurious flow is predicted in both layers, which can be explained simply as follows. When solving the momentum equations in (1) and (2) (i.e., the second, third, fifth, and sixth components of (1) and (2)), the numerical fluxes defined in (A.3) are zero because , where denotes the th component of and ; and the corresponding source terms (see (A.6b) and (A.6c)) are also zero because (1) for a horizontally flat upperlayer surface one obtains from (A.7), and thus the first terms on the righthandside (R.H.S.) of (A.6b) and (A.6c) are zero; (2) the second term on the R.H.S. of (A.6c) is also zero because . Therefore, from (A.1) we have , and thus no spurious flow is predicted when solving the momentum equations. Nevertheless, the predicted interface smears around as can be seen in Figure 10(a). This kind of “diffusion” in the interface is not caused by the discontinuity in the bed profile because the interface shape does not change with time at where the bed profile also contains a discontinuity (see Figure 10(a)). When we look at the continuity equations, i.e., the first and fourth components of (1) and (2), no source terms exist for these two equations. As for the numerical mass fluxes, it is zero for the fourth component of (1) because and (see (A.3)–(A.5)), while the numerical mass flux related to the first component of (1) is nonzero where the reconstructed values satisfy , i.e., . At locations where differs with , the second term on the R.H.S. of (A.5) is nonzero when solving the first component of (1) (note the other terms on the R.H.S. of (A.4a) and (A.5) are zero), and this causes a nonzero numerical mass flux which leads to a numerical solution with smearing interface. It should be noted that the numerically predicted time variation in the interface does not cause spurious flows in later times. It should be also noted that as the interface smears with time, the maximum difference between and and thus the second term on the R.H.S. of (A.5) and the numerical mass flux defined in (A.3) gradually decreases. As seen in Figure 10(a), the difference in the interface profiles at and is indistinguishable. From this test, we see that the proposed 2LHRM is capable of predicting lakeatrest solutions with , but with a smearing interface at locations where the reconstructed layer depth . The numerically predicted smearing interface may not be seen as an error for this particular test because actually this is a onelayer flow test and the interface does not physically exist.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported in part by National Natural Science Foundation of China (Grant No. 51409195), the Central Level, Scientific Research Institutes for Basic R&D Special Fund Business of Yangtze River Scientific Research Institute (Grant No. CKSF2016013/HL), and the Natural Science Foundation of Hubei Province of China (Grant No. 2016CFB385).
References
 A. Berbudez and M. E. Vazquez, “Upwind methods for hyperbolic conservation laws with source terms,” Computers & Fluids, vol. 23, no. 8, pp. 1049–1071, 1994. View at: Publisher Site  Google Scholar
 J. G. Zhou, D. M. Causon, C. G. Mingham, and D. M. Ingram, “The surface gradient method for the treatment of source terms in the shallowwater equations,” Journal of Computational Physics, vol. 168, no. 1, pp. 1–25, 2001. View at: Publisher Site  Google Scholar  MathSciNet
 E. Audusse, F. Bouchut, M.O. Bristeau, R. Klein, and B. Perthame, “A fast and stable wellbalanced scheme with hydrostatic reconstruction for shallow water flows,” SIAM Journal on Scientific Computing, vol. 25, no. 6, pp. 2050–2065, 2004. View at: Publisher Site  Google Scholar  MathSciNet
 X. Lu and S. Xie, “Depthaveraged nonhydrostatic numerical modeling of nearshore wave propagations based on the FORCE scheme,” Coastal Engineering, vol. 114, pp. 208–219, 2016. View at: Publisher Site  Google Scholar
 X. Lu, B. Mao, and B. Dong, “Comment on “An efficient and stable hydrodynamic model with novel source term discretization schemes for overland flow and flood simulations” by Xilin Xia et al,” Water Resources Research, vol. 54, no. 1, pp. 621–627, 2018. View at: Publisher Site  Google Scholar
 Q. Liang and A. G. L. Borthwick, “Adaptive quadtree simulation of shallow flows with wetdry fronts over complex topography,” Computers & Fluids, vol. 38, no. 2, pp. 221–234, 2009. View at: Publisher Site  Google Scholar
 J. Hou, Q. Liang, F. Simons, and R. Hinkelmann, “A stable 2D unstructured shallow flow model for simulations of wetting and drying over rough terrains,” Computers & Fluids, vol. 82, pp. 132–147, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 X. Lu, B. Dong, B. Mao, and X. Zhang, “A twodimensional depthintegrated nonhydrostatic numerical model for nearshore wave propagation,” Ocean Modelling, vol. 96, pp. 187–202, 2015. View at: Publisher Site  Google Scholar
 W. K. Lee, A. G. L. Borthwick, and P. H. Taylor, “On mathematical balancing of a twolayer shallow flow model,” in Proceddings of the 1st IAHR European IAHR Congress, Edinburgh, UK, 2010. View at: Google Scholar
 B. Spinewine, V. Guinot, S. SoaresFrazão, and Y. Zech, “Solution properties and approximate Riemann solvers for twolayer shallow flow models,” Computers & Fluids, vol. 44, no. 1, pp. 202–220, 2011. View at: Publisher Site  Google Scholar
 X. Lu, B. Dong, B. Mao, and X. Zhang, “A robust and wellbalanced numerical model for solving the twolayer shallow water equations over uneven topography,” Comptes Rendus Mecanique, vol. 343, no. 78, pp. 429–442, 2015. View at: Publisher Site  Google Scholar
 R. Abgrall and S. Karni, “Twolayer shallow water system: a relaxation approach,” SIAM Journal on Scientific Computing, vol. 31, no. 3, pp. 1603–1627, 2009. View at: Publisher Site  Google Scholar
 N. Izem, M. Seaid, I. Elmahi, and M. Wakrim, “Discontinuous Galerkin method for twodimensional bilayer shallow water equations,” Journal of Engineering Mathematics, vol. 96, no. 1, pp. 1–21, 2016. View at: Publisher Site  Google Scholar
 A. Chertock, A. Kurganov, A. Kurganov, Z. Qu, and T. Wu, “Threelayer approximation of twolayer shallow water equations,” Mathematical Modelling and Analysis, vol. 18, no. 5, pp. 675–693, 2013. View at: Publisher Site  Google Scholar
 J. T. Frings, An Adaptive Multilayer Model for DensityLayered Shallow Water Flows [Ph.D. thesis], RWTH Aachen University, 2012.
 M. J. Castro, E. D. Fernández, J. M. González, and C. Parés, “Numerical treatment of the loss of hyperbolicity of the twolayer shallowwater system,” Journal of Scientific Computing, vol. 48, no. 13, pp. 16–40, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 M. Dudzinski and M. LukáčováMedvid’ová, “Wellbalanced bicharacteristicbased scheme for multilayer shallow water flows including wet/dry fronts,” Journal of Computational Physics, vol. 235, pp. 82–113, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 K. Arun, M. Kraft, M. LukáčováMedvid’ová, and P. Prasad, “Finite volume evolution Galerkin method for hyperbolic conservation laws with spatially varying flux functions,” Journal of Computational Physics, vol. 228, no. 2, pp. 565–590, 2009. View at: Publisher Site  Google Scholar
 E. F. Toro, ShockCapturing Methods for FreeSurface Shallow Flows, Wiley, 2001.
 C. Swartenbroekx, B. Spinewine, V. Guinot, and Y. Zech, “Hyperbolicity preserving HLL solver for twolayer shallowwater equations applied to dambreak flows,” in Proceedings of the 5th International Conference on Fluvial Hydraulics, pp. 1379–1387, Karlsruhe, Germany, 2010. View at: Google Scholar
 Y. Hu, X. Guo, X. Lu, Y. Liu, R. A. Dalrymple, and L. Shen, “Idealized numerical simulation of breaking water wave propagating over a viscous mud layer,” Physics of Fluids, vol. 24, no. 11, Article ID 112104, 2012. View at: Publisher Site  Google Scholar
 B. Van Leer, “Towards the ultimate conservative difference scheme. V. A secondorder sequel to Godunov’s method,” Journal of Computational Physics, vol. 32, pp. 101–136, 1979. View at: Publisher Site  Google Scholar  MathSciNet
 P. L. Roe and M. J. Baines, “Algorithms for advection and shock problems,” in Proceedings of the 4th GAMM Conference on Numerical Methods in Fluid Mechanics, pp. 281–290, Vieweg, Braunschweig, Germany, 1982. View at: Google Scholar
 B. F. Sanders, “Highresolution and nonoscillatory solution of the St. Venant equations in nonrectangular and nonprismatic channels,” Journal of Hydraulic Research, vol. 39, no. 3, pp. 321–330, 2010. View at: Publisher Site  Google Scholar
 A. Valiani and L. Begnudelli, “Divergence form for bed slope source term in shallow water equations,” Journal of Hydraulic Engineering, vol. 132, no. 7, pp. 652–665, 2006. View at: Publisher Site  Google Scholar
 F. Bouchut and T. Morales de Luna, “An entropy satisfying scheme for twolayer shallow water equations with uncoupled treatment,” ESAIM: Mathematical Modelling and Numerical Analysis, vol. 42, no. 4, pp. 683–698, 2008. View at: Publisher Site  Google Scholar
 W. Lee, A. G. Borthwick, and P. H. Taylor, “A fast adaptive quadtree scheme for a twolayer shallow water model,” Journal of Computational Physics, vol. 230, no. 12, pp. 4848–4870, 2011. View at: Publisher Site  Google Scholar
 A. Kurganov and G. Petrova, “Centralupwind schemes for twolayer shallow water equations,” SIAM Journal on Scientific Computing, vol. 31, no. 3, pp. 1742–1773, 2009. View at: Publisher Site  Google Scholar
 Z. Yang, X.H. Lu, X. Guo, Y. Liu, and L. Shen, “Numerical simulation of sediment suspension and transport under plunging breaking waves,” Computers & Fluids, vol. 158, pp. 57–71, 2017. View at: Publisher Site  Google Scholar  MathSciNet
 X. Ying and S. S. Wang, “Improved implementation of the HLL approximate Riemann solver for onedimensional open channel flows,” Journal of Hydraulic Research, vol. 46, no. 1, pp. 21–34, 2010. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Jiang Lin et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.