Development of a Cell-Centered Godunov-Type Finite Volume Model for Shallow Water Flow Based on Unstructured Mesh

Based on the Godunov-type cell-centered finite volume method, this paper presents a two-dimensional well-balanced shallow water model for simulating flows over arbitrary topography with wetting and drying. The central upwind scheme is used for the computation ofmass andmomentumfluxes on interface.The novel aspect of the presentmodel is a robust and accurate nonnegative water depth reconstructionmethod which is implemented in the unstructuredmesh to achieve second-order accuracy in space and to track the moving wet/dry fronts of the flow over irregular terrain. By defining the bed elevation and primary flow variables at the cell center in the nonstaggered grid system, all computational cells are either fully wet or dry to avoid the problem of being partially wetted. The developed model is capable of being well balanced and preserving the computed water depth to be nonnegative under a certain CFL restriction, which makes it robust and stable. The present model is validated against three benchmark tests and two laboratory dam-break cases. Finally, the good agreement between the numerical results by the established model and measured data of the Malpasset dam break event on a 1/400 scale physical model demonstrates the capability of the model for the real-life applications.


Introduction
The two-dimensional (2D) shallow water equations (SWEs) have been widely used to mathematically describe free surface flows over complex topography, such as river and overland flows, dam-break floods, and estuarine and coastal circulation.Because the exact analytical solutions of SWEs are only available for some simple and specific cases, it is important to develop numerical methods of SWEs with good properties for general applications in hydraulic and coastal engineering.Various numerical methods have been developed to obtain the satisfactory solutions of SWEs, such as the finite difference [1,2], the finite volume [3][4][5][6], and lattice Boltzmann methods [7].The finite volume method, which can provide solutions with good mass and momentum conservations in both structured and unstructured computational meshes, is probably the most popular numerical technique to solve SWEs.
In recent decades, Godunov-type cell-centered finite volume schemes have been applied to solve SWEs numerically due to its robustness [8][9][10][11][12][13][14][15][16].By considering the mesh element as the control volume, Godunov-type schemes can reasonably simulate the most complicated flow phenomena including flows with mix regimes and shock-wave discontinuities [9].However, when modeling the flows over arbitrary topography in realistic engineering applications, there are still challenges of handling the wetting and/or drying moving boundaries and preserving the well-balanced property [17].
Several mathematical and numerical treatments have been proposed to preserve the "lake at rest" steady states, which are known as well-balanced property [17] or Cproperty [18].For example, Bermúdez and Vazquez [18] proposed an upwind discretization method to balance the numerical fluxes and the bed slope source terms.In comparison with earlier models that used the the centred discretization method, this technique significantly improved the accuracy of numerical solution.Zhou et al. [15] adopted the surface gradient method (SGM) to reconstruct the water surface level instead of water depth to preserve the wellbalanced conditions, but special corrections were required for the vertical bed step case [19].Liang and Borthwick [20] derived a new formulation of SWEs to mathematically balance the flux gradient and source terms.However, local modifications of Riemann states are needed near the wet/dry front to prevent the generation of spurious flow in dry areas.Another balancing method based on hydrostatic reconstruction was presented by Audusse et al. [21], which had been found to be robust and effective to simulate coastal hydrodynamics [22].
To track the wetting and drying moving boundary accurately without numerical instability has attracted significant attention in solving of SWEs in recent years.This is because the velocity is calculated by dividing the discharge per unit width by water depth in finite-volume schemes; nonphysically high velocity might be produced near the wet/dry fronts, which yields negative water depths to crash the numerical simulation.Hence, a successful shallow water model should preserve the well-balanced property and the positivity of water depth at the same time.A variety of positivity preserving schemes for SWEs have been proposed and widely applied in numerical models with structure or unstructured meshes [21,[23][24][25][26][27].
A two-dimensional model, which was developed by Bryson et al. [26] to solve SWEs on triangular grids using a second-order central upwind scheme, has been proved to have the properties of being well balanced and preserving the positivity of water depth.However, this model has a main drawback that is not being able to preserve the stationary steady states when dry areas are contained in the computational domain, besides ignoring the effects of the friction terms.The main reason is that the bed elevation was defined at the mesh's corner and the flow variables defined at the mesh's center, which might produce partially wetted cells near the wet/dry fronts as shown in Figure 1(a).In order to handle partially wetted cells, the model simply redefined the slope of water surface as shown in Figure 1(b) to convert the partially wetted cells into the fully wetted ones.Obviously, this treatment violates the well-balanced property of numerical model.
To conquer the difficulties above, this study presents an improved cell-centered finite volume model based on the Godunov-type central upwind scheme to simulate shallow water flow over arbitrary topography in this paper.The model is developed using the unstructured mesh to make it more applicable to the real-life engineering problems with complicated boundaries and bathymetries.The central upwind scheme is adopted to calculate the fluxes of mass and momentum because it is Riemann-problem-solver-free and provides a robust and accurate solution.The nonnegative water depth reconstruction proposed by Liang [10] for Cartesian grids is implemented in our unstructured model to achieve second-order accuracy in space and track rapidly moving wet/dry fronts efficiently and robustly.
The paper is organized as follows.Section 2 presents the governing equations.Section 3 gives the numerical scheme including discretization, the second-order spatial reconstruction, and the treatments of the bed slope and friction terms.Section 4 tests the proposed model against three benchmark cases and two laboratory cases of dambreak flow and then applies the model to a real-life dam failure case.The conclusions are given in Section 5.

Governing Equations
Based on the hydrostatic pressure assumption, the depthaveraged 2D SWEs can be obtained by integrating 3D Navier-Stokes equations over water depth.Neglecting the kinetic and turbulent viscous terms, wind effects, and the Coriolis term, the conservative form of the depth-averaged 2D shallow water equations can be written as [27]  in which where q represents the vector of conserved variables; f and g are the flux vectors associated with the conserved variables in the and -directions, respectively; S represents the vector of source terms;  indicates time;  and  are Cartesian coordinates; ℎ is water depth;   represents bed elevation over the datum;  is the water surface, which can be expressed as  = ℎ +   (see Figure 2). and V are depth-averaged velocity components in the and -directions, respectively;  is the gravitational acceleration;   and   are the friction source terms in the and -directions, respectively.The friction source terms in the equations are estimated using Manning's formula as where   is Manning's coefficient reflecting the roughness of the bed.

Finite Volume Discretization for SWEs over Triangular
Grids.The computational domain is discretized using a number of triangular cells defined as control volumes, as shown in Figure 3.The conserved variables of system are stored at the center of each cell.Integrating (1) over a triangular control volume and applying Green's theorem yield where Ω is the control volume; Γ is the boundary of the control volume; n is the unit outward vector normal to the boundary with the components of (  ,   )  ; F(q) ⋅ n is the flux vector normal to the boundary, expressed as Rewriting the line integral terms in (4) into an algebraic form and using Euler scheme for the time derivative, (4) can be discretized as where the superscript  represents time level, Δ  is the area of cell , and  and  are the index and length of the edge of cell , respectively.

Nonnegative Water Depth
Reconstruction.The limited central difference (LCD) scheme [28], one of the most widely used limiters, is employed to calculate the variables at the center of the edge to achieve the second-order accuracy in space.Let   be the midpoint of the th edge of the cell ,  = 1, 2, 3 (see Figure 3).Then interface value  that is the component of conserved variable q at   can be calculated as where   is the value at the centroid of the cell ; r is the distance vector of midpoint   related to the cell center; ∇  is the limited gradient of the variable  over cell , expressed as where   is the limited function; ∇  is the unlimited gradient, which can be evaluated by the values at the centroid of three neighboring cells.The unlimited gradient ∇  is given by [13] ∇ where 3);   is the centroid value of the cell adjacent to the th edge of the cell .(  ,   ) are the coordinates of the centroid of the cell adjacent to the th edge of the cell .
In the framework of Godunov-type scheme, the determination of the interface fluxes requires the values at both sides of the edge.Thus, the left states of the variables at the midpoint   can be obtained using (7): Similarly, the right states of the variables at the midpoint   can be calculated by In order to obtain the final states, a single value of bed elevation suggested by Audusse et al. [21] should be defined as Then, the reconstructed water depth at the midpoint   should be redefined as Obviously, ( 14) could ensure the reconstructed water depth to be nonnegative.However, the reconstructed water depths may be very small or even zero and lead to local extremely high values of velocity, which might produce negative water depths and affect the stability of the numerical model.In order to prevent this problem, the regularization technique suggested by Kurganov and Petrova [24] is used to calculate the corresponding velocity components at the midpoint   : Water surface where  is an empirically predefined positive number with  being 10 −6 m in all simulations in this paper.It should also be mentioned that  = 10 −6 m is used to distinguish the dry and wet cells.When the water depth at the cell center is less than 10 −6 m, the cell is defined as a dry cell and the velocity is set to be zero without solving momentum equations.Hence, other flow variables at the midpoint   should be recalculated by It is trivial to prove that the reconstruction processes above do not affect the well-balanced property of the numerical scheme when the wet-bed applications are simulated.However, for dry-bed applications, a numerical technique introduced by Liang [10] is naturally borrowed here to preserve the wellbalanced solutions.As shown in Figure 4(a), a wet cell  shares a common edge with a dry cell 1, and the bed elevation of dry cell is higher than the water level at the centroid of cell .Considering a general case of "lake at rest" with  ≡ 0, V ≡ 0, and  ≡ const at wet area, after aforementioned reconstructed method, the left and right states at , which is the midpoint of the common edge, are given by Therefore, the fluxes at the interface between cell  and 1 are calculated using the bed elevation (  )  rather than the actual water surface elevation  in cell .However, fluxes at the other two cell interfaces are evaluated with the actual water surface elevation  in cell , which drives the flow into motion in the cell  and violates the well-balanced property of the scheme.The local bed modification is proposed by Liang [10] to tackle this problem.As shown in Figure 4(b), the difference between the actual and fake water level at midpoint  is calculated by Then the bed elevation and reconstructed water surface elevations at midpoint  are locally modified by subtracting Δ from original values Obviously, (  )  = ()   = ()   =  after the local bed modification above and the well-balanced property of the scheme is regained for dry-bed applications.

Interface Fluxes Calculation.
The Godunov-type centralupwind scheme is adopted in this study to calculate fluxes at interface since this scheme is Riemann-problem-solverfree with good robustness and accuracy (Kurganov and Levy [25]).The interface flux can be calculated by where q    and q    are the reconstructed variables on the left and right sides of the midpoint   , respectively. in  and  out  are the local one-sided speeds of propagation, given by where     and     are the normal velocities at the midpoint   , calculated by 3.4.Treatment of Source Terms.The source terms include slope source terms and friction terms.The slope source terms need to be discretized reasonably to exactly balance the numerical fluxes to guarantee the well-balanced property of the scheme.A well-balanced slope source term discretization proposed by Bryson et al. [26] is adopted in this paper, which had been proved to be an effective way to preserve the Cproperty.Hence, the slope source terms are approximated as where (  )  and (  )  are the components of the limited gradient of the water level  over cell  in the and -directions, respectively.
To increase the stability of the numerical model, the friction terms are evaluated by a semi-implicit method [11] in this work, which are given by 3.5.Stability Criteria.It is crucial to choose an appropriate time step for an explicit scheme to maintain its stability.In this study, the CFL condition is used to estimate the time step on triangular grids (Bryson et al. [26]) as follows: where Δ is the time step;   is the Courant number specified in the range 0 <   < 1/6;   = 0.15 is adopted in this paper.  is the corresponding altitude of the th edge of the th cell.It had been proved that the model could preserve the positivity of the water depth when the time step satisfies (25) (see [26] for mathematical proof).

Numerical Results and Discussions
In this section, the developed model is applied to a wide variety of problems including three benchmark tests, two experimental cases, and a field-scale application to test its performance.All the computational triangular grids are generated using an open grid generation package "Triangle, " developed by Shewchuk [29].

Stationary Flow with Wet/Dry Interface over Uneven
Bed.To verify the capability of the model to preserve the well-balanced property when the dry area is included in the computational domain, a stationary flow with wet/dry front interface over uneven bed is simulated.A frictionless rectangular computational domain [0, 2]×[0, 1], with the bed elevation given by [30] The initial condition is static state with a constant water level of 0.2 m.So the topography is partially submerged.The total simulation is 100 s, using a mesh of 15846 triangular grids (see Figure 5).The time step is controlled by the CFL criteria.Quiescent flow is found in the entire simulation for our model and the simulated water surface and velocity over uneven topography at 100s is plotted in Figure 5(a).We also presented the simulated results using the model developed by Bryson et al. [26] in Figure 5 in Figure 5(b), while the present model gives a better simulation in Figure 5(a).The  2 errors for ℎ and ℎ are also calculated using the expression as where  is the total number of the wet computational cells and ℎ and ℎ denote the exact solutions.The  2 errors for ℎ and ℎ given by the present model are 9.46482 − 018 m and 1.79539 − 016 m 2 /s, respectively.The results show that the stationary state is well maintained and the present model could preserve the well-balanced property up to round-off error.

Accuracy Validation.
The classical subcritical flow over bump (Goutal and Maurel [31]) is used to investigate the scheme accuracy and convergence.The flow occurs in a 25 m × 5 m frictionless channel, and the topography is given by A unit discharge   = 4.42 m 2 /s and   = 0.0 m 2 /s is applied at the upstream boundary and a fixed water depth ℎ =  = 2.0 m is used at the downstream boundary.The initial condition for water level and unit discharge at the -direction are set to be  = 2.0 m and   = 4.42 m 2 /s at each cell.The model runs 200 s to reach the steady state.Four different triangular grids with mesh sizes of Δ = 1.0 m, 0.5 m, 0.25 m, and 0.125 m are tested to evaluate the convergence of grids in this study (see Figure 6).The  2 error for  and ℎ versus the grid size Δ in logarithm scale is plotted in Figure 7.It can be found that the  2 error for  and ℎ vanishes as (Δ 2 ) when the grid is progressively refined.

Two-Dimensional Shallow Water Sloshing in a Frictionless
Parabolic Basin.Periodic flow in parabolic basins proposed by Thacker [32] is investigaed here to verify the capability of the present model to accurately simulate flows with wet/dry fronts over uneven topographies.This case has been recognized as a very difficut problem for numerical models and has been used by many researchers [13,33] to test their models.The bed elevation of parabolic basin in this case is given by where ℎ 0 and  are positive constants, representing the water depth at the center of the basin and the distance from the center to the shoreline with zero water surface elevation, respectively.The analytical solution for the periodic ocillation in parabolic basins can be expressed as where  is a constant which determines the amplitude of the motion. = √2ℎ 0 / = 2/ is the rotation frequency of flow in the parabolic basin and  is the rotation period.
It should be mentioned that the aforementioned analytical solution is only valid when  ≪  2 /2 ℎ 0 .The test is performed in a square domain with a dimension of [−10000 m, 10000 m] × [−10000 m, 10000 m] without friction.The corresponding parameters ℎ 0 = 10 m,  = 8025.5m,  = /10 = 802.55,and the rotation period  = 2/ = 3600 s, which were proposed by Liang [10], are used in this study.The computational domain is discretized by 2 × 200 × 200 structured triangles.Since the flow can never reach the boundaries, the boundary conditions have no effect Analytical Numerical on the results and close boundary conditions are applied.The initial conditions are determined according to the solutions of (30), which are evaluated at  = 0 s.The simulation is carried out until time is equal to 3. Figure 8 shows the comparison between numerical and analytical solutions of contours of water depth at  =  + /2, 2, 2 + /6 and 2 + /3, respectively.The comparison shows a very good agreement and the moving wet/dry fronts are captured accurately without spurious oscillation.The comparison between numerical and analytical solutions of time history of velocity component in -direction at point (5500, 0) and (7900, 0) is illustrated in Figure 9.It is clearly shown that a good agreement is achieved for the point (5500, 0), which remains wet all the time.For the point (7900, 0), which gets wet and dry during the simulation, little deviation between numerical and analytical solutions can be found when wetting and drying occur.In general, the developed model can well simulate the unsteady flow with wet/dry fronts over uneven topography.

Dam-Break Flow in a Converging-Diverging Channel.
A series of dam-break experiments were carried out by Bellos et al. [34] to analyze the behavior of resulting flow.One of the experiments was simulated here to test the accuracy of our model.This experiment was conducted in a convergingdiverging shaped channel with a length of 20.7 m and a uniform slope of 0.6%.The dam was located at the narrowest part of the channel, where the width was 0.6 m, as shown in Figure 10.The initial water level in the upstream of the dam was set to be 0.3 m and the downstream region was dry.The dam was removed instantly at  = 0 s and four probes located along the centerline of the channel were used to record the water depths.In this case, the upstream boundary and two sidewalls are defined as solid walls and the free outlet boundary is applied at the downstream boundary.The total simulation is 60 s, using a mesh of 3963 triangular grids.Manning's coefficient is set to be 0.012 m −1/3 s, which is also used in Xia et al. [35] and Kuiry et al. [36].The simulated and measured water depths at four different points are plotted in Figure 11.It can be found that the calculated water depths agree very well with the measurements at two upstream points (P1 and P2).A slight discrepancy between predicted water depth and measurements is observed at two downstream points (P3 and P4), which is also reported by other researchers [35,36] and is acceptable for dam-break flow simulation.

Laboratory Dam Break over a Triangular
Hump.The present model is applied to reproduce the experimental dambreak flow over a triangular bottom sill, which is a part of IMPACT project [37].All the data are available in [38].
The experiment was conducted in a rectangular channel with a length of 5.6 m and a width of 0.5 m, as shown in Figure 12.The dam was located at 2.39 m away from the end of upstream and the water level of the reservoir in the upstream was 0.111 m.A symmetric triangle hump with the height of 0.065 m and a slope of 0.14 at both sides was installed in the downstream 2.06 m away from the dam.From hump to the end of downstream end closed by a solid wall, a pool was initially filled with 0.02 m water at rest.The dam was removed instantly at  = 0 s and time series of water depth were recorded at three gauges, which were located in the downstream 1.545 m (G1), 2.535 m (G2), and 3.185 m (G3) away from the dam as shown in Figure 12.
The simulation runs for 45 s on a mesh of 8914 triangular cells.A constant Manning's coefficient   = 0.011 is used throughout the entire domain, as recommended in [37].Figure 13 shows the computed water depth at three gauges in comparison with the observed data.It can be found that the wave arrival time and water depth are predicted in a good accuracy at G1.The water depth at gauges G2 and G3 is overpredicted, which is also reported by Hou et al. [39].A possible explanation is that Manning's coefficient provided by [37] is relatively low, which leads to smaller resistance and thus more water flows over the triangle hump.Comparisons between numerical and experimental water surface profile in the centerline around the hump at  = 1.8 s, 3.0 s, and 8.4 s are plotted in Figure 14.Generally, the water profiles at different time step are well reproduced by the present model.The numerical results show that the present model is capable of modeling flow over an irregular bed with a satisfied accuracy.

2D Malpasset Dam Break.
The Malpasset dam was located on the Reyran River in southern France (see Figure 15).It was a 66.5 m high double-curvature arch dam with a crest length of 223 m.The dam failed in 1959 due to extreme rainfall and 421 people died during this catastrophic event.After the dam failure, a physical model with 1/400 scale was set up by Laboratoire Nationale d'Hydraulique (LNH) of EDF in 1964 to measure the maximum water level and flood arrival time at 9 points along the Reyran River valley.The locations of the observation points are shown in Figure 15.
Because of the varying topography and availability of measurement data, this case is often simulated by many researchers [16,27,[39][40][41][42][43][44] to test their numerical models.The computational domain is discretized by 40426 triangular meshes with 20894 nodes.Figure 16 shows triangular cells near the dam.The initial water level in the upstream reservoir is assumed to be 100 m.The downstream floodplain is considered as dry bed since the initial flow rate in the downstream is really small comparing to the dam-break flood.Manning's coefficient is set to   = 0.033 s m −1/3 for the entire computational domain ( [16,27,[39][40][41][42][43]). Figure 17 shows the computational water depth at time  = 1800 s. Figure 18 gives the comparisons between the measured data and simulated maximum water depth and flood arrival time at the gauges.It can be found that the simulated results agree well with the observed data in general, although the simulated propagation speed of flood wave is a little slower than the observation, as also found in [27,39].In fact, Zhang and Wu [44] reported that Manning's coefficient is very sensitive to the flood arrival time in this case.As tested, a reduced

Conclusions
A two-dimensional shallow water flow model has been developed based on cell-centered unstructured triangular grids to simulate the flow over arbitrary topography.The governing equation is discretized by the explicit finite volume method.The intercell fluxes are calculated using a Godunovtype central upwind scheme that is the Riemann-problemsolver-free method for hyperbolic conservation laws.In comparison with some other unstructured models based on central upwind scheme, the present model defines the bed elevation and primary flow variables at the cell center in the nonstaggered grid system to avoid the problem of partially wetted cells.The nonnegative reconstruction of the water depth method is implemented in the unstructured mesh to track the stationary or wet/dry fronts.The friction source term is discretized using a semi-implicit treatment to reduce the risk of numerical instability caused by small water depth near the wet/dry front.With the improvements above, the present model does not need to adjust the wet/dry fronts in the nonstaggered mesh and is able to ensure the static steady states in the presences of dry areas such as buildings, islands, and railroad foundation.
The developed model was tested using three benchmark cases with analytical solutions and two dam-break experiments with measured data.Good agreements were achieved between the numerical results by the established model and the data in the literature, which shows that the present model has the satisfying robustness and accuracy.The developed model was also applied to simulate a real-life dam-break case with complex topography, that is, Malpasset dam failure.The successful prediction demonstrates that the model is capable of simulating extreme flood over irregular topography with wet/dry fronts with well-balanced property and positivity preserving property.The established model can be a useful tool in flood management and engineering practice.

Figure 1 :
Figure 1: Stationary state with dry boundaries (a) real situation and (b) reconstructed situation.

Figure 2 :
Figure 2: Sketch of water surface elevation, water depth, and bed elevation.

Figure 4 :
Figure 4: Local bed modification for dry-bed applications.

Figure 5 :
Figure 5: The bed elevation, free surface, and flow velocity at  = 100 s (a) by our model (b) by the model developed by Bryson et al. [26].

Figure 6 :
Figure 6: Definition of triangular grids with mesh size Δ.

Figure 7 :
Figure 7: Convergence curves in logarithm scale for the water level and the unit discharge in -direction.

Figure 11 :Figure 12 :
Figure 11: Comparisons between the observed and calculated water depths at gauges.