A Well-balanced 2-D Model for Dam-break Flow with Wetting and Drying

This paper presents a two-dimension (2-D) numerical model for simulating dam-break flow involving wet-dry fronts over irregular topography. The Central-upwind scheme is chosen to calculate interface fluxes for each cell edge. A second order spatial linear reconstruction with multidimensional limiter and second order TVD RungeKutta scheme are chosen to acquire high order accuracy in space and time. Non-negative water reconstruction of variables at cell interfaces and compatible discretization of slope source term lead to stable and well-balanced scheme for hydraulics over irregular topography. The friction term is discretized with a semi-implicit scheme for numerical stability when very small water depth exists. An accurate and effective technique is presented for tracking wetting-drying interfaces during the process of wave front propagation on dry bed. The capacity and accuracy of current model are verified by several benchmark tests as well as a real dam-break case, and good performances are achieved in tests.


Nomenclature t
The simulating time h The water depth x The horizontal coordinate y The horizontal coordinate u The depth-averaged velocity components in x direction v The depth-averaged velocity components in y direction Z The bed elevation g The gravitational acceleration ν t The eddy viscosity coefficient n b The Manning's roughness coefficient S f x The friction slope terms in the x direction S f y The friction slope terms in the y direction η The water surface level The components of unit normal vector in the  direction The components of unit normal vector in the  direction l ik The length of the kth edge of control volume i F ik The numerical convective flux across l ik F ̃The diffusive flux across l ik in the  direction G ̃The diffusive flux across l ik in the  direction U R The one-sided local speed of propagation on right side of kth edge U L The one-sided local speed of propagation on left side of kth edge S R The one-sided local speeds of propagation on right side of kth edge S L The one-sided local speeds of propagation on left side of kth edge The threshold for defining wet and dry cells

Introduction
Dam-break flows over irregular bed often experience the cases as transcritical flows, steep bed slope, very small water depth, cells' wetting and drying, wave propagation.Inappropriate hydraulic simulation may leads to inaccurate and unstable prediction.Hence, it's necessary to apply a robust, accurate and effective hydraulic model to predict Hydraulic parameters of dam-break flow.
By using Godunov-type scheme, complicated shallow water flow phenomena such as transcritical flows, shock-type flows and moving wet-dry interface of a water wave front can be appropriately simulated.To solve a Godunov-type scheme, approximate Riemann solver is normally adopted to estimate the numerical fluxes of hyperbolic system and various numerical schemes have been proposed to solve Riemann problems.As explained in the next section, the Central-upwind method which is applied in current study, does not require computationally expensive decomposition of numerical flux on the basis of eigenvalues and furthermore, only an estimation of largest and smallest eigenvalues of the Jacobian matrix, which leads to a significant reduction of computational cost.
Flow over initially dry bed is very common in dam break flows, which involves complicated boundary conditions.For irregular topography, both positive and negative bed slopes generally exist, which may leads to cell drying and wetting with moving fronts which could not be easily solved by horizontal boundary condition.Some techniques have been developed in using finite volume method and shallow water equations.Zhao et al. [19] and Sleigh et al. [14] introduced two similar schemes to track the wetting and drying fronts, in which cells are divided into wet, dry and partially dry types according to two tolerances.Brufau et al. [2,3] proposed a technique using unsteady wetting and drying conditions in flows and claimed that the method gave zero mass error, which is valid for an FVM with only first-order accuracy.Falconer et al. [5] and Falconer et al. [6] developed a wetting and drying method for regular grid finite difference model, which is recently refined for triangular grids by Xia et al. [17].In present study, a technology for tracking wet-dry front is developed combing with method of Brufau et al. [3] to achieve zero mass error.
It is well-known that accuracy is the most important aspect for flow solver since it has a direct influence on the number of computational cells required.This means that a higher-order implementation involving a piecewise linear reconstruction is necessary.Higher order schemes often produce nonphysical oscillations which can be effectively suppressed by using limiters.In the present study, the multidimensional limiter proposed by Jawahar et al. [9] is adopted to calculate the limited gradient for reconstruction of variables.Moreover, a second order accuracy in time could be obtained in current model by apply Runge-Kutta scheme.
Based on the efficient divergence form of the slope source term proposed by Valiani et al. [16], Hou et al. [8]developed a novel slope source term treatment which is devised to transform the slope source of a cell into a flux form.By splitting the integral of the bed slope source term over a cell into those of the sub-cells, higher accuracy can be achieved by the novel treatment than that proposed by Valiani et al. [16].This method can strictly preserve the wellbalanced property and can be conveniently employed with second order or even higher order schemes.In addition, this treatment is able to handle the occurrence of wet-dry fronts, in conjunction with the non-negative water depth reconstruction.

Numerical Scheme
The 2-D shallow water is adopted in the current study as governing equations which can be described as: The 2-D shallow water equations constitute a hyperbolic system which can be presented in the following vector form: ̃= ( in which t = time, h = the water depth, x and y are horizontal coordinates, u and v are the depth-averaged velocity components in x and y directions respectively, Z is the bed elevation, g is the gravitational acceleration, ν t is eddy viscosity coefficient given by x and S f y are friction slope terms in the x, y directions, respectively, which can be determined by conventional formulas involving Manning roughness coefficient Besides, in this work, water surface level η is used in the second order spatial reconstruction and the non-negative water depth reconstruction.η can be calculated as η = h + Z.

Discretization of Flow and Sediment Governing Equations
The coupled system is discretized on an unstructured triangular grid by finite volume method.
in which   and   are components of unit normal vector in the  and  directions, l ik is the length of the kth edge of control volume i, F ik is numerical convective flux across l ik , (F ̃nx + G ̃ny ) is the diffusive flux across l ik .

Central-upwind scheme for the interface flux
Central-upwind schemes on general triangular grids for solving two-dimensional systems of conservation laws are developed by Kurganov et al. [11], which enjoy the main advantages of the Godunov-type central schemes, i.e. simplicity, universality and robustness and can be applied to problems with complicated geometries.The triangular centralupwind schemes are based on the use of the directional local speeds of propagation and are a generalization of the central-upwind schemes on rectangular grids, introduced in [10].
Applying the Central-upwind scheme, the convective fluxes in Eqs. ( 7) could be estimated by: in which F ⃗⃗ (U L ) ik and F ⃗⃗ (U R ) ik are normal fluxes on the left and right sides of the kth edge, respectively.S R and S L are one-sided local speeds of propagation on right and left sides of kth edge, respectively, and can be determined by with ik being the N eigenvalues of the Jacobian matrix ∂F n (U)/ ∂U using reconstructed variables (U L ) ik or (U R ) ik .
It can be deduced that, when a cell is dry or nearly dry, S R + S L ≈ 0 exists in denominator of equation Eqs.(8) which may cause numerical instability.In order to handle this situation that both S R and S L are zero (or very closeto zero), following the suggestion in (Bryson et al., 2011), the scheme Eqs.(8) reduces to: +   < 10 −8 (10)

Spatial linear Reconstruction
In this paper, the following 2-D linear reconstruction is employed:  ˚(, ) =   + (  ̆) ( −   ) + (  ̆) ( −   ), (, ) ∈   (11) in which U ˚(x, y) is the reconstructed value of variables at point (x, y) inside of cell i, (U x ̆)i and (U y ̆)i are component-wise approximation of numerical derivatives which are computed via a nonlinear limiter, used to minimize the oscillations of the reconstructions.
In current study, the multidimensional limiter proposed by Jawahar et al. [9] is adopted to calculated limited gradient ∇U i ̆ within a cell by taking the weighted average of three representative unlimited gradients.
In order to preserve the well-balanced property for second order schemes, as suggested by Audusse et al. [1], surface levels η ˚L M , η ˚R M , water depths h ˚L M , h ˚R M , flow discharges (hu) ˚L M , (hu) ˚R M ,(hv) ˚L M and (hv) ˚R M are reconstructed at the midpoint M of considered edges.The reconstructed bed levels at M are given by  ˚  =  ˚  − ℎ ˚  ,  ˚  =  ˚  − ℎ ˚  (12) Besides, flow velocities required at M are computed as  ˚  = (ℎ) ˚  /ℎ ˚  ,  ˚  = (ℎ) ˚  / ℎ ˚  ,  ˚  = (ℎ) ˚  /ℎ ˚  ,  ˚  = (ℎ) ˚  /ℎ ˚ A threshold  will be introduced for defining wet and dry cells.The second order reconstructions are only applicable to the wet cells.

Non-negative water depth reconstruction
A robust and efficient reconstruction approach to preserve non-negative water depth which is suggested by Liang et al. [12] is adopted in present study.The Z M bed elevation at the midpoint of the considered edge is calculated as: Then the non-negative water depth values on both sides are reconstructed as: The discharges at M are in turn reconstructed as:

Treatment of the source term 2.5.1. Well-balanced treatment of the slope source term
Hou et al. [8] developed a novel slope source term treatment which is devised to transform the slope source of a cell into a flux form, which can strictly preserve the well-balanced property and can be conveniently employed with second order or even higher order schemes.In addition, this treatment is able to handle the occurrence of wet-dry fronts, in conjunction with the non-negative water depth reconstruction.The vector of source term F nb (U) ik at the considered faces k in cell i becomes where   and ℎ   are obtained by spatial linear reconstruction.

Calculation of derivatives in source terms
The approach adopted by Mohammadian et al., [13] is used to calculate the unlimited gradient of variables.A similar approach may be also applied to calculate the diffusive terms

Friction source term treatment
In current study, a simple semi-implicit treatment suggested by Yoon et al. [18] is adopted to deal with very shallow water depth.
in which the superscript n denotes the time step.

Time integration and CFL condition
In order to obtain the second-order accuracy in time and retain the stability of the model, the twostage explicit TVD Runge-Kutta method is applied in current study.The time step is limited by the Courant-Friedrichs-Lewy (CFL) condition.

Treatment of wetting and drying fronts
Wetting and drying fronts need to be specially treated to retain the numerical stability when very tiny water depth is introduced.In current model, a scheme of wetting and drying treatment is proposed which can be summarized as: 1.A tolerance water depth  is introduced to classify the wet and dry cell.In present model,  = 10 −6  is used.In the dry cell, first order reconstruction will be applied to maintain the numerical stability.
2. For dry cell with ℎ   < , only continuity equation will be solved to deal with potentially wetting.
3. To deal with the cell-drying, if a cell with ℎ  +1 > , it's treated as completely dry cell in next time-step and all hydraulic parameters are set to zero in this cell.

Numerical Tests 3.1. Quiescent water around a hump with sediment deposition
The first test is applied to verify the wellbalanced property preserving of present model involving wet-dry interfaces.To implement this test, the elevation of a hump on flat bed is defined as: The hump is located at the center of a 8m×8m computational domain, the height of the bump is 2m.A quiescent lake around the dump is defined with initial water surface elevation of 1m, so the test could involve the wet-dry interface.Fig. 1 shows the 3d views of simulated bed profile and still water surface at t=50s.The undisturbed water surfaces are observed through the whole simulating process.

2-D shorelines tracking in a parabolic bowl
This test is adopted here to investigate the hydraulic model and the accuracy of tracking the wetting and drying fronts.The bottom topography with the center (x 0 ,y 0 ) is defined as: in which h 0 is the water depth at the center of the domain, a is the distance from the center to the edge of the shoreline.Since the topography is set to be frictionless in this test, as described in [15], the periodic analytical solution of the evolutions of surface elevation, water depth and velocities can be computed using following equations:  Fig. 3 shows the comparison between simulated water surface profiles and analytical solutions at t=3T and 3.5T, respectively, where T represents the circulation period of the pool.It can be observed that the calculated free surface from current model agrees well with the analytical solution, no obvious distortion is observed near the shorelines, the treatment of wetting and drying boundaries successfully handle the task of tracking moving wet-dry fronts.Fig. 4 shows the velocity of circulating pool at t=3.5T, no unphysical high velocity is observed near the shoreline where very shallow water exists.The around area preserve zero velocity when it becomes dry.

2-D partially dam-break flow on initially dry bed
For examining the numerical performance of the present scheme, a 2-d dam break problem with rapid varying unsteady flow is chosen as test cases.This test case is firstly introduced by Fennema et al. [7] in a numerical method study, which has been widely used by many researchers.The original case is mostly applied on an initially wet and fixed bed.It is adopted here for the aim of testing the capacity of current model to simulate the front wave propagation over an initial dry bed with wet-dry interface tracking, with a particular attention to the 2D aspects of the flow motion; and the model performance in prediction of fast erosion under high energetic flow in large scale.
The modeling domain area is a 200m×200m basin over a flat and dry bed.A 10m-thick dam splits the basin into two equal-sized regions.The water depths are 10m and 0m on the left and right sides of the dam, respectively.At t = 0s, a 75m wide breach centered at y = 125m is assumed to form instantaneously.The duration of the simulation is 12s.The initial velocity of the whole modeling area is 0m/s and the outlet boundary at x = 200m is specified with a free out flow boundary condition, meanwhile all other boundaries are set to be standard wall condition.The simulating domain is discretized into 6806 triangular cells, as shown in Fig. 5.  6 shows 3D views of the water front-wave propagations and predicted bed erosion at varying times (first row).It can be seen that a shock wave forms and propagates downstream and a depression wave spreads upstream during the simulating process, no unphysical high velocities and oscillations are generated at the front-wave locations where wetting and drying interfaces are existing.The velocity field with contour of water depth is present in second row.

Conclusion
A robust two-dimensional numerical model for dam-break flow with wetting and drying is proposed in current study, based on finite volume method using unstructured triangular grids.The Central-upwind scheme can accurately estimate the convective fluxes.Current model can strictly preserve the well-balanced property with vector-form discretization of slope source term in conjunction with nonnegative water depth reconstruction.The semi-implicit method on friction slope term maintains the stability of present model.The adopted wetting and drying scheme could efficiently track the wet-dry fronts.The testing results confirm the capacity and accuracy of current model in dealing with various cases of dam-break flow over irregular bed in practical conditions.

Figure 1 .
Figure 1.Quiescent lake around a hump with initial sediment concentration and wet-dry boundaries at t=50s

Figure 2 .
Figure 2. Simulated grids and bed profile of a parabolic bowl

Figure 3 .Figure 4 .
Figure 3. Simulated shoreline profiles and analytical solution at different times