Next Article in Journal
Evapotranspiration Response to Climate Change in Semi-Arid Areas: Using Random Forest as Multi-Model Ensemble Method
Next Article in Special Issue
Experimental Study on Landslides in Terraced Fields in the Chinese Loessial Region under Extreme Rainfall
Previous Article in Journal
Spatio-Temporal Analysis of Surface Water Quality in Mokopane Area, Limpopo, South Africa
Previous Article in Special Issue
A Study on Interaction between Overfall Types and Scour at Bridge Piers with a Moving-Bed Experiment
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Key Laboratory of Mountain Hazards and Earth Surface Process, Chinese Academy of Sciences, Chengdu 610041, China
2
Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu 610041, China
3
University of Chinese Academy of Sciences, Beijing 100049, China
4
Key Laboratory of Land Surface Pattern and Simulation, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China
5
The Engineering & Technical College of Chengdu University of Technology, Leshan 614000, China
6
Sichuan Highway Planning, Survey, Design and Research Institute Ltd., Chengdu 610041, China
*
Author to whom correspondence should be addressed.
Water 2021, 13(2), 221; https://doi.org/10.3390/w13020221
Submission received: 9 November 2020 / Revised: 12 January 2021 / Accepted: 14 January 2021 / Published: 18 January 2021
(This article belongs to the Special Issue Soil–Water Conservation, Erosion, and Landslide)

Abstract

:
When utilizing a finite volume method to predict outburst flood 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 first 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 flux. 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 verifies the effectiveness of the new approach. This method could be applied to acquire wetting and drying processes during flood evolution on structured meshes. Furthermore, a new settlement introduces few modification steps, so it could be easily applied to matrix calculations. The new method proposed in this study can facilitate the simulation of flood routing in real terrain.

1. 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 × 104 m3/s [15] and 12.4 × 104 m3/s [8], respectively; the Yarlung Tsangpo gorge was blocked twice in 2018, with a peak discharge of 3.2 × 104 m3/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.

2. Governing Equations and Schemes

2.1. 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]:
U , t + F , x + G , y = S ,
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:
η h u h v , t + h u h u 2 + g ( η 2 2 η Z ) / 2 h u v , x + h u h u v h v 2 + g ( η 2 2 η Z ) / 2 , y = 0 τ b x / ρ g η Z , x τ b y / ρ g η Z , y ,
τ b x = ρ g n 2 u u 2 + v 2 h 1 / 3 ,
τ b y =   ρ g n 2 v u 2 + v 2 h 1 / 3 ,
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, τ b x and τ b y are the bottom shear stress in the x and y directions, g is gravity acceleration, and n is the Manning coefficient.

2.2. 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:
t ε U d Ω + ε ( F x + G y ) d Ω = ε S d Ω .
By using Green’s formula, Equation (6) can be described as:
t ε U d Ω + L ( F + G ) d L = ε S d Ω ,
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:
L F i d L + L G j d L = ( F i + 1 / 2 F i 1 / 2 ) Δ y + ( G j + 1 / 2 G j 1 / 2 ) Δ x ,
U i , j n + 1 = U i , j n Δ t Δ x F i + 1 / 2 , j F i 1 / 2 , j Δ t Δ y G i , j + 1 / 2 G i , j 1 / 2 + Δ t S i ,
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.

2.3. 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.

2.4. Slope Limiter

The face value of variables required for the MUSCL-Hancock reconstruction step and for the time updating step is:
U i + 1 / 2 = U i + r U i ,
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:
U i + 1 / 2 = U i + φ r r U i ,
where φ r is a limiter function. We adopted the Minmod limiter in case tests. Special gradients of variables were predicted by:
r i , j = η i + F n , j + G n η i , j η i , j η i F n , j G n h u i + F n , j + G n h u i , j h u i , j h u i F n , j G n h v i + F n , j + G n h v i , j h v i , j h v i F n , j G n ,
where ri,j is slope in mesh (i, j), which includes two directions’ values. If intercell interpolation is in the x direction, Fn = 1 and Gn = 0; if intercell interpolation is in the y direction, Fn = 0 and Gn = 1.

2.5. 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 i n + 1 / 2 , and the flux is calculated based on cell face reconstruction, which is predicted by the cell slope limiter:
U M i + 1 / 2 n = U i n + 1 2 φ ( r ) U i n U i 1 n ,
where U M i + 1 / 2 n is the reconstructed cell boundary vector. The predicted cell center value is calculated by:
U i t + 1 / 2 = U i t + k x F U M i + 1 / 2 n F i + 1 / 2 U M i 1 / 2 n + k y G U M j + 1 / 2 n G U M j 1 / 2 n + Δ t 2 S i
k x = Δ t 2 Δ x ;   k y = Δ t 2 Δ y .
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:
U i + 1 / 2 L = U i n + 1 / 2 + 1 2 φ ( r ) U i n U i 1 n .
Riemann states in another direction to use the same method.

2.6. 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:
Δ t = m i n C Δ x u i + g h i , C Δ y v i + g h i , Δ T ,
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.

3. 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:
Z i + 1 / 2 = Z i + Z i + 1 / 2 ,
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:
Z i + 1 / 2 = m i n η i , η i + 1 .
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:
Z i + 1 / 2 = m a x Z i , Z i + 1 .
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:
g ( η 2 2 η Z ) / 2 , x = g η i 2 2 η i Z i 1 / 2 η i 2 2 η i Z i + 1 / 2 2 Δ x g η i Z , x .
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:
η i = Z i + 1 / 2 + Z i 1 / 2 2 ,
and the slope of the surface elevation of the dry cell was calculated as:
r i , j η i = Z i + 1 / 2 Z i 1 / 2 2 Δ x ,
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:
g ( η 2 2 η Z ) / 2 , x = g Z i + 1 / 2 2 Z i 1 / 2 2 2 Δ x = g Z i + 1 / 2 + Z i 1 / 2 Z i + 1 / 2 Z i 1 / 2 2 Δ x = g η i Z , x
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:
g ( η 2 2 η Z ) / 2 , x = g η i 1 / 2 2 + η i + 1 / 2 2 + 2 η i + 1 / 2 Z i + 1 / 2 2 η i 1 / 2 Z i 1 / 2 2 Δ x = g η i + 1 / 2 Z i 1 / 2 Z i + 1 / 2 Δ x = g η i Z , x .
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.

4. Results and Discussion

4.1. 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:
Z x , y = m a x 0 , Z B 1 , Z B 2 ,
Z B 1 = 2000 0.00032 x 3000 2 + y 5000 2 Z B 2 = 900 0.000144 x 5000 2 + y 3000 2 .
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.

4.2. 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 wet-dry 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:
Z x , y = h 0 x 2 + y 2 a 2 1 ,
where h0 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 = 8 g h 0 / a , in which T is the time of one cycle. In the analytical solution for the process, the moving range is small:
η x , y , t = m a x Z x , y , h 0 1 A 2 1 A c o s ω t x 2 + y 2 a 2 1 A 2 1 A c o s ω t 2 1 1 ,
where A = a 4 r 0 4 / a 4 + r 0 4 , and r0 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, h0 = 0.1 m, and r0 = 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).

4.3. 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).

4.4. 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:
a = 1 1 8 x 30 2 + y 6 2 b = 1 1 8 x 30 2 + y 24 2 c = 3 3 10 x 47.5 2 + y 15 2 Z x , y = m a x 0 , a , b , c ,
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].

5. Conclusions

We propose a new approach to process dry cells and wet-dry front cells via a Godunov-type 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.

Author Contributions

D.L. (Dingzhu Liu) and H.W. design this research and draft the manuscript. J.T. gave many suggestions on testing cases. Y.C. help with coding the program. N.A.B. help with the language. H.C. and D.L. (Daochuan Liu) help with find checking data. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant no. 41941017; the Applied Fundamental Research Program of Sichuan Province, grant no. 2019JY0387; Key Research Program of Frontier Sciences, Chinese Academy of Sciences, grant no. QYZDY-SSW-DQC006 and the National Natural Science Foundation of China.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study is contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Capps, D.M.; Clague, J.J. Evolution of glacier-dammed lakes through space and time; Brady Glacier, Alaska, USA. Geomorphology 2014, 210, 59–70. [Google Scholar] [CrossRef] [Green Version]
  2. Cook, S.J.; Kougkoulos, I.; Edwards, L.A.; Dortch, J.M.; Hoffmann, D. Glacier change and glacial lake outburst flood risk in the Bolivian Andes. Cryosphere 2016, 10, 2399–2413. [Google Scholar] [CrossRef] [Green Version]
  3. Chen, C.; Zhang, L.; Xiao, T.; He, J. Barrier lake bursting and flood routing in the Yarlung Tsangpo Grand Canyon in October 2018. J. Hydrol. 2020, 583. Available online: https://www.sciencedirect.com/science/article/abs/pii/S0022169420300639 (accessed on 16 January 2021). [CrossRef]
  4. Hu, K.; Zhang, X.; You, Y.; Hu, X.; Liu, W.; Li, Y. Landslides and dammed lakes triggered by the 2017 Ms6.9 Milin earthquake in the Tsangpo gorge. Landslides 2019, 16, 993–1001. [Google Scholar] [CrossRef]
  5. Wei, R.; Zeng, Q.; Davies, T.; Yuan, G.; Wang, K.; Xue, X.; Yin, Q. Geohazard cascade and mechanism of large debris flows in Tianmo gully, SE Tibetan Plateau and implications to hazard monitoring. Eng. Geol. 2018, 233, 172–182. [Google Scholar] [CrossRef]
  6. Fan, X.; Xu, Q.; Alonsorodriguez, A.; Subramanian, S.S.; Li, W.; Zheng, G.; Dong, X.; Huang, R. Successive landsliding and damming of the Jinsha River in eastern Tibet, China: Prime investigation, early warning, and emergency response. Landslides 2019, 16, 1003–1020. [Google Scholar] [CrossRef]
  7. Li, B.; Feng, Z.; Wang, G.; Wang, W. Processes and behaviors of block topple avalanches resulting from carbonate slope failures due to underground mining. Environ. Earth. Sci. 2016, 75, 694. [Google Scholar] [CrossRef]
  8. Liu, W.; Carling, P.A.; Hu, K.; Wang, H.; Zhou, Z.; Zhou, L.; Liu, D.; Lai, Z.; Zhang, X. Outburst floods in China: A review. Earth Sci. Rev. 2019, 197, 102895. [Google Scholar] [CrossRef]
  9. Cui, P.; Dang, C.; Zhuang, J.Q.; You, Y.; Chen, X.Q.; Scott, K.M. Landslide-dammed lake at Tangjiashan, Sichuan province, China (triggered by the Wenchuan Earthquake, May 12, 2008): Risk assessment, mitigation strategy, and lessons learned. Environ. Earth. Sci. 2012, 65, 1055–1065. [Google Scholar] [CrossRef]
  10. Wang, G.Q.; Fan, L. Simulation of dam breach development for emergency treatment of the Tangjiashan Quake Lake in China. Sci. China Ser. E Technol. Sci. 2008, 51, 82–94. [Google Scholar] [CrossRef]
  11. Yan, Y.; Cui, Y.F.; Xin, T.; Hu, S.; Guo, J.; Wang, Z.; Yin, S.Y.; Liao, L.F. Seismic Signal Recognition and Interpretation of the 2019 “7.23” Shuicheng Landslide by Seismogram Stations. Landslides 2020, 17, 1206–1911. [Google Scholar] [CrossRef]
  12. Zhang, L.; Xiao, T.; He, J.; Chen, C. Erosion-based analysis of breaching of Baige landslide dams on the Jinsha River, China, in 2018. Landslides 2019, 16, 1965–1979. [Google Scholar] [CrossRef]
  13. Costa, J.E.; Schuster, R.L. The formation and failure of natural dams. Geol. Soc. Am. Bull. 1988, 100, 1054–1068. [Google Scholar] [CrossRef]
  14. Zhou, J.; Cui, P.; Hao, M. Comprehensive analyses of the initiation and entrainment processes of the 2000 Yigong catastrophic landslide in Tibet, China. Landslides 2016, 13, 39–54. [Google Scholar] [CrossRef]
  15. Ma, D.T. Study on Influences of Mountain Hazards in Yigong Zangbu River Basin to Mitigation and Reconstruction of Sichuan-Tibetan Highway Line. Ph.D. Thesis, The Graduate School of Chinese Academy of Sciences, Beijing, China, 2006. [Google Scholar]
  16. Cui, P.; Zhu, Y.Y.; Han, Y.S.; Chen, X.Q.; Zhuang, J.Q. The 12 May Wenchuan earthquake-induced landslide lakes: Distribution and preliminary risk evaluation. Landslides 2009, 6, 209–223. [Google Scholar] [CrossRef]
  17. Carling, P.A. Freshwater megaflood sedimentation: What can we learn about generic processes? Earth Sci. Rev. 2013, 125, 87–113. [Google Scholar] [CrossRef]
  18. Carling, P.A.; Fan, X. Particle comminution defines megaflood and superflood energetics. Earth Sci. Rev. 2020, 204, 103087. [Google Scholar] [CrossRef]
  19. Turzewski, M.D.; Huntington, K.W.; LeVeque, R.J. The Geomorphic Impact of Outburst Floods: Integrating Observations and Numerical Simulations of the 2000 Yigong Flood, Eastern Himalaya. J. Geophys. Res Earth. 2019, 124, 1056–1079. [Google Scholar] [CrossRef]
  20. Teller, J.T.; Leverington, D.W.; Mann, J.D. Freshwater outbursts to the oceans from glacial Lake Agassiz and their role in climate change during the last deglaciation. Quat. Sci. Rev. 2002, 21, 1–887. [Google Scholar] [CrossRef]
  21. Garcia-Castellanos, D.; Estrada, F.; Jimenez-Munt, I.; Gorini, C.; Fernandez, M.; Verges, J.; De Vicente, R. Catastrophic flood of the Mediterranean after the Messinian salinity crisis. Nature. 2009, 462, 778–781. [Google Scholar] [CrossRef]
  22. Burr, D.M.; Carling, P.A.; Baker, V.R. Megaflooding on Earth and Mars; Cambridge University Press: Cambridge, UK, 2009. [Google Scholar]
  23. Anacona, P.I.; Mackintosh, A.; Norton, K. Reconstruction of a glacial lake outburst flood (GLOF) in the Engano Valley, Chilean Patagonia: Lessons for GLOF Risk management. Sci. Total Environ. 2015, 527–528, 1–11. [Google Scholar] [CrossRef] [PubMed]
  24. Bohorquez, P.; Cañada-Pereira, P.; Jimenez-Ruiz, P.J.; del Moral-Erencia, J.D. The fascination of a shallow-water theory for the formation of megaflood-scale dunes and antidunes. Earth Sci. Rev. 2019, 193, 91–108. [Google Scholar] [CrossRef]
  25. George, D.L. Adaptive finite volume methods with well-balanced Riemann solvers for modeling floods in rugged terrain: Application to the Malpasset dam-break flood (France, 1959). Int. J. Numer. Methods Fluids 2011, 66, 1000–1018. [Google Scholar] [CrossRef]
  26. Swartenbroekx, C.; Zech, Y.; Soares-Frazão, S. Two-dimensional two-layer shallow water model for dam break flows with significant bed load transport. Int. J. Numer. Methods Fluids 2013, 73, 477–508. [Google Scholar] [CrossRef]
  27. Hou, J.; Liang, Q.; Simons, F.; Hinkelmann, R. A 2D well-balanced shallow flow model for unstructured grids with novel slope source term treatment. Adv. Water Resour. 2013, 52, 107–131. [Google Scholar] [CrossRef]
  28. Ma, D.J.; Sun, D.J.; Yin, X.Y. Solution of the 2D shallow water equations with source terms in surface elevation splitting form. Int. J. Numer. Methods Fluids 2007, 55, 431–454. [Google Scholar] [CrossRef]
  29. Liang, Q.; Borthwick, A.G.L. Adaptive quadtree simulation of shallow flows with wet-dry fronts over complex topography. Comput. Fluids 2009, 38, 221–234. [Google Scholar] [CrossRef]
  30. Rogers, B.; And, M.F.; Borthwick, A.G.L. Adaptive Q-tree Godunov-type scheme for shallow water equations. Int. J. Numer. Methods Fluids 2001, 3, 247–280. [Google Scholar] [CrossRef]
  31. Moukalled, F.; Mangani, L.; Darwish, M. The Finite Volume Method in Computational Fluid Dynamics. An Advanced Introduction with OpenFoam® and Matlab®; Springer: Berlin, Germany, 2016. [Google Scholar]
  32. Toro, E.F. Shock-Capturing Methods for Free-Surface Shallow Flows; Wiley: Hoboken, NJ, USA, 2001. [Google Scholar]
  33. Toro, E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction; Springer: Berlin, Germany, 2013. [Google Scholar]
  34. Van Leer, B. Towards the ultimate conservative difference scheme, V: Asecond-order sequel to Godunov’s method. J. Comput. Phys. 1979, 32, 101–136. [Google Scholar] [CrossRef]
  35. Harten, A.; Lax, P.D.; van Leer, B. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev. 1983, 25, 35–61. [Google Scholar] [CrossRef]
  36. Cea, L.; Puertas, J.; Vazquezcendon, M. Depth Averaged Modelling of Turbulent Shallow Water Flow with Wet-Dry Fronts. Arch. Comput. Method E 2007, 14, 303–341. [Google Scholar] [CrossRef]
  37. Hu, P.; Cao, Z.; Pender, G.; Tan, G. Numerical modelling of turbidity currents in the Xiaolangdi reservoir, Yellow River, China. J. Hydrol. 2012, 464, 41–53. [Google Scholar] [CrossRef]
  38. Liu, L.; Liu, L.; Yang, G.W. Cache performance optimization of irregular sparse matrix multiplication on modern multi-core CPU and GPU. High Technol. Lett. 2013, 339–345. [Google Scholar] [CrossRef]
  39. Wu, W. Computational River Dynamics; Crc Press: Boca Raton, FL, USA, 2007. [Google Scholar]
  40. Cao, Z.; Pender, G.; Wallis, S.; Carling, P. Computational dam-break hydraulics over erodible sediment bed. J. Hydraul. Eng. 2004, 130, 689–703. [Google Scholar] [CrossRef]
  41. Liang, Q. Flood Simulation Using a Well-Balanced Shallow Flow Model. J. Hydraul. Eng. 2010, 136, 669–675. [Google Scholar] [CrossRef]
  42. Song, L.; Zhou, J.; Li, Q.Q.; Yang, X.L.; Zhang, Y.C. An unstructured finite volume model for dam-break floods with wet/dry fronts over complex topography. Int. J. Numer. Methods Fluids 2011, 67, 960–980. [Google Scholar] [CrossRef]
  43. Hiver, J. Adverse-slope and slope (bump). In Proceedings of the Concerted Action on Dam Break Modelling: Objectives, Project Report, Test Cases, Civil Engineering Department, Hydraulic Division, Universitè Catholique de Lille, Lille, France, 26–28 June 2000. [Google Scholar]
  44. Liao, C.; Wu, M.S.; Liang, S.J. Numerical simulation of a dam break for an actual river terrain environment. Hydrol. Processes 2007, 21, 447–460. [Google Scholar] [CrossRef]
  45. Rebollo, T.C.; Nieto, E.D.; Marmol, M.G. A flux-splitting solver for shallow water equations with source terms. Int. J. Numer. Methods Fluids 2003, 42, 23–55. [Google Scholar] [CrossRef]
  46. Zhou, J.G.; Causon, D.M.; Minghan, C.G. Numerical prediction of dam-break flows in general geometries with complex bed topography. J. Hydraul. Eng. 2004, 130, 332–340. [Google Scholar] [CrossRef]
  47. Kawahara, M.; Umetsu, T. Finite element method for moving boundary problems in river flow. Int. J. Numer. Methods Fluids 1986, 6, 365–386. [Google Scholar] [CrossRef]
  48. Brufau, P.; Vázquez-Cendón, M.E.; García-Navarro, P. A numerical model for the flooding and drying of irregular domains. Int. J. Numer. Methods Fluids 2002, 39, 247–275. [Google Scholar] [CrossRef]
Figure 1. The traditional elevation modification of wet-dry front. (a) Modify elevation to the same as wet cell; (b,c) Two times elevation modification of one dry cell.
Figure 1. The traditional elevation modification of wet-dry front. (a) Modify elevation to the same as wet cell; (b,c) Two times elevation modification of one dry cell.
Water 13 00221 g001
Figure 2. Classification of interface property. (a) Type A: wet cells at the left and right side, hL > 0, hR > 0, hL and hR are flow depth in the left and right side of intercell respectively; (b) Type B: dry cells at the left and right side of the intercell face; (c) Type C: wet and dry cells are connected through the intercell face, the free surface elevation of the wet cell is higher than the dry cell; (d) Type D: wet and dry cells are connected between the intercell face, and the free surface elevation of the wet cell is lower than the dry cell.
Figure 2. Classification of interface property. (a) Type A: wet cells at the left and right side, hL > 0, hR > 0, hL and hR are flow depth in the left and right side of intercell respectively; (b) Type B: dry cells at the left and right side of the intercell face; (c) Type C: wet and dry cells are connected through the intercell face, the free surface elevation of the wet cell is higher than the dry cell; (d) Type D: wet and dry cells are connected between the intercell face, and the free surface elevation of the wet cell is lower than the dry cell.
Water 13 00221 g002
Figure 3. Modification of the intercell elevation. (a,b) The intercell does not need modification, which is related to Type A and Type B; (c) the intercell elevation is modified to the wet cell’s elevation, which is related to Type D; (d) the sharp slope cell is modified to the dry cell’s bed elevation.
Figure 3. Modification of the intercell elevation. (a,b) The intercell does not need modification, which is related to Type A and Type B; (c) the intercell elevation is modified to the wet cell’s elevation, which is related to Type D; (d) the sharp slope cell is modified to the dry cell’s bed elevation.
Water 13 00221 g003
Figure 4. The calculated parameters of a shallow water equation with no special treatment.
Figure 4. The calculated parameters of a shallow water equation with no special treatment.
Water 13 00221 g004
Figure 5. The dry cell’s center elevation is calculated by the average of two intercell elevations. Intercell elevations are predicted from the latest two steps that are based on the intercell type and the real conditions. (a,c) One side is a wet-dry front and the one side is dry-dry; (b) both sides are a wet-dry front; (d) both sides are sharp slopes; (e) both sides are dry cells; (f) flow chart of the method.
Figure 5. The dry cell’s center elevation is calculated by the average of two intercell elevations. Intercell elevations are predicted from the latest two steps that are based on the intercell type and the real conditions. (a,c) One side is a wet-dry front and the one side is dry-dry; (b) both sides are a wet-dry front; (d) both sides are sharp slopes; (e) both sides are dry cells; (f) flow chart of the method.
Water 13 00221 g005
Figure 6. C-property checking for a static lake. (a) Lake geometry; (b) results after 8000 s.
Figure 6. C-property checking for a static lake. (a) Lake geometry; (b) results after 8000 s.
Water 13 00221 g006
Figure 7. Simulated results compared with real analytical results. (a) Geometry of the calculation area and the initial condition; (be) comparison between the simulated results and the analytical solution at T/6, T/3, T/2, T.
Figure 7. Simulated results compared with real analytical results. (a) Geometry of the calculation area and the initial condition; (be) comparison between the simulated results and the analytical solution at T/6, T/3, T/2, T.
Water 13 00221 g007aWater 13 00221 g007b
Figure 8. Flume test setup of the experiment.
Figure 8. Flume test setup of the experiment.
Water 13 00221 g008
Figure 9. Free surface elevation during the flood evolution in the experiment. (a) At 1 s, flood flows at the dry bed; (b) at 8 s, the flood flows up to the obstacle and has an influence downstream; (c) at 16 s, all the upstream water flows to the obstacle and a run-up forms; (d) at 40 s, the flow downstream remains static, with waves at the upstream lake.
Figure 9. Free surface elevation during the flood evolution in the experiment. (a) At 1 s, flood flows at the dry bed; (b) at 8 s, the flood flows up to the obstacle and has an influence downstream; (c) at 16 s, all the upstream water flows to the obstacle and a run-up forms; (d) at 40 s, the flow downstream remains static, with waves at the upstream lake.
Water 13 00221 g009
Figure 10. Monitored data compared with the simulated results at the four locations. (a) 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.
Figure 10. Monitored data compared with the simulated results at the four locations. (a) 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.
Water 13 00221 g010
Figure 11. Simulated flood evolution on a complex three-hump condition. (a) The flood starts to reach the first two low humps at 2 s; (b) the flood flows over the two low humps at 6 s; (c) the flood flows downstream of the high humps at 12 s; (d) the flood forms a higher flow depth downstream at 30 s; (e) there is some weak flow in the tank at 100 s; (f) the tank maintains a static condition at 300 s.
Figure 11. Simulated flood evolution on a complex three-hump condition. (a) The flood starts to reach the first two low humps at 2 s; (b) the flood flows over the two low humps at 6 s; (c) the flood flows downstream of the high humps at 12 s; (d) the flood forms a higher flow depth downstream at 30 s; (e) there is some weak flow in the tank at 100 s; (f) the tank maintains a static condition at 300 s.
Water 13 00221 g011aWater 13 00221 g011b
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Liu, D.; Tang, J.; Wang, H.; Cao, Y.; Bazai, N.A.; Chen, H.; Liu, D. A New Method for Wet-Dry Front Treatment in Outburst Flood Simulation. Water 2021, 13, 221. https://doi.org/10.3390/w13020221

AMA Style

Liu D, Tang J, Wang H, Cao Y, Bazai NA, Chen H, Liu D. A New Method for Wet-Dry Front Treatment in Outburst Flood Simulation. Water. 2021; 13(2):221. https://doi.org/10.3390/w13020221

Chicago/Turabian Style

Liu, Dingzhu, Jinbo Tang, Hao Wang, Yang Cao, Nazir Ahmed Bazai, Huayong Chen, and Daochuan Liu. 2021. "A New Method for Wet-Dry Front Treatment in Outburst Flood Simulation" Water 13, no. 2: 221. https://doi.org/10.3390/w13020221

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop