# A Finite Volume Method for Modeling Shallow Flows with Wet-Dry Fronts on Adaptive Cartesian Grids

Authors: Sheng Bi; Jianzhong Zhou; Yi Liu; Lixiang Song
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A second-order accurate, Godunov-type upwind finite volume method on dynamic refinement grids is developed in this paper for solving shallow-water equations. The advantage of this grid system is that no data structure is needed to store the neighbor information, since neighbors are directly specified by simple algebraic relationships. The key ingredient of the scheme is the use of the prebalanced shallow-water equations together with a simple but effective method to track the wet/dry fronts. In addition, a second-order spatial accuracy in space and time is achieved using a two-step unsplit MUSCL-Hancock method and a weighted surface-depth gradient method (WSDM) which considers the local Froude number is proposed for water depths reconstruction. The friction terms are solved by a semi-implicit scheme that can effectively prevent computational instability from small depths and does not invert the direction of velocity components. Several benchmark tests and a dam-break flooding simulation over real topography cases are used for model testing and validation. Results show that the proposed model is accurate and robust and has advantages when it is applied to simulate flow with local complex topographic features or flow conditions and thus has bright prospects of field-scale application.


---

## Body

## 1. Introduction

Solution to the two-dimensional (2D) hyperbolic conservation laws of the shallow-water equations (SWE) is relevant to many real-life hydrodynamic problems, such as free surface flows in shallow lakes, dam breaks, dyke breaches, wide rivers, flash floods, tidal bores, and coastal inundation, among others. Meanwhile, a variety of methods can be used to solve the nonlinear shallow-water equations, such as variational principle [1, 2], finite element method [3], and finite volume method (FVM) [4–8]. HE [1] deduced a semi-inverse method for establishment of a generalized variational principle which applies to any conservation equations of hydrodynamics, and this method is helpful to study the analytical solutions which are useful in shallow-water equations. In recent years, with the development of numerical methods, more and more attentions have been paid to the FVM which are based on the conservative form of the governing equations. In particular, the upwind Godunov-type schemes via the FVM are very popular because of their robustness and accuracy in capturing discontinuities.When solving the shallow-water equations in the context of a Godunov-type framework, it is generally necessary and efficient to use the well-balanced methods, which employ a special treatment for the bed slope source terms that balanced the source terms and flux gradients. This well-balanced concept is also known by exact conservation property (C-property) [4]. Several authors have conducted this study; see, for example, [4, 6, 7]. Among these, the prebalanced SWEs derived by choosing the water level and discharge as dependent variables are derived by Liang and Borthwick [6]. However, near the wet/dry fronts, negative water depth would be produced which could cause the mass error. To deal with this problem, Liang and Marche [8] then proposed the nonnegative reconstruction of Riemann states and compatible discretization of the slope source term to maintain stable and well-balanced property for dry-bed simulations. Begnudelli and Sanders [9] adopted flux correction terms to preserve the well-balanced property in a fully wet cell; however, the correction terms were later proved to be a potential source of large unphysical velocity [10].Moreover, it is important that the numerical shallow-flow models are able to simulate accurately the evolution of relatively small-scale features, such as fronts, within large complex domains. Therefore, various approaches to adaptive mesh resolution for improving computational efficiency have been suggested for the shallow-water equations in the published literatures. Hu et al. [11] briefly reviewed the different mesh patching methods and presented a new patching technique for solving the SWEs using a finite volume Godunov-type scheme. Rogers et al. [12] and Liang and Borthwick [6] presented typical adaptive quadtree solvers of the shallow-water equations to achieve local grid refinement. Sleigh et al. [13] developed a finite volume method model based on adaptive unstructured triangular grids, which are widely used in the domains with highly complicated boundary configurations due to their capability of robustness. However, it is noteworthy that most of the above adaptive grids generation techniques require data structure to store their neighbour information that increases the computational cost. All in all, it is reasonable to anticipate that the ideally grid generators should be robust and efficient and fit the flow boundaries correctly.Various numerical methods developed for general systems of hyperbolic conservation laws have been applied to achieve second-order accuracy in space, on structured and unstructured meshes. These methods provide efficient approaches to accommodate sub-, super-, and transcritical flows as well as flows with discontinuities. The conventional depth gradient method (DGM) approach which evaluates intercell water depths starting from the extrapolation of the same conserved variable can obtain a more robust behavior in the front tracking. However, if a centered discretization is used for the bed slope source terms, the depth gradient method may not maintain the static condition. Zhou et al. [7] proposed surface gradient method (SGM), which is suitable for both steady and unsteady shallow-water flows and is easy to implement. A central-upwind scheme using the surface elevation instead of the water depth has been used by Kurganov and Levy [14]. The SGM is efficient, but near the wet/dry fronts on uneven topographies the SGM reconstruction can lead to very small or negative predicted water depths that may cause numerical instabilities. Therefore, there is a need to modify the numerical scheme in order to avoid unphysical results.It is significant to design a numerical approach to handle wetting and drying, which is one of the central challenges that researchers are tackling. In Godunov-type finite volume schemes, the depth-averaged velocity components are normally determined by dividing the discharge per unit width by the local water depth which can inevitably lead to predictions of unacceptable negative water depth and numerical instability. Liang and Borthwick [6] presented a Godunov-type shallow-flow model to solve the well-balanced SWEs with wet/dry fronts by a local bed slope modification method. Brufau et al. [15] presented a numerical technique that temporarily modifies ground elevation to approximate wetting and drying for both steady and unsteady flows. Begnudelli and Sanders [9] tracked fluid volume and the free surface elevation in partially submerged cells for wetting and drying and developed a high-resolution, unstructured finite volume model for shallow-water flows over arbitrary topography. The fully implicit treatments of friction source terms are employed by Yoon and Kang [16] to alleviate the problems associated with numerical instabilities due to small water depths near the wet/dry fronts. Moreover, near the wet/dry interface, where the depth is vanishing, the friction terms may cause the stiff problem and significantly affect the stability of the numerical model [17]. Therefore, it is necessary to employ a robust scheme, such as the semi-implicit scheme [9, 17], for friction terms when modeling the real-world shallow flow over irregular terrain.In this paper, we describe details of a computationally efficient numerical approach to simulating complex shallow flow. A structured, but nonuniform, Cartesian grid system which allows a simulation to be carried out in dynamic refinement grids is used for solving the SWEs. A weighted surface-depth gradient method (WSDGM) is adopted with the aim of combining the advantages of SGM and DGM approaches and avoiding their drawbacks. An efficient and robust technique is used to track the stationary or move the wet/dry fronts. The friction terms are treated by a semi-implicit scheme, which may relieve the problem of high velocity when the water depth is small and does not invert the direction of velocity components and thus would benefit model robustness. This paper is organized as follows. In Section2, the prebalanced shallow-water equation is derived. Section 3 presents the nonuniform gridding technique. Section 4 describes the numerical solution to the equations on dynamic refinement grids using a second-order Godunov-type finite volume scheme. Finally, the results of model validation tests and a simulation of a flood wave are drawn in Section 5.


## 2. Shallow-Water Equations

The two-dimensional depth-integrated shallow-water equations are obtained by integrating the Navier-Stokes equations over the flow depth with the following assumptions: hydrostatic pressure distribution, small bottom slope, uniform velocity distribution in the vertical direction, incompressible fluid, and constant fluid density with free surface. In matrix form, the conservation laws of the two-dimensional nonlinear shallow-water equations are(1)∂U∂t+∂E∂x+∂G∂y=S,
where t is the time; x and y are Cartesian coordinates; U, E, G, and S are vectors of the conserved flow variables, fluxes in the x- and y-directions, and source terms, respectively. Traditionally, the vector terms are given by
(2)U=[hhuhv],E=[huhu2+gh22huv],G=[hvhuvhv2+gh22],S=[0-ghSfx-gh∂zb∂x-ghSfy-gh∂zb∂y],
where h is the water depth; u and v are the depth-averaged velocity components in the x- and y-directions, respectively; g=9.81m/s2 is the gravity acceleration; ∂zb/∂x and ∂zb/∂y are bed slopes in the x- and y-directions, respectively; Sfx and Sfy are friction slopes in the x- and y-directions, respectively. In this study, the friction slopes are estimated by using the Manning formulae
(3)Sfx=n2uu2+v2h4/3,Sfx=n2vu2+v2h4/3,
in which n is Manning’s roughness coefficient. In cases involving wetting and drying, any grid cell with h<10-6m is assumed to be dry.Note thatη=h+zb, in which η is defined as the water surface elevation above datum. Then, -(g/2)(∂h2/∂x)-gh(∂zb/∂x) can be rewritten as -(g/2)(∂(η-zb)2/∂x)-g(η-zb)(∂zb/∂x). Simple algebraic manipulation leads to
(4)-g2∂h2∂x-gh∂zb∂x=-g2∂(η2-2ηzb)∂x-gη∂zb∂x.
A similar analysis can be applied in the y-direction momentum equation. Based on this, (2) will be transformed into another formulation constituted by (1) and
(5)U=[ηhuhv],E=[huhu2+g(η2-2ηzb)2huv],G=[hvhuvhv2+g(η2-2ηzb)2],S=S0+Sf=[0-gη∂zb∂x-gη∂zb∂y]+[0-ghSfx-ghSfy].Consider thex-direction momentum equation at the initial steady state in which η is a constant for a wet-bed application with h≠0. After eliminating terms with zero velocities, the momentum equation reduces to
(6)g2∂(η2-2ηzb)∂x=-gη∂zb∂x.
When using HLLC approximate Riemann solver based Godunov-type numerical scheme to calculate the interface fluxes, the term on the left-hand side of (6) is approximated by (fE-fW)/Δx, where fE and fW are fluxes through the east and west cell interfaces, respectively. Since η is a constant and h≠0, it is easy to prove that no matter which calculation formula in HLLC solver is used, the fluxes are
(7)fE=g2(η2-2ηzbE),fW=g2(η2-2ηzbW).
Based on (7), the flux gradient is
(8)fE-fWΔx=-gηzbE-ηzbWΔx.
If the bed slope term is discretized as -(ηzbE-ηzbW)/Δx, (6) is balanced. A similar analysis can be applied in the y-direction momentum equation. Therefore, (1) and (5) are in flux gradient and source terms balanced form without any need to use special numerical treatment for the source terms.


## 3. Computational Grid

### 3.1. Structured but Nonuniform Grids

Because of the complexity of computational domain, a simple structured rectangular mesh requires a large number of cells to resolve the detailed flow pattern near the dam breaking location, navigation channels, and in-stream structures. To optimize the use of computational resources, we use the multiple-level nonuniform rectangular mesh with local refinement proposed by Liang [18].In this mesh system, various levels of fine cells are placed close to the dam-break locations and in-stream structures or inlet areas where the flow gradients are high, while coarse grids are used in the low-gradient regions. To simplify the mesh, a cell is refined by splitting into four equal child cells. Corresponding to this refining, the computational mesh will be simpler because any cell has one or two faces on each of its neighbors, as shown in Figure1.Figure 1
Typical structured but nonuniform grids.To obtain a structured but nonuniform grid is straightforward and only involves the following four simple steps.(1)
Discretize the computational domain using a coarse uniform grid.(2)
Divide the unit square into four equal quadrant cells according to specific subdivision criteria.(3)
Check each cell in turn and subdivide it if necessary.(4)
Regularize the grid to ensure that an arbitrary cell is at most two times bigger or smaller than any of its neighbors (2 : 1 rule).The first step is simply to determine the number of cells in thex- and y-directions, respectively. In the second step of subdividing the unit uniform grid (root grid) into four quadrants (child cells), seeding points are adopted or certain criteria may be used in this work. Then, each cell (including child cells) is checked in turn and a prescribed highest subdivision level is allocated to those cells containing seeding points. The procedure essentially creates an irregular grid, where a cell may have a neighbor that is many times bigger or smaller than its own. The stability and accuracy of the solution may be affected if no action is taken. In order to reduce grid irregularity, the 2 : 1 rule is satisfied and provides a trade-off between grid flexibility and accuracy.Figure1 shows the different levels of cells (i-1,j), (i,j), and (i+1,j), where i and j represent the cell indexes in the x- and y-directions, respectively. Considering cell (i,j), as Figure 1 shows, the subgrid is indexed by is=1,Ms and js=1,Ms, where Ms is the number of subgrid cells in the x- or y- direction. Ms=2lev with lev being the subdivision level of cell (i,j). The advantage of this nonuniform grid system is that no data structure is needed to store the neighbor information, neighbors are directly specified by simple algebraic relationships, and details on implementing the neighbor interpolation are well documented in the literature [18].


### 3.2. Dynamic Grid Adaptation

In this work, the grid adaptation criterion is based on the averaged water surface gradient and the critical value ofϕsub and ϕcoar should be prescribed for grid refinement and coarsening. The averaged water surface gradient at cell ic is defined as [6]
(9)ϕic=(∂η∂x)ic2+(∂η∂y)ic2.
The following procedures are implemented for grid refinement and coarsening.(1)
Check every cell in turn; if any cell hasϕic>ϕsub and the subdivision level is less than the maximum, the cell will be marked for further subdivision and its subdivision level will be increased tolev +1. When a cell is marked for subdivision, its eight neighbors are also checked and subdivided if necessary in order to meet the 2 : 1 rule.(2)
For grid coarsening, if the cell hasϕic<ϕcoar and the subdivision level is greater than 0, the cell will be marked for coarsening and its subdivision level will be changed tolev −1. Again the 2 : 1 rule must remain satisfied after grid coarsening.It is worth noting that the above adaptation procedure is performed at every time step during a simulation, and the flow information of a newly created subcell is obtained by the aforementioned interpolation method.


## 3.1. Structured but Nonuniform Grids

Because of the complexity of computational domain, a simple structured rectangular mesh requires a large number of cells to resolve the detailed flow pattern near the dam breaking location, navigation channels, and in-stream structures. To optimize the use of computational resources, we use the multiple-level nonuniform rectangular mesh with local refinement proposed by Liang [18].In this mesh system, various levels of fine cells are placed close to the dam-break locations and in-stream structures or inlet areas where the flow gradients are high, while coarse grids are used in the low-gradient regions. To simplify the mesh, a cell is refined by splitting into four equal child cells. Corresponding to this refining, the computational mesh will be simpler because any cell has one or two faces on each of its neighbors, as shown in Figure1.Figure 1
Typical structured but nonuniform grids.To obtain a structured but nonuniform grid is straightforward and only involves the following four simple steps.(1)
Discretize the computational domain using a coarse uniform grid.(2)
Divide the unit square into four equal quadrant cells according to specific subdivision criteria.(3)
Check each cell in turn and subdivide it if necessary.(4)
Regularize the grid to ensure that an arbitrary cell is at most two times bigger or smaller than any of its neighbors (2 : 1 rule).The first step is simply to determine the number of cells in thex- and y-directions, respectively. In the second step of subdividing the unit uniform grid (root grid) into four quadrants (child cells), seeding points are adopted or certain criteria may be used in this work. Then, each cell (including child cells) is checked in turn and a prescribed highest subdivision level is allocated to those cells containing seeding points. The procedure essentially creates an irregular grid, where a cell may have a neighbor that is many times bigger or smaller than its own. The stability and accuracy of the solution may be affected if no action is taken. In order to reduce grid irregularity, the 2 : 1 rule is satisfied and provides a trade-off between grid flexibility and accuracy.Figure1 shows the different levels of cells (i-1,j), (i,j), and (i+1,j), where i and j represent the cell indexes in the x- and y-directions, respectively. Considering cell (i,j), as Figure 1 shows, the subgrid is indexed by is=1,Ms and js=1,Ms, where Ms is the number of subgrid cells in the x- or y- direction. Ms=2lev with lev being the subdivision level of cell (i,j). The advantage of this nonuniform grid system is that no data structure is needed to store the neighbor information, neighbors are directly specified by simple algebraic relationships, and details on implementing the neighbor interpolation are well documented in the literature [18].


## 3.2. Dynamic Grid Adaptation

In this work, the grid adaptation criterion is based on the averaged water surface gradient and the critical value ofϕsub and ϕcoar should be prescribed for grid refinement and coarsening. The averaged water surface gradient at cell ic is defined as [6]
(9)ϕic=(∂η∂x)ic2+(∂η∂y)ic2.
The following procedures are implemented for grid refinement and coarsening.(1)
Check every cell in turn; if any cell hasϕic>ϕsub and the subdivision level is less than the maximum, the cell will be marked for further subdivision and its subdivision level will be increased tolev +1. When a cell is marked for subdivision, its eight neighbors are also checked and subdivided if necessary in order to meet the 2 : 1 rule.(2)
For grid coarsening, if the cell hasϕic<ϕcoar and the subdivision level is greater than 0, the cell will be marked for coarsening and its subdivision level will be changed tolev −1. Again the 2 : 1 rule must remain satisfied after grid coarsening.It is worth noting that the above adaptation procedure is performed at every time step during a simulation, and the flow information of a newly created subcell is obtained by the aforementioned interpolation method.


## 4. The Numerical Method

### 4.1. Finite Volume Method

Integrating (1) over an arbitrary cell, the integral form thereof is
(10)∂∂t∫UdΩ+∫(∂F∂x+∂G∂y)dΩ=∫ΩSdΩ,
in which all the vector terms are given by (5). By applying Green’s theorem to the second term in (10),
(11)∂∂t∫ΩUdΩ+∫sF·ndS=∫ΩSdΩ,
where S is the boundary of Ω and F is the vector of flux functions through S given by
(12)F·n=Fnx+Gny,
in which nx and ny are the Cartesian components of the unit normal vector n. The FVM is based on (11), the discretization of which ensures that mass and momentum are conserved across a discontinuity. Discretizing (11), the equation becomes
(13)Uin+1=Uin-ΔtΔx(FE-FW)-ΔtΔy(GN-GS)+ΔtSi,
where superscript n is the time level, subscript i represents the cell index, Δt is the time step, FE and FW are the fluxes through the east and west cell interfaces, respectively, and GN and GS are the fluxes through the south and north cell interfaces, respectively.


### 4.2. HLLC Approximate Riemann Solver

To capture shock waves, reduce the numerical damping effect and overcome the numerical oscillation problems; the Harten-Lax-van Leer-Contact (HLLC) Riemann solver is employed to solve the local Riemann problems and calculate the interface fluxes. An interface flux, such asFE, is computed [6]:
(14)FE={FLifs1≥0,F*Lifs1<0≤s2,F*Rifs2<0≤s3,FRifs3≤0,
where FL=F(UL) and FR=F(UR) are calculated from the left and right Riemann states UL and UR, respectively, for a local Riemann problem and F*L and F*R are the numerical fluxes in the left and right sides of the middle region of the Riemann solution, respectively, given by
(15)F*L=[f*1f*2vLf*1],F*R=[f*1f*2vRf*1],
where vL and vR represent the left and right tangential velocity components of the Riemann states, respectively, and f*1 and f*2 are the first two components of the normal flux F*, which is calculated using the HLLC solver by
(16)F*=s3FL-s1FR+s1s3(UR-UL)s3-s1,
where s1, s2, and s3 represent the speeds of the left, contact, and right waves, respectively, which are determined by (17). The third flux component in the middle flux vector F* is calculated by f*3=vLf*1 or vRf*1:
(17)s1={uR-2ghRifhL=0,min⁡(uL-ghL,u*-gh*)ifhL>0,s3={uL+2ghLifhR=0,max⁡(uR+ghR,u*+gh*)ifhR>0,s2=s1hR(uR-s3)-s3hL(uL-s1)hR(uR-s3)-hL(uL-s1),
where hL, hR, uL, and uR are the components of the left and right initial Riemann states for a local Riemann problem, with the middle states h* and u* calculated by
(18)u*=12(uL+uR)+ghL-ghR,h*=1g[12(ghL+ghR)+14(uL-uR)]2.
Similar formulae are used to calculate FW, GN, and GS.


### 4.3. High-Order Reconstruction and WSDGM

To achieve higher-order accuracy, which is the most important aspect for any flow solver, better reconstruction techniques than a piecewise constant function should always be used. Overall second-order accuracy in space and time is achieved using a two-step unsplit MUSCL-Hancock method. The predictor and corrector steps are given by (19) and (20), respectively,
(19)Uin+1/2=Uin-Δt2Δx(FE-FW)-Δt2Δy(GN-GS)+Δt2Si,(20)Uin+1=Uin+1/2-ΔtΔx(FE-FW)-ΔtΔy(GN-GS)+ΔtSi.
No Riemann solutions are required at the predictor step and the interface fluxes are essentially evaluated within each cell using interpolated values at the extremities of the inner interfaces [6].In the corrector step, the fluxes through the cell interfaces are evaluated by the HLLC approximate Riemann solver. Then, the Riemann states can be reconstructed using the WSDGM, taking the eastern interfaceUi+1/2,j as an example:
(21)Ui+1/2,jL=Ui,jn+1/2+12ϕ(r)(Ui,jn-Ui-1,jn),Ui+1/2,jR=Ui+1,jn+1/2-12ϕ(r)(Ui+1,jn-Ui,jn),
where Ui,jn+1/2 and Ui+1,jn+1/2 are the predicted cell-centre values at cell Ci and its eastern neighbor obtained in the predictor step; the slope limiter ϕ(r) is minmod limiter function;r is the ratio of gradients given by
(22)r=Ui+1,jn-Ui,jnUi,jn-Ui-1,jn.
Using WSDGM, the first component of (21) is reconstructed as
(23)hi+1/2,jL=wi,jhi+1/2,jL,D+(1-wi,j)hi+1/2,jL,S,hi+1/2,jL=wi,jhi+1/2,jL,D+(1-wi,j)hi+1/2,jL,S,
where wi,j is a weighting parameter which is defined as
(24)wi,j={00≤Fri,j≤11Fri,j>1,
and hi+1/2,jL,D, hi+1/2,jL,S, hi+1/2,jR,D, and hi+1/2,jR,S are the DGM and SGM reconstructions, respectively. They are labeled as
(25)hi+1/2,jL,D=hi,jn+12ϕ(r)(hi,jn-hi-1,jn),hi+1/2,jR,D=hi+1,jn-12ϕ(r)(hi+1,jn-hi,jn),hi+1/2,jL,S=ηi+1/2,jL-zbi+1/2,j=ηi,jn+12ϕ(r)(ηi,jn-ηi-1,jn)-zbi+1/2,j,hi+1/2,jR,S=hi+1/2,jR,S-zbi+1/2,j=ηi+1,jn-12ϕ(r)(ηi+1,jn-ηi,jn)-zbi+1/2,j,
where bi+1/2,j=max⁡(ηi+1/2,jL-hi+1/2,jL,ηi+1/2,jR-hi+1/2,jR) represents an estimation of the intercell bed elevation if a linear variation of the bed profile is assumed.


### 4.4. Friction Terms Treatments

When considering a case with complex topography and involving wetting and drying, the friction term actually dominates the stability of a numerical model. In particular, when the water depth becomes very small, the Manning formulation of friction may lead to an exaggerated force that can even reverse the flow, which is obviously physically incorrect. To deal with this problem, in this work, the friction terms are solved by the semi-implicit scheme which is given by(26)un+1=11+Δtτnun,vn+1=11+Δtτnvn,
where τ=gn2(un)2+(vn)2(hn)-4/3; un,vn, and hn are the initial states for the friction source term treatment. This scheme does not invert the direction of velocity and prevent instability effectively because (26) means un+1·un≥0 and vn+1·vn≥0. Moreover, when the water depth becomes very small, this method may relieve the problem of high velocity. For example, when hn=10-2, un=500m/s, vn=-800m/s, and n=0.018m/s1/3, then the velocities computed by (26) are un+1=3.56m/s and vn+1=-5.7m/s. This indicates that the semi-implicit scheme (26) can well benefit the computational stability.


### 4.5. Time Stepping and Wet/Dry Treatments

The numerical model is explicit, and its stability is controlled by the Courant-Friedrichs-Lewy (CFL) condition. For a two-dimensional Cartesian grid, the CFL criterion is used for predicting an appropriate time step [6]:
(27)Δt=Cmin⁡(Δtx,Δty),
with Δtx=min⁡i(Δxi/(|ui|+ghi)) and Δty=min⁡i(Δyi/(|vi|+ghi)), where 0≤C≤1 is the Courant number and C = 0.8 is adopted in this study. It is worth mentioning that the local time stepping method [18] can be implemented to achieve higher computational efficiency.When considering the application to wetting and drying processes on a complex terrain, there is a simple but effective method to capture wet/dry fronts. After the whole solution is updated by (20), the computed water depths lower than zero are set at zero; besides, the velocities of the dry cell and its neighbors are also set at zero. This simple method which prevents the fake fluxes through the dry cell interfaces could sharpen the wet/dry fronts of the new conserved flow variables Un+1.


## 4.1. Finite Volume Method

Integrating (1) over an arbitrary cell, the integral form thereof is
(10)∂∂t∫UdΩ+∫(∂F∂x+∂G∂y)dΩ=∫ΩSdΩ,
in which all the vector terms are given by (5). By applying Green’s theorem to the second term in (10),
(11)∂∂t∫ΩUdΩ+∫sF·ndS=∫ΩSdΩ,
where S is the boundary of Ω and F is the vector of flux functions through S given by
(12)F·n=Fnx+Gny,
in which nx and ny are the Cartesian components of the unit normal vector n. The FVM is based on (11), the discretization of which ensures that mass and momentum are conserved across a discontinuity. Discretizing (11), the equation becomes
(13)Uin+1=Uin-ΔtΔx(FE-FW)-ΔtΔy(GN-GS)+ΔtSi,
where superscript n is the time level, subscript i represents the cell index, Δt is the time step, FE and FW are the fluxes through the east and west cell interfaces, respectively, and GN and GS are the fluxes through the south and north cell interfaces, respectively.


## 4.2. HLLC Approximate Riemann Solver

To capture shock waves, reduce the numerical damping effect and overcome the numerical oscillation problems; the Harten-Lax-van Leer-Contact (HLLC) Riemann solver is employed to solve the local Riemann problems and calculate the interface fluxes. An interface flux, such asFE, is computed [6]:
(14)FE={FLifs1≥0,F*Lifs1<0≤s2,F*Rifs2<0≤s3,FRifs3≤0,
where FL=F(UL) and FR=F(UR) are calculated from the left and right Riemann states UL and UR, respectively, for a local Riemann problem and F*L and F*R are the numerical fluxes in the left and right sides of the middle region of the Riemann solution, respectively, given by
(15)F*L=[f*1f*2vLf*1],F*R=[f*1f*2vRf*1],
where vL and vR represent the left and right tangential velocity components of the Riemann states, respectively, and f*1 and f*2 are the first two components of the normal flux F*, which is calculated using the HLLC solver by
(16)F*=s3FL-s1FR+s1s3(UR-UL)s3-s1,
where s1, s2, and s3 represent the speeds of the left, contact, and right waves, respectively, which are determined by (17). The third flux component in the middle flux vector F* is calculated by f*3=vLf*1 or vRf*1:
(17)s1={uR-2ghRifhL=0,min⁡(uL-ghL,u*-gh*)ifhL>0,s3={uL+2ghLifhR=0,max⁡(uR+ghR,u*+gh*)ifhR>0,s2=s1hR(uR-s3)-s3hL(uL-s1)hR(uR-s3)-hL(uL-s1),
where hL, hR, uL, and uR are the components of the left and right initial Riemann states for a local Riemann problem, with the middle states h* and u* calculated by
(18)u*=12(uL+uR)+ghL-ghR,h*=1g[12(ghL+ghR)+14(uL-uR)]2.
Similar formulae are used to calculate FW, GN, and GS.


## 4.3. High-Order Reconstruction and WSDGM

To achieve higher-order accuracy, which is the most important aspect for any flow solver, better reconstruction techniques than a piecewise constant function should always be used. Overall second-order accuracy in space and time is achieved using a two-step unsplit MUSCL-Hancock method. The predictor and corrector steps are given by (19) and (20), respectively,
(19)Uin+1/2=Uin-Δt2Δx(FE-FW)-Δt2Δy(GN-GS)+Δt2Si,(20)Uin+1=Uin+1/2-ΔtΔx(FE-FW)-ΔtΔy(GN-GS)+ΔtSi.
No Riemann solutions are required at the predictor step and the interface fluxes are essentially evaluated within each cell using interpolated values at the extremities of the inner interfaces [6].In the corrector step, the fluxes through the cell interfaces are evaluated by the HLLC approximate Riemann solver. Then, the Riemann states can be reconstructed using the WSDGM, taking the eastern interfaceUi+1/2,j as an example:
(21)Ui+1/2,jL=Ui,jn+1/2+12ϕ(r)(Ui,jn-Ui-1,jn),Ui+1/2,jR=Ui+1,jn+1/2-12ϕ(r)(Ui+1,jn-Ui,jn),
where Ui,jn+1/2 and Ui+1,jn+1/2 are the predicted cell-centre values at cell Ci and its eastern neighbor obtained in the predictor step; the slope limiter ϕ(r) is minmod limiter function;r is the ratio of gradients given by
(22)r=Ui+1,jn-Ui,jnUi,jn-Ui-1,jn.
Using WSDGM, the first component of (21) is reconstructed as
(23)hi+1/2,jL=wi,jhi+1/2,jL,D+(1-wi,j)hi+1/2,jL,S,hi+1/2,jL=wi,jhi+1/2,jL,D+(1-wi,j)hi+1/2,jL,S,
where wi,j is a weighting parameter which is defined as
(24)wi,j={00≤Fri,j≤11Fri,j>1,
and hi+1/2,jL,D, hi+1/2,jL,S, hi+1/2,jR,D, and hi+1/2,jR,S are the DGM and SGM reconstructions, respectively. They are labeled as
(25)hi+1/2,jL,D=hi,jn+12ϕ(r)(hi,jn-hi-1,jn),hi+1/2,jR,D=hi+1,jn-12ϕ(r)(hi+1,jn-hi,jn),hi+1/2,jL,S=ηi+1/2,jL-zbi+1/2,j=ηi,jn+12ϕ(r)(ηi,jn-ηi-1,jn)-zbi+1/2,j,hi+1/2,jR,S=hi+1/2,jR,S-zbi+1/2,j=ηi+1,jn-12ϕ(r)(ηi+1,jn-ηi,jn)-zbi+1/2,j,
where bi+1/2,j=max⁡(ηi+1/2,jL-hi+1/2,jL,ηi+1/2,jR-hi+1/2,jR) represents an estimation of the intercell bed elevation if a linear variation of the bed profile is assumed.


## 4.4. Friction Terms Treatments

When considering a case with complex topography and involving wetting and drying, the friction term actually dominates the stability of a numerical model. In particular, when the water depth becomes very small, the Manning formulation of friction may lead to an exaggerated force that can even reverse the flow, which is obviously physically incorrect. To deal with this problem, in this work, the friction terms are solved by the semi-implicit scheme which is given by(26)un+1=11+Δtτnun,vn+1=11+Δtτnvn,
where τ=gn2(un)2+(vn)2(hn)-4/3; un,vn, and hn are the initial states for the friction source term treatment. This scheme does not invert the direction of velocity and prevent instability effectively because (26) means un+1·un≥0 and vn+1·vn≥0. Moreover, when the water depth becomes very small, this method may relieve the problem of high velocity. For example, when hn=10-2, un=500m/s, vn=-800m/s, and n=0.018m/s1/3, then the velocities computed by (26) are un+1=3.56m/s and vn+1=-5.7m/s. This indicates that the semi-implicit scheme (26) can well benefit the computational stability.


## 4.5. Time Stepping and Wet/Dry Treatments

The numerical model is explicit, and its stability is controlled by the Courant-Friedrichs-Lewy (CFL) condition. For a two-dimensional Cartesian grid, the CFL criterion is used for predicting an appropriate time step [6]:
(27)Δt=Cmin⁡(Δtx,Δty),
with Δtx=min⁡i(Δxi/(|ui|+ghi)) and Δty=min⁡i(Δyi/(|vi|+ghi)), where 0≤C≤1 is the Courant number and C = 0.8 is adopted in this study. It is worth mentioning that the local time stepping method [18] can be implemented to achieve higher computational efficiency.When considering the application to wetting and drying processes on a complex terrain, there is a simple but effective method to capture wet/dry fronts. After the whole solution is updated by (20), the computed water depths lower than zero are set at zero; besides, the velocities of the dry cell and its neighbors are also set at zero. This simple method which prevents the fake fluxes through the dry cell interfaces could sharpen the wet/dry fronts of the new conserved flow variables Un+1.


## 5. Numerical Results and Discussion

In all the simulations undertaken in this section,g=9.81m/s2 and ρ=1000kg/m3.


### 5.1. Still Water Test with Wet/Dry Front

This test is chosen to verify that the present model indeed preserves the stationary steady flow with wet/dry fronts over an uneven bottom. The numerical experiment takes place in a square pool1m×1m. At the centre of the pool, a symmetric bump at the bottom is mathematically defined as the form
(28)z(x,y)=max⁡[0,0.25-5*((x-0.5)2+(y-0.5)2)].The still water as the initial condition covering partially the bump in a frictionless domain is fully closed by vertical walls; the water surface elevation is 0.1 m and 0.3 m. A simulation oft=200 s is performed on a grid of 50 × 50 cells. Figure 2 shows the computed 3D water surface in the different conditions. Similarly, Figure 3 presents the water level in cells versus the x-coordinate of the cell centroid along the line y=0.5 m. In terms of [19, 20], the errors,
(29)L1(η)=1Nc∑i=1Nc|ηi-ηi¯|,L1(q)=1Nc∑i=1Nc|qi-qi¯|,
where q=hu2+v2, ηi¯, and qi¯ denote the analytical solutions and Nc is the number of cells with positive water depth, are listed in Table 1. The errors are sufficiently small to the level of round-off error of a double accuracy calculation. In addition, ε1=maxt((Vt-V0+ΔVt)/V0)=7.9×10-16, ε2=0.0 represents the maximum relative error in global mass conservation over the simulation duration under the two conditions, in which V0 and Vt are the initial volume of water and the computed volume of water at the time t, respectively, and ΔV represents the net volume transported out of the domain from the beginning of the simulation period to the time t.Table 1
Still water test: error norms for free surface elevation and unit discharge.
Cases
L
1
(
η
)
L
1
(
q
)
Water level = 0.1 m
5.43
×
1
0
-
16 m
6.0
×
1
0
-
16 m2/s
Water level = 0.3 m
1.22
×
1
0
-
18 m
2.2
×
1
0
-
17 m2/sStill water with wet/dry front: the 3D water surface aftert=300 s at different cases.
(a)
(b)Still water with wet/dry front: the water surface elevation along the liney=0.5 m at different cases. The results are obtained at t=300 s.
(a)
(b)The results show that the stationary solutions are well maintained, and the current numerical scheme is able to preserve the lake at rest solution involving wet/dry interface and hence is well balanced.


### 5.2. Moving Shorelines in a Parabolic Bowl with Friction

Thacker [21] provides a good benchmark test case for the present shallow-flow model in dealing with wetting and drying and bed slope as well as the friction source term. This paper firstly deduced the analytical solution of the test case.Consider the case where the motion of shallow water in a basin is governed by the equations [21]
(30)∂ζ∂t+∂(zb+ζ)∂x+∂(zb+ζ)∂y=0,(31)∂u∂t+u∂u∂x+v∂u∂y+τu+g∂ζ∂x=0,(32)∂v∂t+u∂v∂x+v∂v∂y+τv+g∂ζ∂y=0,
where ζ(x,y,t) is the height of the water surface above mean water level, zb is the bottom surface, u(x,y,t) is the averaged component of the water current to the east, v(x,y,t) is the depth-averaged velocity component of the water current to the north, τ is the bottom friction parameter which is considered to be a constant (related to the bed friction coefficient by Cf=hτu2+v2), g is the acceleration due to gravity, and t is the time.Assuming thatu was function of t only, u=u0(t), and v=0, then (31) and (32) imply that
(33)ζ(x,y,t)=ζx=0(t)+xζ1(t),(34)ζ1(t)=-1g(du0(t)dt+τu0(t)).The bed profile of the domain is defined by(35)zb(x)=h0(xa)2,
where h0 and a are constants so that flow takes place above parabolic bottom topography.Substituting (33), (34), (35), u=u0(t), and v=0 in (30) gives
(36)dζ0(t)dt+xdζ1(t)dt-2u0(t)h0xa2+u0(t)ζ1(t)=0.
Following Thacker, equate the time-varying coefficients of the linearly independent terms, 1 and x, respectively,
(37)dζ0(t)dt+u0(t)ζ1(t)=0,(38)xdζ1(t)dt-2u0(t)h0xa2=0.
Substituting (34) in (38),
(39)d2u0(t)dt2+τdu0(t)dt+2gh0u0(t)a2=0.
Substituting (34) in (37),
(40)dζ0(t)dt-1gu0(t)du0dt-τgu0(t)2=0.
The auxiliary equation for (39) is
(41)λ2+τλ+2gh0a2=0.
The roots of this equation are
(42)λ=-τ±τ2-p22,p=8gh0a2.In this test case, the computational domain has planned dimensions of-5000m×5000m; the coefficients are a=3000m, h0=10m, and τ=0.001s-1. So, considering the case with τ<p, a solution of (39) is
(43)u0(t)=Be-τt/2sinst,s=p2-τ22,
where B=5m/s is a constant, obtained by using given values for u0(t) and u0′(t). Substituting (43) in (40) and integrating with respect to t give(44)ζ0(t)=a2B2e-τt8g2h0(-sτsin2st+(τ24-s2)cos⁡2st)-B2e-τt4g.
Substituting (43) in (34) gives
(45)ζ1(t)=-e-τt/2g(Bscos⁡st+τB2sinst).
Substituting (44) and (45) in (33) gives
(46)ζ(t)=-e-τt/2g(Bscos⁡st+τB2sinst)x+a2B2e-τt8g2h0[-sτsin2st+(τ24-s2)cos⁡2st]-B2e-τt4g.
Thus the analytical water level η(x,t)=max⁡(h0+ζ,zb(x)) is given by
(47)η(x,t)=max⁡([-sτsin2st+(τ24-s2)cos⁡2st]h0-e-τt/2g(Bscos⁡st+τB2sinst)xmmax+a2B2e-τt8g2h0[-sτsin2st+(τ24-s2)cos⁡2st]mmax-B2e-τt4g,zb(x)[-sτsin2st+(τ24-s2)cos⁡2st]).We compared the numerical solution of the model in this paper and the analytical solution; the simulation lasted 6000 s. Figure4 shows a good agreement between the numerical predictions and analytical solutions of the free surface profile at different times; the moving shoreline is accurately simulated with no signs of amplified spurious oscillations and due to the friction effect, the sloshing amplitude of water level is damping with increasing simulation time. Figure 5 shows the comparison between numerical and analytical solutions for water depth at points −2750 m, 50 m, and 2750 m. The numerical predictions are found to match perfectly with those given by the analytical solutions and no obvious distortion is detected at the wet-dry interface. The results validate the wetting and drying method of the present model, as well as the source terms treatments.Sloshing motions in a vessel with parabolic bottom topography at different output times: (a)t=1000 s; (b) t=2000 s; (c) t=3000 s; (d) t=4000 s; (e) t=5000 s; (f) t=6000 s.
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
Moving shorelines in a parabolic bowl with friction: comparison between numerical and analytical results for water depth at points −2750 m, 50 m, and 2750 m.


### 5.3. 1D Steady Flows over a Steep Bump with WSDGM

The goal of this test case is to investigate the ability of the present model to converge to a nonstationary steady state. The test concerns 1D steady flows in a[-10m≤x≤10m], which is discretized with 0.1 m grid spacing to make the test more severe, frictionless channel. The bottom profile, characterized by the presence of a steep bump, is described by the following equation [22]:
(48)zb(x)={0.8(1-x24)-2m≤x≤2m0elsewhere.In order to demonstrate the capability of the WSDGM reconstruction, the results of supercritical flow, subcritical flow, and transcritical flow with a shock are displayed in Figures6, 7, and 8. The three cases were simulated using a mesh of 400 (Δx=0.05m) cells to steady states.1D steady flow over a steep bump-supercritical flow: comparison between reference solutions and numerical results.
(a)
(b)
(c)1D steady flow over a steep bump-subcritical flow: comparison between reference solutions and numerical results.
(a)
(b)
(c)1D steady flow over a steep bump-transcritical flow with a shock: (a) water surface elevation and (b) discharge.
(a)
(b)For the case of supercritical flow over the bump, we impose a constant discharge per unit widthq=1.5m2/s,h=0.25m at the upstream boundary and the downstream boundary set to be transmissive [22]; the Fr range from 2.31 to 3.83. Figure 6, in the subdomain [-3m≤x≤3m], shows the numerical results of the DGM (Frlim⁡=0, (a)), SGM (Frlim⁡→∞, (b)), and WSDGM (Frlim⁡=1, (c)) reconstruction, respectively. Figure 6(b) shows that the numerical scheme with Frlim⁡→∞ converges toward the steady solution worse than Frlim⁡=1. The results verified that pure SGM reconstruction is unsuitable when the flow is supercritical.On the other hand, if the flow is subcritical (e.g.,q=1.0m2/s at the upstream and h=1.7 at the downstream), then the Fr range from 0.14 to 0.41; Figure 7(a) shows a dip in the water elevation upstream and downstream of the bump, where the bottom profile is not smooth. Moreover, the pure DGM reconstruction does not satisfy the C-property.In Figures6(c), 7(c), and 8, excellent agreements are observed between the numerical predictions and analytical solutions under different conditions. All in all, the results prove that WSDGM performs better than the schemes based on pure DGM or SGM reconstructions. For this reason, the parameter Frlim⁡ was set at 1 also in the following tests.


### 5.4. 2D Circular Dam-Break Problem

This problem is based on the hypothetical test case studied by Alcrudo and Garcia-Navarro [23], which involves the breaking of a circular dam. In this problem, a circular discontinuity over a 200 m × 200 m is investigated. The circle’s radius is 50 m and it is centered at the middle of the square. The initial states are 10 m water depth inside the circle, and the water is at rest. The initial water depth outside the circle is considered for two cases, that is, 5 m and dry bed. The wall is then assumed to be removed completely and no slope or friction is considered. It is assumed that the circular dam is broken instantaneously and water surface profile, 2 s after dam failure, is studied.In these examples, the computational grid consists of 100 × 100 cells with each cell 4 m2 as the background grid. During the simulation, the adaptation criteria ϕsub=0.08 and ϕcoar=0.05 are used. Figure 9 shows the computed water surface profile 2 s after dam failure outside water depth equal to 5 m and dry, respectively. Obviously, the circular dam-break waves spread and propagate radically and symmetrically as the water drains from the deepest region. Figure 10 presents the final numerical outputs at 2 s after the dam failed in terms of adapted grid. From the 3D water surface and depth contours, it is conspicuous that the sharp shock front is well captured and the depth contours show that most of the depth contour lines appear to be circular. The case outside the water depth is equal to 5 m, and the flow regime changes from subcritical to supercritical. This case can be a good test for validation of a numerical scheme. In the case of dry-bed simulation, it can be seen that no bore forms, but instead a rarefaction wave extends into the dry region (Figure 11). The scheme is capable of handling the dry-bed problem. Similar results were obtained by Alcrudo and Garcia-Navarro [23] and Zoppou and Roberts [5].Circular dam break: (a) 3D water surface; (b) depth contours for the circular dam break att=2 s with the downstream water depth set at 5 m.
(a)
(b)Figure 10
2D circular dam break: adapted grid att=2 s.Circular dam break: (a) 3D water surface; (b) depth contours att=2 s with the downstream water depth set dry.
(a)
(b)In order to perform the advantage of the adaptive grid-based model, the simulations are also run on uniform grids with different resolutions. The results are shown in Figure12 in terms of water surface profile. The structured adaptive grid-based prediction matches closely the one produced on the finest uniform grid with 200 × 200 cells, which indicates similar accuracy of the two simulations. On the adaptive grid system, the mesh is refined to those areas with sharp surface gradients. Therefore, the computational efficiency is greatly improved due to the fact that the total number of cells is reduced. As listed in Table 2, the new adaptive grid can save more than twice of the computational cost, compared with the simulation on a uniform grid with 200 × 200 cells.Table 2
Circular dam break: CPU time required by different simulations.
Grid
Uniform
Adaptive
Number of cells
2500
10000
40000
16562
Δ
x (m)
4.0
2.0
1.0
CPU time (s)
0.375
1.247
6.621
2.510Figure 12
Circular dam break: water surface profiles resulting from simulations on different grids.


### 5.5. Bore Wave Past a Dip

This test gives unsteady results of a right-going bore wave past a downward hump. The Froude number isFr=2, and the initial conditions to the left and the right of the bore are
(49)hR=1m,uR=vR=0ms-1,hL=hR2(1+8Fr2-1)m,uL=Fr(1-hRhL)ghRm/s,vL=0m/s.
The initial discontinuity is located at x=-5.5m. The bottom function b(x,y) is defined by
(50)b(x,y)=-0.4e0.2(12.5-x2-y2).
The computational domain is a square [-6m,10m]×[-10m,10m]. We use inflow boundary conditions at the west boundary and outflow boundary conditions for the other boundaries. Simulation is carried out for 2 s when the bore has passed the dip and travels towards the eastern boundary.Instead of using a uniform fine mesh to decompose the entire computational domain, adaptive grid on the coarse background grid 160 × 200 with one level local refinement is used.ϕsub=0.35 and ϕcoar=0.15 are prescribed for grid adaptation, as Figure 13(a) shows. The complex wave pattern is captured by the adapted grid with steep surface gradient; the dynamically adaptive grid ranges from 32000 to 44684. The results of the simulation are presented in Figures 13(b) and 13(c) in terms of 3D water surface and depth contours, which are compared very well with those presented by Liang [18] and Tang [24]. Comparing Figure 13(a) with Figure 13(c), it is obvious that the wave fronts can be partly predicted by the adapted grid.Bore past a dip: (a) adapted grid att=2 s; (b) 3D water surface; (c) depth contours.
(a)
(b)
(c)


### 5.6. Dam-Break Flooding over Three Humps in a Closed Channel

This test of a 2D dam-break wave traveling over an initially dry and rough floodplain with three hillocks has been widely accepted as a standard benchmark to assess the adequacy of a numerical model for realistic flood modeling applications. The dam break occurs in a 75 m long and 30 m wide closed container, with the bed topography defined by(51)zb(x,y)=max⁡[0,1-18(x-30)2+(y-6)2,max⁡1-18(x-30)2+(y-24)2,max⁡3-310(x-47.5)2+(y-15)2].The dam is initially placed atx=16 m to hold a tranquil water body with a surface elevation of 1.875 m. The rest of the domain is dry and a global Manning coefficient is set to be n=0.018. In this test, grid enrichment and coarsening are undertaken according to ϕsub=0.07 and ϕcoar=0.05.Figure14(a) illustrates the bed topography within the computational domain and the upstream still water retained by the dam. Figures 14(b), 14(c), 14(d), 14(e), and 14(f) illustrate the 3D view of free surface elevation at t=2 s, 6 s, 12 s, 30 s, and 300 s, respectively. As illustrated in Figure 14(c), for t=6 s, while the violent dam-break flow continues on moving downstream and interacts with the three hills, a reflected shock is developed after the depression wave hit the western boundary wall and propagating downstream. At t=12 s, the flood water is passing either side of the large hill and starts to flood the lee of the hill. The upstream directed reflection shock waves from the small hills have coalesced into a nearly straight bore. The main flood wave has reflected against the large hill. By t=300 s, steady state is achieved, with the peaks of the small hills no longer submerged. Similar results were obtained by [6, 17]. Figure 15 shows the evolution of the number of cells in this adaptive grid, which initially increases from 2364 to 6126 and then decreases to a steady value of about 2700 in virtue of the flow becoming steadier, this essentially indicates the efficiency by using dynamically adaptive grids for flood inundation problems.Dam-break wave propagating over three humps: 3D water surface elevation at different time.
(a)
(b)
(c)
(d)
(e)
(f)Figure 15
Dam-break wave propagating over three humps: time evolution of the dynamically adaptive grid.


### 5.7. Flood Flows over Complex Topography

In order to test the applicability of the present model to real case studies characterized by an irregular topography, a practical application related to river flood over realistic topography is presented. The study area is a river in southern China. A free outgoing condition is imposed along all boundaries and the Manning coefficient is chosen to be 0.019 over the whole domain.Figure16(a) illustrates the plan view of the bed elevation together with the locations of the gauges; Figures 16(b) and 16(c) show that the peak-flood discharge is imposed at the inflow boundary and the field-measured water level is specified at the outflow boundary, respectively. As t=0 s, the velocity field is set to u=v=0, whereas the initial water level is set to 17.366 m in the entire computational domain.(a) Bed elevation with the locations of the surveyed points; (b) the predefined discharge at inflow boundary; (c) the field-measured water level at outflow boundary.
(a)
(b)
(c)Simulations are executed for 300 hours after the flood as mentioned previously. Figure17 presents the velocity field of regions A, B, and C. The close-up shows that the major stream flow is confined along the main river channel, which indicates that the simulated velocity field is quite reasonable. It is worth noting that this is a sandbar in region C where the bed elevation is much higher than the water surface elevation of its nearby area; thus an obviously detour flow phenomenon is shown in Figure 17(c), validating the wetting and drying method in the present model. Figure 18 shows the comparison of the measured and simulated water flow at three gauges, and the computed results match the measured water flow very well. Besides, the suddenly rising phenomena of water level are shown obviously at the start of the flood tide and then the water level began to drop. All velocities were monitored during the simulation, and we noticed that no unphysical velocities (bigger than 25 m/s) happened.The velocity field of regions A, B, and C.
(a)
(b)
(c)Comparisons of the field-measured and the simulated water levels at gauges A, B, and C.
(a)
(b)
(c)


## 5.1. Still Water Test with Wet/Dry Front

This test is chosen to verify that the present model indeed preserves the stationary steady flow with wet/dry fronts over an uneven bottom. The numerical experiment takes place in a square pool1m×1m. At the centre of the pool, a symmetric bump at the bottom is mathematically defined as the form
(28)z(x,y)=max⁡[0,0.25-5*((x-0.5)2+(y-0.5)2)].The still water as the initial condition covering partially the bump in a frictionless domain is fully closed by vertical walls; the water surface elevation is 0.1 m and 0.3 m. A simulation oft=200 s is performed on a grid of 50 × 50 cells. Figure 2 shows the computed 3D water surface in the different conditions. Similarly, Figure 3 presents the water level in cells versus the x-coordinate of the cell centroid along the line y=0.5 m. In terms of [19, 20], the errors,
(29)L1(η)=1Nc∑i=1Nc|ηi-ηi¯|,L1(q)=1Nc∑i=1Nc|qi-qi¯|,
where q=hu2+v2, ηi¯, and qi¯ denote the analytical solutions and Nc is the number of cells with positive water depth, are listed in Table 1. The errors are sufficiently small to the level of round-off error of a double accuracy calculation. In addition, ε1=maxt((Vt-V0+ΔVt)/V0)=7.9×10-16, ε2=0.0 represents the maximum relative error in global mass conservation over the simulation duration under the two conditions, in which V0 and Vt are the initial volume of water and the computed volume of water at the time t, respectively, and ΔV represents the net volume transported out of the domain from the beginning of the simulation period to the time t.Table 1
Still water test: error norms for free surface elevation and unit discharge.
Cases
L
1
(
η
)
L
1
(
q
)
Water level = 0.1 m
5.43
×
1
0
-
16 m
6.0
×
1
0
-
16 m2/s
Water level = 0.3 m
1.22
×
1
0
-
18 m
2.2
×
1
0
-
17 m2/sStill water with wet/dry front: the 3D water surface aftert=300 s at different cases.
(a)
(b)Still water with wet/dry front: the water surface elevation along the liney=0.5 m at different cases. The results are obtained at t=300 s.
(a)
(b)The results show that the stationary solutions are well maintained, and the current numerical scheme is able to preserve the lake at rest solution involving wet/dry interface and hence is well balanced.


## 5.2. Moving Shorelines in a Parabolic Bowl with Friction

Thacker [21] provides a good benchmark test case for the present shallow-flow model in dealing with wetting and drying and bed slope as well as the friction source term. This paper firstly deduced the analytical solution of the test case.Consider the case where the motion of shallow water in a basin is governed by the equations [21]
(30)∂ζ∂t+∂(zb+ζ)∂x+∂(zb+ζ)∂y=0,(31)∂u∂t+u∂u∂x+v∂u∂y+τu+g∂ζ∂x=0,(32)∂v∂t+u∂v∂x+v∂v∂y+τv+g∂ζ∂y=0,
where ζ(x,y,t) is the height of the water surface above mean water level, zb is the bottom surface, u(x,y,t) is the averaged component of the water current to the east, v(x,y,t) is the depth-averaged velocity component of the water current to the north, τ is the bottom friction parameter which is considered to be a constant (related to the bed friction coefficient by Cf=hτu2+v2), g is the acceleration due to gravity, and t is the time.Assuming thatu was function of t only, u=u0(t), and v=0, then (31) and (32) imply that
(33)ζ(x,y,t)=ζx=0(t)+xζ1(t),(34)ζ1(t)=-1g(du0(t)dt+τu0(t)).The bed profile of the domain is defined by(35)zb(x)=h0(xa)2,
where h0 and a are constants so that flow takes place above parabolic bottom topography.Substituting (33), (34), (35), u=u0(t), and v=0 in (30) gives
(36)dζ0(t)dt+xdζ1(t)dt-2u0(t)h0xa2+u0(t)ζ1(t)=0.
Following Thacker, equate the time-varying coefficients of the linearly independent terms, 1 and x, respectively,
(37)dζ0(t)dt+u0(t)ζ1(t)=0,(38)xdζ1(t)dt-2u0(t)h0xa2=0.
Substituting (34) in (38),
(39)d2u0(t)dt2+τdu0(t)dt+2gh0u0(t)a2=0.
Substituting (34) in (37),
(40)dζ0(t)dt-1gu0(t)du0dt-τgu0(t)2=0.
The auxiliary equation for (39) is
(41)λ2+τλ+2gh0a2=0.
The roots of this equation are
(42)λ=-τ±τ2-p22,p=8gh0a2.In this test case, the computational domain has planned dimensions of-5000m×5000m; the coefficients are a=3000m, h0=10m, and τ=0.001s-1. So, considering the case with τ<p, a solution of (39) is
(43)u0(t)=Be-τt/2sinst,s=p2-τ22,
where B=5m/s is a constant, obtained by using given values for u0(t) and u0′(t). Substituting (43) in (40) and integrating with respect to t give(44)ζ0(t)=a2B2e-τt8g2h0(-sτsin2st+(τ24-s2)cos⁡2st)-B2e-τt4g.
Substituting (43) in (34) gives
(45)ζ1(t)=-e-τt/2g(Bscos⁡st+τB2sinst).
Substituting (44) and (45) in (33) gives
(46)ζ(t)=-e-τt/2g(Bscos⁡st+τB2sinst)x+a2B2e-τt8g2h0[-sτsin2st+(τ24-s2)cos⁡2st]-B2e-τt4g.
Thus the analytical water level η(x,t)=max⁡(h0+ζ,zb(x)) is given by
(47)η(x,t)=max⁡([-sτsin2st+(τ24-s2)cos⁡2st]h0-e-τt/2g(Bscos⁡st+τB2sinst)xmmax+a2B2e-τt8g2h0[-sτsin2st+(τ24-s2)cos⁡2st]mmax-B2e-τt4g,zb(x)[-sτsin2st+(τ24-s2)cos⁡2st]).We compared the numerical solution of the model in this paper and the analytical solution; the simulation lasted 6000 s. Figure4 shows a good agreement between the numerical predictions and analytical solutions of the free surface profile at different times; the moving shoreline is accurately simulated with no signs of amplified spurious oscillations and due to the friction effect, the sloshing amplitude of water level is damping with increasing simulation time. Figure 5 shows the comparison between numerical and analytical solutions for water depth at points −2750 m, 50 m, and 2750 m. The numerical predictions are found to match perfectly with those given by the analytical solutions and no obvious distortion is detected at the wet-dry interface. The results validate the wetting and drying method of the present model, as well as the source terms treatments.Sloshing motions in a vessel with parabolic bottom topography at different output times: (a)t=1000 s; (b) t=2000 s; (c) t=3000 s; (d) t=4000 s; (e) t=5000 s; (f) t=6000 s.
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
Moving shorelines in a parabolic bowl with friction: comparison between numerical and analytical results for water depth at points −2750 m, 50 m, and 2750 m.


## 5.3. 1D Steady Flows over a Steep Bump with WSDGM

The goal of this test case is to investigate the ability of the present model to converge to a nonstationary steady state. The test concerns 1D steady flows in a[-10m≤x≤10m], which is discretized with 0.1 m grid spacing to make the test more severe, frictionless channel. The bottom profile, characterized by the presence of a steep bump, is described by the following equation [22]:
(48)zb(x)={0.8(1-x24)-2m≤x≤2m0elsewhere.In order to demonstrate the capability of the WSDGM reconstruction, the results of supercritical flow, subcritical flow, and transcritical flow with a shock are displayed in Figures6, 7, and 8. The three cases were simulated using a mesh of 400 (Δx=0.05m) cells to steady states.1D steady flow over a steep bump-supercritical flow: comparison between reference solutions and numerical results.
(a)
(b)
(c)1D steady flow over a steep bump-subcritical flow: comparison between reference solutions and numerical results.
(a)
(b)
(c)1D steady flow over a steep bump-transcritical flow with a shock: (a) water surface elevation and (b) discharge.
(a)
(b)For the case of supercritical flow over the bump, we impose a constant discharge per unit widthq=1.5m2/s,h=0.25m at the upstream boundary and the downstream boundary set to be transmissive [22]; the Fr range from 2.31 to 3.83. Figure 6, in the subdomain [-3m≤x≤3m], shows the numerical results of the DGM (Frlim⁡=0, (a)), SGM (Frlim⁡→∞, (b)), and WSDGM (Frlim⁡=1, (c)) reconstruction, respectively. Figure 6(b) shows that the numerical scheme with Frlim⁡→∞ converges toward the steady solution worse than Frlim⁡=1. The results verified that pure SGM reconstruction is unsuitable when the flow is supercritical.On the other hand, if the flow is subcritical (e.g.,q=1.0m2/s at the upstream and h=1.7 at the downstream), then the Fr range from 0.14 to 0.41; Figure 7(a) shows a dip in the water elevation upstream and downstream of the bump, where the bottom profile is not smooth. Moreover, the pure DGM reconstruction does not satisfy the C-property.In Figures6(c), 7(c), and 8, excellent agreements are observed between the numerical predictions and analytical solutions under different conditions. All in all, the results prove that WSDGM performs better than the schemes based on pure DGM or SGM reconstructions. For this reason, the parameter Frlim⁡ was set at 1 also in the following tests.


## 5.4. 2D Circular Dam-Break Problem

This problem is based on the hypothetical test case studied by Alcrudo and Garcia-Navarro [23], which involves the breaking of a circular dam. In this problem, a circular discontinuity over a 200 m × 200 m is investigated. The circle’s radius is 50 m and it is centered at the middle of the square. The initial states are 10 m water depth inside the circle, and the water is at rest. The initial water depth outside the circle is considered for two cases, that is, 5 m and dry bed. The wall is then assumed to be removed completely and no slope or friction is considered. It is assumed that the circular dam is broken instantaneously and water surface profile, 2 s after dam failure, is studied.In these examples, the computational grid consists of 100 × 100 cells with each cell 4 m2 as the background grid. During the simulation, the adaptation criteria ϕsub=0.08 and ϕcoar=0.05 are used. Figure 9 shows the computed water surface profile 2 s after dam failure outside water depth equal to 5 m and dry, respectively. Obviously, the circular dam-break waves spread and propagate radically and symmetrically as the water drains from the deepest region. Figure 10 presents the final numerical outputs at 2 s after the dam failed in terms of adapted grid. From the 3D water surface and depth contours, it is conspicuous that the sharp shock front is well captured and the depth contours show that most of the depth contour lines appear to be circular. The case outside the water depth is equal to 5 m, and the flow regime changes from subcritical to supercritical. This case can be a good test for validation of a numerical scheme. In the case of dry-bed simulation, it can be seen that no bore forms, but instead a rarefaction wave extends into the dry region (Figure 11). The scheme is capable of handling the dry-bed problem. Similar results were obtained by Alcrudo and Garcia-Navarro [23] and Zoppou and Roberts [5].Circular dam break: (a) 3D water surface; (b) depth contours for the circular dam break att=2 s with the downstream water depth set at 5 m.
(a)
(b)Figure 10
2D circular dam break: adapted grid att=2 s.Circular dam break: (a) 3D water surface; (b) depth contours att=2 s with the downstream water depth set dry.
(a)
(b)In order to perform the advantage of the adaptive grid-based model, the simulations are also run on uniform grids with different resolutions. The results are shown in Figure12 in terms of water surface profile. The structured adaptive grid-based prediction matches closely the one produced on the finest uniform grid with 200 × 200 cells, which indicates similar accuracy of the two simulations. On the adaptive grid system, the mesh is refined to those areas with sharp surface gradients. Therefore, the computational efficiency is greatly improved due to the fact that the total number of cells is reduced. As listed in Table 2, the new adaptive grid can save more than twice of the computational cost, compared with the simulation on a uniform grid with 200 × 200 cells.Table 2
Circular dam break: CPU time required by different simulations.
Grid
Uniform
Adaptive
Number of cells
2500
10000
40000
16562
Δ
x (m)
4.0
2.0
1.0
CPU time (s)
0.375
1.247
6.621
2.510Figure 12
Circular dam break: water surface profiles resulting from simulations on different grids.


## 5.5. Bore Wave Past a Dip

This test gives unsteady results of a right-going bore wave past a downward hump. The Froude number isFr=2, and the initial conditions to the left and the right of the bore are
(49)hR=1m,uR=vR=0ms-1,hL=hR2(1+8Fr2-1)m,uL=Fr(1-hRhL)ghRm/s,vL=0m/s.
The initial discontinuity is located at x=-5.5m. The bottom function b(x,y) is defined by
(50)b(x,y)=-0.4e0.2(12.5-x2-y2).
The computational domain is a square [-6m,10m]×[-10m,10m]. We use inflow boundary conditions at the west boundary and outflow boundary conditions for the other boundaries. Simulation is carried out for 2 s when the bore has passed the dip and travels towards the eastern boundary.Instead of using a uniform fine mesh to decompose the entire computational domain, adaptive grid on the coarse background grid 160 × 200 with one level local refinement is used.ϕsub=0.35 and ϕcoar=0.15 are prescribed for grid adaptation, as Figure 13(a) shows. The complex wave pattern is captured by the adapted grid with steep surface gradient; the dynamically adaptive grid ranges from 32000 to 44684. The results of the simulation are presented in Figures 13(b) and 13(c) in terms of 3D water surface and depth contours, which are compared very well with those presented by Liang [18] and Tang [24]. Comparing Figure 13(a) with Figure 13(c), it is obvious that the wave fronts can be partly predicted by the adapted grid.Bore past a dip: (a) adapted grid att=2 s; (b) 3D water surface; (c) depth contours.
(a)
(b)
(c)


## 5.6. Dam-Break Flooding over Three Humps in a Closed Channel

This test of a 2D dam-break wave traveling over an initially dry and rough floodplain with three hillocks has been widely accepted as a standard benchmark to assess the adequacy of a numerical model for realistic flood modeling applications. The dam break occurs in a 75 m long and 30 m wide closed container, with the bed topography defined by(51)zb(x,y)=max⁡[0,1-18(x-30)2+(y-6)2,max⁡1-18(x-30)2+(y-24)2,max⁡3-310(x-47.5)2+(y-15)2].The dam is initially placed atx=16 m to hold a tranquil water body with a surface elevation of 1.875 m. The rest of the domain is dry and a global Manning coefficient is set to be n=0.018. In this test, grid enrichment and coarsening are undertaken according to ϕsub=0.07 and ϕcoar=0.05.Figure14(a) illustrates the bed topography within the computational domain and the upstream still water retained by the dam. Figures 14(b), 14(c), 14(d), 14(e), and 14(f) illustrate the 3D view of free surface elevation at t=2 s, 6 s, 12 s, 30 s, and 300 s, respectively. As illustrated in Figure 14(c), for t=6 s, while the violent dam-break flow continues on moving downstream and interacts with the three hills, a reflected shock is developed after the depression wave hit the western boundary wall and propagating downstream. At t=12 s, the flood water is passing either side of the large hill and starts to flood the lee of the hill. The upstream directed reflection shock waves from the small hills have coalesced into a nearly straight bore. The main flood wave has reflected against the large hill. By t=300 s, steady state is achieved, with the peaks of the small hills no longer submerged. Similar results were obtained by [6, 17]. Figure 15 shows the evolution of the number of cells in this adaptive grid, which initially increases from 2364 to 6126 and then decreases to a steady value of about 2700 in virtue of the flow becoming steadier, this essentially indicates the efficiency by using dynamically adaptive grids for flood inundation problems.Dam-break wave propagating over three humps: 3D water surface elevation at different time.
(a)
(b)
(c)
(d)
(e)
(f)Figure 15
Dam-break wave propagating over three humps: time evolution of the dynamically adaptive grid.


## 5.7. Flood Flows over Complex Topography

In order to test the applicability of the present model to real case studies characterized by an irregular topography, a practical application related to river flood over realistic topography is presented. The study area is a river in southern China. A free outgoing condition is imposed along all boundaries and the Manning coefficient is chosen to be 0.019 over the whole domain.Figure16(a) illustrates the plan view of the bed elevation together with the locations of the gauges; Figures 16(b) and 16(c) show that the peak-flood discharge is imposed at the inflow boundary and the field-measured water level is specified at the outflow boundary, respectively. As t=0 s, the velocity field is set to u=v=0, whereas the initial water level is set to 17.366 m in the entire computational domain.(a) Bed elevation with the locations of the surveyed points; (b) the predefined discharge at inflow boundary; (c) the field-measured water level at outflow boundary.
(a)
(b)
(c)Simulations are executed for 300 hours after the flood as mentioned previously. Figure17 presents the velocity field of regions A, B, and C. The close-up shows that the major stream flow is confined along the main river channel, which indicates that the simulated velocity field is quite reasonable. It is worth noting that this is a sandbar in region C where the bed elevation is much higher than the water surface elevation of its nearby area; thus an obviously detour flow phenomenon is shown in Figure 17(c), validating the wetting and drying method in the present model. Figure 18 shows the comparison of the measured and simulated water flow at three gauges, and the computed results match the measured water flow very well. Besides, the suddenly rising phenomena of water level are shown obviously at the start of the flood tide and then the water level began to drop. All velocities were monitored during the simulation, and we noticed that no unphysical velocities (bigger than 25 m/s) happened.The velocity field of regions A, B, and C.
(a)
(b)
(c)Comparisons of the field-measured and the simulated water levels at gauges A, B, and C.
(a)
(b)
(c)


## 6. Conclusions

A Godunov-type numerical model based on the finite volume method has been developed and tested to simulate shallow-water flows over a wet or dry bed with complex domain topography. The physical domain has been discretized on uniform and nonuniform structured grids with dynamic refinement. The numerical model solves the two-dimensional shallow-water equations derived as conservation laws with the flux and source terms conveniently balanced so that there is no need for further numerical treatment. The weighted surface-depth gradient method (WSDGM) computes water depth at the cell boundaries through a weighted average, based on the local Froude number of the extrapolated values deriving from DGM and SGM reconstructions. The method is simple, accurate, robust, and easy to implement and can be used to solve both steady and unsteady shallow-water problems. The present results demonstrate the accuracy of the new finite volume method and its capability to simulate water flows in the hydrodynamic regimes considered.



---

*Source: 209562-2014-07-10.xml*