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.
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á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. 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 two-layer shallow-water 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 upper-layers, respectively; and ) are the discharges per unit width in the - and -directions, respectively; () denote the layer-averaged 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 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.
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 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 , 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 lake-at-rest solution, a reconstruction scheme should satisfy, e.g., in the -direction, Eqs. (5) and (7) lead to 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 -direction, under general conditions, where denotes the reconstructed inter-cell “bed elevation” for the upper layer, and under a QFC with a constant upper-layer 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 first-order piecewise constant or piecewise linear method). Then, based on the reconstructed values, a single-valued inter-cell 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 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 of , , , and are interpolated from those at the cell centers, e.g., in the -direction, by using the monotonic upstream-centered 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 first-order piecewise constant reconstruction scheme is recovered.
Step 2. Reconstruction of the lower-layer 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 cell-interface velocities by
Step 4. Calculation of the single-valued inter-cell 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, lake-at-rest 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 lower-layer 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 upper-layer does not exist (i.e., and ).
To gain a better understanding of the two-layer 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 single-valued 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 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 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 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 () 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 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., ) 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 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 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 . 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 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 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)).
(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. 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 lake-at-rest solution is achieved.
(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á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 grid-independence study.
Figures 7(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 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 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 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 order-of-accuracy 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 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)).
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 . 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 Gaussian-shaped 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 second-order in this test (see also the results in the 1D internal dam-break test in Section 3.2).
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., . The relevant discussions in numerically predicting lake-at-rest solutions in Sections 2 and 3 are also for flows with . In Appendix B, we tested the proposed method in predicting lake-at-rest solutions with , which shows that the proposed method is capable of predicting lake-at-rest solutions with , but with a smearing interface at locations where the reconstructed interface value of is unequal to . To improve the present two-layer reconstruction scheme to be able to predict lake-at-rest 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 time-step 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 two-step 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 Lake-at-Rest Solutions with = 1
For the two-layer shallow-water flows, a QFC is different whether is equal to 1 or not. When , the initially quiescent two-layer flows can be maintained at rest only when both the upper- and lower-layer surface elevations are constants (for wet regions). However, when , it is only required that the upper-layer 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 two-layer lake-at-rest 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 upper-layer surface one obtains from (A.7), and thus the first terms on the right-hand-side (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 lake-at-rest 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 one-layer 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).