A New Method for Wet-Dry Front Treatment in Outburst Flood Simulation

: When utilizing a ﬁnite volume method to predict outburst ﬂood evolution in real geometry, the processing of wet-dry front and dry cells is an important step. In this paper, we propose a new approach to process wet-dry front and dry cells, including four steps: (1) estimating intercell properties; (2) modifying interface elevation; (3) calculating dry cell elevations by averaging intercell elevations; and (4) changing the value of the ﬁrst term of slope limiter based on geometry in dry cells. The Harten, Lax, and van Leer with the contact wave restored (HLLC) scheme was implemented to calculate the ﬂux. By combining the MUSCL (Monotone Upstream–centred Scheme for Conservation Laws)-Hancock method with the minmod slope limiter, we achieved second-order accuracy in space and time. This approach is able to keep the conservation property (C-property) and the mass conservation of complex bed geometry. The results of numerical tests in this study are consistent with experimental data, which veriﬁes the effectiveness of the new approach. This method could be applied to acquire wetting and drying processes during ﬂood evolution on structured meshes. Furthermore, a new settlement introduces few modiﬁcation steps, so it could be easily applied to matrix calculations. The new method proposed in this study can facilitate the simulation of ﬂood routing in real terrain.


Introduction
Glacier avalanche [1,2], debris flow [3][4][5], and landslide [6][7][8] in mountain areas could trigger the occurrence of river blocking [9][10][11][12]. Some of this blocking produces large-scale lakes, which leads to back flooding upstream and may inundate roads and villages. Most dammed lakes breach in a short time after their formation, causing massive water to be released catastrophically [9,13]. Yigong Lake was blocked by catastrophic landslides in 1902 and 2000 [14] and formed outburst floods with peak discharges of around 18.9 × 10 4 m 3 /s [15] and 12.4 × 10 4 m 3 /s [8], respectively; the Yarlung Tsangpo gorge was blocked twice in 2018, with a peak discharge of 3.2 × 10 4 m 3 /s in the second outburst flood [3,4]. This kind of dynamic process can impose catastrophic damage to downstream people and infrastructure [16]. Outburst floods may also have significant geomorphic and geologic impacts; they have substantial erosive and transport capacity that can rapidly transform river channels and bedforms [17][18][19], and may even lead to climate change [20] and a global sea level decrease [21]. Outburst floods and their impacts even appear in the myths and stories of many civilizations, such as the Bible and the Koran [22].
Back analysis of outburst flood is an impressive method to determine risk, which has been used to reconstruct large-scale geomorphological dynamic processes that occurred ten thousand years ago. In general, the submerge area and related velocity determine the risk of outburst floods, and a shallow water dynamic model is a widely used and reliable method to predict it [23][24][25][26].
Shallow water equations are popular in long-wave hydrodynamic simulation [27] and are an effective way to analyze outburst flood routing. The Godunov-type finite volume method is an effective and convenient method to calculate flood evolution in complex geometry and is widely used in structured cells and unstructured cells [27]. There are two popular forms for shallow water equations: (1) not consider gravity source term in advection terms [28] and (2) consider the geometry in advection terms [29,30].
A TVD (total variation diminishing) scheme is used to limit numerical oscillations near discontinuity [31][32][33]. Slope limiters such as the minmod limiter, double limiter, and van-Leer limiter are popularly used to keep the solving scheme that has a TVD property [33]. By using a slope limiter, a monotone upstream-centered scheme for conservation laws (MUSCL) reconstruction in the cell center provides second-order accuracy in space [34,35]. The MUSCL method is one of the most successful high-resolution schemes for hyperbolic conservation laws and is applied widely [24,29,33].
Wet-dry front treatment is a key problem when applying shallow water equations to real geometry. Sharp slope geometry especially can over-predict flux and generate negative flow depth [27,29]. Specific treatments during calculation have been applied to limit flux and the gravity source term or to modify geometry [27,29,36], thus or avoiding extremely high flux in intercells and velocity in the cell center. In the process of variable modifications, the limiter's value of the dry cell would equal zero after modifying the local geometry [27,29,36,37].
Many traditional treatments to the wet-dry front change the elevation of the dry cell equal to the wet cell's free surface elevation as shown in Figure 1a [27,29]. If the dry cell is surrounded by four wet cells with different free surface elevations, four elevation modifications are necessary to achieve a balanced flux in the surrounded four cells (Figure 1b,c), and it is very hard to achieve a matrix calculation during simulation as well. A matrix calculation and less cell modification save time because matrix operators are faster than cell loops [38]. In order to apply shallow water equations to a river with a complex geometry and avoid more elevation modifications, we propose processing dry cells by adopting the first term of the slope limiter function in dry cells to solve the wet-dry front problem and accomplish matrix simulation in the whole calculation area. This method can avoid modifications in the dry cell's elevation and achieve a matrix calculation. This method was tested with many cases and is applicable to a complex geometry for outburst flood analysis.

Governing Equations
Two-dimensional shallow water equations are integral forms of Reynolds-averaged Navier-Stokes equations. This equation presumptively neglects vertical momentum exchange and sets the pressure distribution as hydrostatic [39]: where t represents time direction, x and y are two Cartesian coordinates, U is a variable with vector form, F and G are fluxes vectors at two directions, and S is a vector represents source term. The equation is a conserved equation. For general use, the conserved equation is written as: where η = Z + h is the elevation of the flood free surface, where the specific treatment to initial shallow water equations adds geometry information to the advections [29], Z is the elevation of the river bed, h is the flow depth, u is the flow velocity in the x direction, v is the flow velocity in the y direction, τ bx and τ by are the bottom shear stress in the x and y directions, g is gravity acceleration, and n is the Manning coefficient.

Finite Volume Method
The finite volume method has been used in many areas to solve partial equations [40]. The method is implemented by integrating partial equations over the space area for an arbitrary grid. In this study, shallow water equations are hyperbolic equations, which can be integrated as follows: By using Green's formula, Equation (6) can be described as: where L is the mesh boundary of the integral line, and ε is the integral area, which is a rectangular grid here. By using the integral form equation at mesh (i, j), the second term becomes: where n is the time, and i + 1/2 and j + 1/2 are the predicted flux at the interface, predicted by two Riemann states.

HLLC Riemann Solver for Fluxes Prediction
In order to solve the Riemann problem approximately, Harten Lax and van Leer proposed the famous HLL Riemann solver in 1983, which is widely used by researchers to solve shallow water equations today. The scheme requires estimations for the fastest signal velocities from the discontinuity at the interface, resulting in a two-wave model including shock waves, rarefaction waves, and discontinuity. Toro modified the scheme to a three-wave model [33], and the solver was suited to calculate cases involving a wet-dry front, so the HLLC (Harten, Lax and van Leer) approximate Riemann solver by Toro is used in this paper.

Slope Limiter
The face value of variables required for the MUSCL-Hancock reconstruction step and for the time updating step is: where r is the distance vector, and ∇U i is the gradient vector of variable in space. In order to avoid numerical oscillations, we adopt a single slope limiter in this study. The formula becomes: where ϕ(r) is a limiter function. We adopted the Minmod limiter in case tests. Special gradients of variables were predicted by: where r i,j is slope in mesh (i, j), which includes two directions' values. If intercell interpolation is in the x direction, F n = 1 and G n = 0; if intercell interpolation is in the y direction, F n = 0 and G n = 1.

MUSCL-Hancock Method
In the MUSCL-Hancock reconstruction step, the calculation is limited in a single cell. Thus, it does not use the HLLC Riemann solver to predict the flux at the intercell. The corrected value in the cell center is U n+1/2 i , and the flux is calculated based on cell face reconstruction, which is predicted by the cell slope limiter: where U M n i+1/2 is the reconstructed cell boundary vector. The predicted cell center value is calculated by: As for the Riemann flux calculation, we use results from the MUSCL-Hancock step to reconstruct the value around the interface. The slope limiter is the same as the MUSCL-Hancock reconstruction step. The formula is: Riemann states in another direction to use the same method.

Stability Criteria
The numerical scheme is explicit. The stability is defined by the Courant-Friedrichs-Lewy (CFL) criterion. Since this is a two-dimensional calculation case, the time step is limited by local real-time results: where C is the Courant number, ranging between 0 and 1. In some cases, a stable ∆T could give a more stable result. If the export results include a specific time point, ∆T should be modified to a smaller time step to match the predicted time point.

Intercell Bed Elevation and Dry Cell
Since the flux calculation should follow the real physics law in the real world, the interface property determines the flux calculation during flow routing in real river geometry. We classified the interface property into four types based on flow depth and surface elevation (as shown in Figure 2): (1) Two cells' flow depth is higher than 0, which would generate flux in these specific two cells. (2) Two cells between the interface are dry cells such that both flow depths are equal to zero. (3) One is a wet cell and another is a dry cell, but the elevation of the wet cell is higher than the dry cell. (4) One is a wet cell and another is a dry cell, but the dry cell is higher than the wet cell.
Based on the physical property, the interface in the first and third type should consider mass and momentum exchanges between the two cells during calculation. It is not necessary to consider this effect for the cell interface in Type B and Type D. Local modification of Z at the intercell is adopted. The modification is used based on the physical property of the real condition (as shown in Figure 3); e.g., (1) the reflection boundary would stop the flow from moving forward; (2) the dry cell has no flux. The intercell property in Types A, B and C do not need modifications, and the intercell bed elevation is: where Z i+1/2 is the elevation at the intercell; Z i and Z i+1 are cell center elevations at the ith and (i + 1)th cell. Type D of the intercell face's elevation is modified as: In the Type C intercell property, a sharp slope would produce an overpredicted flux in the intercell. Based on the intercell property, the intercell bed elevation was modified as: Momentum needs to be modified while the intercell property is Type D. The velocity component that is perpendicular and the limiter of the three variables of the shallow water equations should be set to zero. For rectangular cell simulation, the calculation area could be treated as a matrix. Many simulations are based on circulation to calculate the whole simulated area, and they include a step that checks for cells that do not need flux calculations. We want to skip this step due to the running circulation cost time. The specific form of the shallow water equation includes η, and the unbalanced flux would be predicted during our simulation which formed by a complex real geometry if the matrix is used directly, for example ( Figure 4): If the dry cell's slope limiter function, Equation (11), is zero, the calculated flux would be unbalanced: In order to achieve a matrix calculation and an automatic flux balance during simulation, we adopted the "zero" slope-limiter function and modified the first term based on the geometry. The elevation of dry cell was modified to: and the slope of the surface elevation of the dry cell was calculated as: where r i,j(η i ) is the value of the first term of the slope limiter function, and ∆x is the cell length in the x direction. If the flow depth in the dry cell is zero, η i+1/2 = Z i+1/2 in the interface, and cell center's value is given by Equation (20). The specific treatment to the dry cell is shown below (as shown in Figure 5): The balance in the dry cell is automatically reached: In the reflection boundary, where a higher left dry cell and a lower right wet cell surround the intercell, η i+1/2 = η i−1/2 = η i and η i−1/2 = Z i−1/2 . The flux balance is reached automatically: If the flow velocity at all described cells is zero, the flux balance is controlled by the wet-dry boundary and the dry cells. All the steps of this method are summarized in Figure 5f.

Steady Condition Calculation of Flood
A test case was used to test the numerical scheme's C-property. A static lake is kept steady, and there is no disturbance. The calculation area is an 8000 m × 8000 m. In the dry bed, there are two bumps: The lake elevation is 1000 m, and the lower bump is submerged by the lake. The mesh size is a rectangular mesh of 1 m × 1 m. The calculation time step is 1 s. The finish time is 8000 s.
After 8000 s, the lake remained static, the results in Figure 6 show that this approach follows a C-property, the static keep balance automatically.

Two-Dimensional Smooth River Bed Test
A two-dimensional smooth bed test was adopted here. The case has an analytical solution smooth bed. This test was adopted by many researchers to test their algorithm's wetdry treatment and calculation accuracy [27,29,41,42]. The calculation area is a 4 m × 4 m, and the origin of the coordinates is in the center of the calculation area. The mesh size is 0.1 m × 0.1 m. The bed is a parabola rotation: where h 0 is the initial flow depth of the origin of the coordinates, a is the distance between the origin and the elevation equal to zero, and x and y are coordinate variables. Under this condition, water flows on the smooth bed and cannot stop. The frequency of flow is ω = 2π/T = 8gh 0 /a, in which T is the time of one cycle. In the analytical solution for the process, the moving range is small: where A = a 4 − r 4 0 / a 4 + r 4 0 , and r 0 is the farthest distance to the center. In the simulation test, we consider the same parameter treatments as Song et al. [42], a = 1 m, h 0 = 0.1 m, and r 0 = 0.8 m. We adopted a mesh size of 0.01 × 0.01 m. The initial condition is the same as the analytical solutions in T/6, T/3, T/2 and T (Figure 7).

Dam Breach over a Thump
This test case is a dam break flow over a thump. The experiment was carried at the University of Brussels, Belgium [43]. Many researchers have used this case to test their model on complex geometries [44,45].
The test simulated a sudden dam breach of flood flowing over a triangular hump. The calculation area is a 38 × 1.75 m flume. A hump was set at 15.5 m, and a barrier lake was formed upstream (as shown in Figure 8). The static lake's flow depth is 0.75 m. The peak of the triangular thump is at 28.5 m, with a height and bottom width of 0.4 and 6 m, respectively. In the tail of the obstacle, there is a 0.15 m high gate, where flow depth is also 0.15 m. Downstream, the first gate is the dry bed. Roughness of the calculation area is n = 0.0125 s × m −1/3 . Four downstream monitoring locations were set, named G1, G2, G3, and G4, and the measured data is the flow depth, located at 19.5, 25.5, 26.5, and 28.5 m respectively. The mesh size for the calculation is 0.1 m × 0.1 m.  Figure 9 shows four representative moments of simulation. After 1 s, the flood front arrives at the 19 m point. At 8 s, the flood flows over the obstacle, which causes backwater and imposes disturbance on the tail lake. At 16 s, a higher run-up upstream lake formed at the front of the obstacle, with waves upstream of the hump. A distinct hydraulic jump develops at the tail lake. At 40 s, the water surface before the obstacle is dominated by strong waves, while the tail lake becomes static. The flow upstream cannot flow over the obstacle. We extracted surface elevation data from the simulation results for comparison. Simulated results at G4 and G13 fit the monitored data very well, but the predicted water surface at G10 and G11 is slightly lower than the monitored data, G20 is slightly higher than measured data, which has been captured in many cases [44]. At the lower stage, the simulated results were similar to simulated results later. The short-term-simulated higher flow depth did not influence the real flood evolution at a later stage. Compared with the same simulated work did by Tomas and Liao [44,45], our simulated results show similar result in G10, G11, and G20. In G4 and G13, our result is closer to measured data compared with their results, which shows better results (as shown in Figure 10). The simulated result is similar to the measured data at G4; (b) initially, the simulated results at G10 is lower but did not influence successive results; (c) the same higher simulated flow depth at G11 is similar to G10, a short-term lower elevation; (d) the simulated results fit well with the measured data at G13; (e) the simulated results fit well with the measured data at G20.

Dam Break Wave Propagating over Three Humps
The three humps test is a very famous test case proposed by Kawahara in 1986 [46,47]. Initially, the case was adopted to test the finite element model, which is wildly used. The calculation area in this study is a 75 × 30 m flume, which has three humps. The boundary is a fixed reflection boundary. The centers of the humps are A (30 m, 6 m), B (30 m, 24 m), and C (47.5 m, 15 m). The maximum height of the humps is 1, 1, and 3 m, respectively. In the upstream of x = 16 m, there is a lake with a depth of 1.875 m. The bed roughness is n = 0.018 sm −1/3 . The calculation geometry was calculated from the formulas below: where a and b are geometric functions of the two lower humps, c is the geometry function of the higher humps, and the elevation of the bed bottom is the maximum value of a, b, and c. The mesh size is 0.5 m × 0.5 m. Figure 11 shows the simulated results of six important moments. At 2 s, the water reached two lower humps and started to flow over them. At 6 s, the flood flowed over the two lower humps and started to reached the higher hump. At 12 s, the flood bypassed the higher hump because it could not completely inundate the higher hump. At 30 s, the flood occupied the calculation area. The formed higher flow depth downstream caused backflow. At 100 s, there was still weak flow in the tank. At 300 s, the flow almost stopped and formed a static lake in the tank, and the peaks of all three humps did not submerge. The numerical model properly simulated complex wetting and drying processes and produced similar results to those of other researchers [29,48].

Conclusions
We propose a new approach to process dry cells and wet-dry front cells via a Godunovtype finite volume prediction method of flood evolution. Shallow water equations automatically balance the gravity source term. The modification includes four steps: (1) identify four types of intercells based on flow depth and surface elevation difference; (2) based on the physical properties of the intercells, modify the bed elevation of the intercell, so as to avoid non-physical flux predictions and gravity balance; (3) modify the dry cell's center elevation to equal the averaged elevation of the two surrounding intercell elevations; (4) change the first term of the slope limiter at the dry cell equal to the ratio of the elevation difference between two intercell bed elevations dividing two times of mesh size. This method was applied to a second-order MUSCL-Hancock-HLLC scheme in time and space for flux and variable prediction in a real geometry. The intercell flux predicted by the reconstructed method remained balanced with the gravity source term automatically, which was proved by mathematical derivations. Four simulated cases showed that the method has a C-property in a complex geometry and achieves the same results as those of many other researchers. Results in the analytical case and the experiment monitoring cases fit each other very well. During all the processing steps, modification could be finished in one step, such that cells did not need to be checked through circulation. This new method can increase the convenience and efficiency of matrix calculations and has a potential for faster GPU (Graphics Processing Unit) simulation and parallel computing. It could be used in real world outburst flood simulation with high efficiency.