*Article* **Numerical Investigation on a Flash Flood Disaster in Streams with Confluence and Bifurcation**

**Qingyuan Yang 1,2, Xiekang Wang 3,\*, Yi Sun 4, Wengang Duan <sup>1</sup> and Shan Xie <sup>5</sup>**


**Abstract:** On 20 August 2019, a flash flood occurred in Sanjiang Town, Sichuan, China, and caused great damage to people living there. The town lies at the junction of five streams, with streams A, B, and C combining at the town and further dividing into streams D and E. The slope of streams A, B, and C is about 3~5%, while the slope of streams D and E is around 0.3%. The Sanjiang Town actually lies in the transition from supercritical slope to subcritical slope. During the flood, huge sediments were released to streams A, B, and C, and further transported to stream E. Due to the rapid change of velocity, only few sediments deposited at the supercritical slope parts of the stream, while plenty of them sedimented at the streams with subcritical slope. In order to simulate the flood with a hydrodynamic model, a field investigation was carried out to collect high DEM (digital elevation model) data, flood marks, sediment grading, etc., after the flood. The discharge curve of the flood was also obtained by the hydrometric station near Sanjiang Town. For the inlet sediment concentrations of streams A, B, and C, we made a series of assumptions and utilized the case which best fits the flood marks to set the inlet sediment concentration. Based on these data, we adopted a depth-averaged two-dimensional hydrodynamic model coupled with a sediment transport model to simulate the flash flood accident. The results revealed that the flash flood enlargement in confluence streams is mainly induced by the inflows, and the flash flood enlargement in bifurcation streams is largely affected by the sediment deposition. The bifurcation of flows can decrease the peak discharge of each branch, but may increase the flooded area near the streams. Flow in the supercritical slope runs at a very fast velocity, and seldom deposits sediment in the steep channel. Meanwhile, most sediment is transported to the streams with flat hydraulic slopes. Due to the functioning of the reservoir, the transition region from supercritical slope to subcritical slope has a much larger probability of being submerged during the flood.

**Keywords:** flash flood; bifurcation; confluence; shallow-water models

#### **1. Introduction**

Flash floods are usually induced by rapidly gathered rain. The flash flood runs down the hill and moves quickly, affects a large region, and has little time for early warning. In mountainous areas, people usually reside in flat areas along the river bank, and facilities such as villages, roads, railway, etc., are all in these areas. The flash flood often leads to great damage to people living nearby.

In order to reduce the threat of flash floods, pre-warning systems are set by many countries. The key point of a warning system is creating a standard for the identification

**Citation:** Yang, Q.; Wang, X.; Sun, Y.; Duan, W.; Xie, S. Numerical Investigation on a Flash Flood Disaster in Streams with Confluence and Bifurcation. *Water* **2022**, *14*, 1646. https://doi.org/10.3390/w14101646

Academic Editor: Marco Franchini

Received: 28 April 2022 Accepted: 19 May 2022 Published: 21 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of disasters, and using the standard to classify the possible impact area of flash flood inundation. These warning systems should respond in a rather short time. In the past, the computation resources were limited, and the modeling tools for flash flood propagation were simplified to save in calculation time. These systems usually adopted some simple models, such as statistical models [1] and simplified hydrodynamic models [2,3] focusing on water levels, discharge, rainfall, etc., but the propagation process of flow was seldom mentioned [4–6]. The warning standard obtained by these models was rough and inaccurate for residential areas with complex channels and buildings. With the development of computers and in order to study the warning standard of flash floods in detail, many researchers have studied flash floods with depth-averaged two-dimensional hydrodynamic models, which are powerful tools to capture flow dynamic behavior and save computational time for a spatially large-scale flow domain [7,8]. These models solve the full governing equations, including the rainfall and infiltration sections. The discretization equations using the Godunov method [9] can accurately capture the shock wave and deal with the sharp change of the bed form. The calculation time is acceptable, and this is a very promising too [7,10–15].

Furthermore, flash flood enlargement, where the peak water level or peak discharge is increased or decreased in areas with rapid bed form obstruction [16] will largely impact the flash flood warning standard. Many studies are mainly focused on the natural issues of the flash floods, such as sediment transport, landslide dams [17], confluence of rivers, etc. Yang et al. [18] studied the impact of sediment transport on the peak level and peak discharge of flash floods with numerical models, and Chen et al [4] investigated the effect of landslide dams on the flash flood peak level with a large-scale physical model. Wang et al. [19] made progress in the flash flood propagation in the river confluence zone with both numerical and experimental techniques. Chen et al [20] analyzed the flash flood wave propagation in the confluence of open channels. Hackney et al. [21] made a dent in discharge variation of the flash flood in the river confluence region with field investigations.

Generally speaking, flash floods in mountainous areas occur in channels with very steep slopes, and huge sediments can be transported far away due to the high-speed flow. However, if there is a reservoir setting on the stream, the slope will be slowed by the highly raised water level. A transition region from the supercritical slope to the subcritical slope will be produced by the sedimentation. Flash floods in streams with confluence and bifurcation, especially in the transition region from supercritical slope to subcritical slope, are seldom studied in detail. In this study, we use a depth-average two-dimensional hydrodynamic model to conduct a numerical inversion of a flash flood that occurred on 20 August 2019 in Sanjiang Town, Sichuan, China.

#### **2. Study Area**

The study area covers Sanjiang Town, which lies in Sichuan Province of China. Figure 1 shows the satellite image of Sanjiang Town. The town lies at the confluence of three streams, A, B, and C. Firstly, stream B meets stream C at the upper part of the town, then they both join stream A. The distance between the two junctions is about 100 m. Right after the confluence, stream A is further divided into streams D and E. For stream D, there is a milldam near the bifurcation, and the milldam can discharge part of the flash flood through stream D. Stream E has twice the width of stream D, and releases most of the flood. Streams D and E meet again at the Sanjiang reservoir. The central bar surrounded by streams D and E is the major developing part of the town in recent years. At the end of stream E, a hydrometric station collects hydraulic data during the floods. An overflow dam lies at the end of the Sanjiang reservoir, and gates of the dam are all opened to discharge floods.

**Figure 1.** Satellite image of Sanjing Town, Sichuan, China.

#### **3. Methodology**

In this study, we adopted the depth-averaged, two-dimensional, shallow-water equations, coupled with sediment transport and bed variation equations [22], to simulate the flash flood and sediment transport during the disaster.

Continuity equation:

$$\frac{\partial}{\partial t}(h) + \frac{\partial}{\partial x}(hu) + \frac{\partial}{\partial y}(hv) = -\frac{\partial Z\_{\emptyset}}{\partial t} \tag{1}$$

Momentum equation:

$$\frac{\partial}{\partial t}(hu) + \frac{\partial}{\partial x}\left(hu^2 + \frac{1}{2}gh^2\right) + \frac{\partial}{\partial y}(huv) = gh\left(S\_{bx} - S\_{fx}\right) + hv\_l\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) - \frac{\Lambda \rho g h^2}{2\rho\_l \rho\_m}\frac{\partial S\_T}{\partial x} + \frac{\rho\_0 - \rho\_m}{\rho\_m}u\frac{\partial Z\_b}{\partial t} \tag{2}$$

$$\frac{\partial}{\partial t}(hv) + \frac{\partial}{\partial y}(huv) + \frac{\partial}{\partial y}\left(hv^2 + \frac{1}{2}gh^2\right) = gh\left(S\_{by} - S\_{fy}\right) + hv\_l\left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) - \frac{\Delta \rho g h^2}{2\rho\_l \rho\_m}\frac{\partial S\_T}{\partial y} + \frac{\rho\_0 - \rho\_m}{\rho\_m}v\frac{\partial Z\_b}{\partial t} \tag{3}$$

Equation of bed load transport:

$$\frac{\partial}{\partial t}(hq\_b) + \frac{\partial}{\partial \mathbf{x}}(huq\_b) + \frac{\partial}{\partial y}(hvq\_b) = -a\_b\omega\_b(q\_b - q\_{b\*}) = -\rho' \frac{\Delta Z\_b}{\Delta t} \tag{4}$$

Equation of bed deformation:

$$\frac{\Delta Z\_b}{\Delta t} = \frac{a\_b \omega\_b (q\_b - q\_{b\*})}{\rho'} \tag{5}$$

in which *t* denotes time, *h* means water depth, *u* and *v* present velocity components in the x and y directions, respectively, *g* is gravitational acceleration, *vt* is the turbulent viscosity coefficient, and Δ*ρ* = *ρ<sup>s</sup>* − *ρw*, in which *ρ<sup>s</sup>* is the sediment density, and *ρ<sup>w</sup>* is the water density.

$$
\rho\_{\rm w} = (1 - S\_{\rm v})\rho\_{\rm w} + S\_{\rm v}\rho\_{\rm s} = \left(1 - \frac{S}{\rho\_{\rm s}}\right)\rho\_{\rm w} + S\_{\rm T} \tag{6}
$$

where *Sv* is the volumetric sediment concentration, and *ST* is the total concentration of graded sediments (kg/m3). In this study, only the bed load transport is calculated, so *ST* equals *qb*.

$$
\rho\_0 = \left(1 - \frac{\rho'}{\rho\_s}\right) \rho\_w + \rho' \tag{7}
$$

in which *ρ*<sup>0</sup> and *ρ* denote the density of saturated and dry bed material, respectively. The bed slope terms (*Sbx*, *Sby*) and friction slope terms (*Sf x*, *Sf y*) are written as *Sbx* = −*∂Zb*/*∂x*, *Sby* = −*∂Zb*/*∂y*, and *Sf x* = *<sup>n</sup>*2*<sup>u</sup>* <sup>√</sup>*u*<sup>2</sup> <sup>+</sup> *<sup>v</sup>*2/*h*4/3, *Sf y* <sup>=</sup> *<sup>n</sup>*2*<sup>v</sup>* <sup>√</sup>*u*<sup>2</sup> <sup>+</sup> *<sup>v</sup>*2/*h*4/3 in the x and y directions, respectively, where *Zb* is bed elevation and *n* is Manning's roughness coefficient.

$$q\_{b0} = \frac{K\_b}{C\_0^2} \frac{\rho\_s \rho\_m}{\rho\_s - \rho\_m} (\mathcal{U} - \mathcal{U}\_c) \frac{\mathcal{U}^3}{\mathcal{g}\omega\_b} \tag{8}$$

$$q\_{b\*} = q\_{b0}/(h\mathcal{U}) = \frac{\mathcal{K}\_b}{\mathcal{C}\_0^2} \frac{\rho\_{\mathcal{S}} \rho\_m}{\rho\_{\mathcal{S}} - \rho\_m} (\mathcal{U} - \mathcal{U}\_{\mathcal{E}}) \frac{\mathcal{U}^2}{h \mathcal{G} \omega \omega\_b} \tag{9}$$

where *qb* in kg/m<sup>3</sup> is the amount of bed load in a unit volume of water, *ω<sup>b</sup>* is the settling velocity of bed load, *qb*<sup>0</sup> is the transport capacity of bed load in a unit volume of water, in kg/m3, and *α<sup>b</sup>* is the non-equilibrium adaptation coefficient of bed load. *Kb* is an empirical coefficient, *Uc* is the incipient velocity of bed load, calculated by *Uc* = 0.265 log(11*h*/*d*) (*ρs*−*ρw*) *<sup>ρ</sup><sup>w</sup> gd*, *<sup>C</sup>*<sup>0</sup> is the dimensionless Chézy coefficient, and *qb*<sup>∗</sup> is the value of bed load in a unit volume, and is obtained according to *qb*<sup>∗</sup> = *qb*0/(*hU*).

The governing Equations (1)–(4) can be rewritten in the form:

$$\frac{\partial \Pi}{\partial t} + \frac{\partial E}{\partial x} + \frac{\partial G}{\partial y} = S \tag{10}$$

in which

$$\begin{aligned} \,^t U = \begin{bmatrix} h \\ hu \\ hv \\ hq\_b \end{bmatrix}, \,^t E = \begin{bmatrix} h \\ hu \\ hv \\ hq\_b \end{bmatrix}, \,^G G = \begin{bmatrix} h \\ hu \\ hv \\ hq\_b \end{bmatrix} \end{aligned} $$

*S* is the source section, and is equal to the rest of the sections of each governing equation.

$$\mathcal{U}\_{i}^{n+1} = \mathcal{U}\_{i}^{n} + \sum (En\_{x} + Gu\_{y}) + S \tag{11}$$

The governing equations are solved by the method of finite volume with hybrid triangular and quadrangular meshes. The variables (*h*, *u*, *v*, *qb*, *Zb*) in each mesh are updated in a time-marching manner (Equation (11)). *U<sup>n</sup> <sup>i</sup>* is the variables' value at time *<sup>n</sup>*, *<sup>U</sup>n*+<sup>1</sup> *<sup>i</sup>* is the variables' value at time *n* + 1, and *Enx* + *Gny* denotes the flux crossing the edges of the mesh. In order save the calculation time, a NVIDIA GeForce RTX 3050 GPU (graphics processing unit) was adopted to speed up the calculation. The flow chart of the calculation is shown in Figure 2.

The calculation codes used in this study have been successfully used in the former studies [18,23]. The simulation region covers most parts of Sanjiang Town, and the number of simulation meshes is about 110,000. In order to capture the flash flood propagation accurately with a moderate quantity of mesh, a non-uniform (triangular and quadrangular hybrid) mesh is adopted. The size is about 2.0 m in regions near the stream channel and bank, and about 5.0 m in regions far from the streams (Figure 3).

We collected high-resolution (5.0 m) DEM (digital elevation model) data (Figure 1) after the flash flood in 2019, and the flow rates of stream E (Qt(CS14)) at the hydrometric station of Sanjiang Town. The flow rates at the hydrometric station were gauged at an interval of 20 min. The flash flood lasted about 120 h (from 20 to 24 August in 2019). The 120 h-long flash flood can be divided into four stages, and each of them has one discharge peak (Figure 4a, stage 1, 2, 3, and 4, respectively). We deduced the discharge of streams A (*QA*), B (*QB*), and C (*QC*) based on their upstream catchment area (Figure 4b). For example, since the milldam in stream D was closed during the flood, the total discharge of streams A, B, and C should be equal to the discharge at CS14. Supposing the catchment areas of streams A, B, and C are *AA*, *AB*, and *AC*, respectively, and the total area of

them *AT* = *AA* + *AB* + *AC*, then *QA* = Qt*AA*/*AT*, *QB* = Qt*AB*/*AT*, and *QC* = Qt*AC*/*AT*, respectively. All these data contribute to the numerical inversion of the flash flood.

**Figure 2.** Flow chart of the calculation.

**Figure 3.** Recorded cross-sections (top) and part of the calculation meshes (bottom).

**Figure 4.** Monitored (**a**) and deduced (**b**) discharges of the flash flood.

In this study, the roughness of the bed was set with the Manning's roughness, *n*, of 0.025, which is suggested by hydraulics manuals [24] and verified by our former studies [18]. The bed elevations at each mesh were initiated by interpolating from the high-resolution terrain data, and updated by the simulation codes at each timestep based on Equation (5). In order to study the unsteady process of the flash flood, hydraulic data of 18 cross-sections (Figure 3) were saved at every timestep, and data of the full simulation region were saved at every 100 timesteps during the calculation.

The inlet discharges for streams A, B, and C were consistent with the deduced rate curves in Figure 4b, and flow rate of the outlet boundary was automatically calculated according to the empirical formula for weir flow. The formula is *Q* = *mnb*2*gH*1.5, in which *Q* = discharge, *m* = discharge coefficient, about 0.502, *n* = number of gates opened, *b* = width of each gate, and *H* = the difference between water level and elevation of the weir top. Once the discharge is calculated, the velocity at the outlet is set based on the discharge.

The milldam was closed during the flash flood in 2019 (Figure 5), and released little water to stream D. The water in streams A, B, C, and E performs as the confluence flow, but the flow will move as bifurcation flow when the milldam is open. Therefore, in order to investigate the impact of the milldam on the flash flood, scenarios with both a closed and an opened milldam were simulated in this study.

**Figure 5.** Image of the milldam after the flash flood in 2019.

As the sediment characteristics were not monitored by the hydrometric station during the flash flood disaster, it is difficult to figure out the amount of bed load at inlet boundaries of streams A, B, and C. We assumed the inlet sediment concentrations with different percentages of transport capacity of bed load at the inlets. Transport capacity will change with velocity and depth. For example, if the transport capacity of bed load at the inlets at time A is 100 kg/m3, Qs-0.05 means the inlet sediment concentration is set as <sup>100</sup> × 0.05 = 5 kg/m<sup>3</sup> at time A, and Qs-0 means the inlet sediment concentration is zero all the time. For the diameters of sediment, we collected some sediment of the river bed after the flash flood in 2019, and obtained a median size of 0.02 m. The flow conditions for the simulation scenarios are listed in Table 1.


**Table 1.** Flow conditions for the simulation scenarios.

#### **4. Results**

#### *4.1. Model Calibration*

The purpose of this section is to figure out the scenario which best fits the monitored data. We analyzed the simulation results of the first five scenarios in Table 1. In these scenarios, the milldam in stream D was closed during the flood.

Figure 6a shows the simulated and monitored discharge profiles at CS14, and Figure 6b presents a detailed view of stage 3. It can be seen that the simulated discharge profile agrees well with the monitored data at CS14, except for that of case Qs-0.125 in stage 3. Figure 7 shows the distribution of unit width flux of cases Qs-0.125 and Qs-0.1. The flooded area of case Qs-0.125 was a bit larger than that of case Qs-0.1, and some flow travelled downstream through CS17 instead of CS14, which contributed to the bigger difference in Figure 6b.

**Figure 6.** Simulated and monitored discharge profiles at CS14: (**a**) whole flood and (**b**) stage 3 of the flood.

**Figure 7.** Flooded region of cases Qs-0.1 (**a**) and Qs-0.125 (**b**).

Figure 8a shows the simulated water level of each scenario. For cases Qs-0 and Qs-0.05, the maximum level occurred in stage 2 of the flood, while for cases Qs-0.075, Qs-0.1, and Qs-0.125, the maximum level occurred in stage 3. Figure 8b shows the simulated bed elevation of each scenario. The bed elevation increased rapidly in flood stages 2 and 3, and the larger the inlet sediment concentration, the more the bed elevation increased. The lifted bed elevation further contributed to the occurrence time of the maximum level being delayed from flood stage 2 to 3. The sediment deposition can change the amplitude and occurrence time of the peak level.

**Figure 8.** Simulated water level (**a**) and bed elevation (**b**) at CS14.

According to the analysis above, the discharge curves were nearly coincident with each other among the first four simulation scenarios. However, water level and bed elevation were much different from each other for cases with different inlet sediment concentrations. The different water level and bed elevation further impacted the flooded area during the flood. Based on the field investigation carried out right after the flash flood in 2019, the simulation results of case Qs-0.1 fit better concerning the flood marks, final bed elevation, flood area, etc. Therefore, we used the 10% of transport capacity (Qs-0.1) as the inlet sediment concentration.

#### *4.2. Flash Flood in Confluence Streams (Case Qs-0.1)*

As the gates of the milldam were closed in this case, all of the flow coming from streams A, B, and C went to the reservoir through stream E. Figure 9 presents the peak discharge along the streams. In the four stages of the flash flood, the peak discharges increased rapidly at CS11, and showed little change at other parts of the streams.

**Figure 9.** Peak discharge along the streams (Case Qs-0.1, (**a**) stage 1 and 2, (**b**) stage 3 and 4).

Figure 10a shows the peak level at each cross-section. The water level dropped sharply along the channel of streams A, B, and C, with a minimum slope of 2.5%. On the contrary, in stream E, the peak water level slowly decreased along the channel, and the slope was about 0.3%. Once again, this proves that the town lies in the transition region from a steep channel to a gentle slope.

Peak water levels were generally coincident with the flood marks obtained by field investigations, and it was difficult to compare them as the slope was steep (Figure 10a). Figure 10b presents the relative peak water level at all the cross-sections. The relative peak water level was obtained by subtracting the peak water level of each flood stage from that of stage 1, so the relative peak water levels of stage 1 were all zero. For cross-sections in streams A, B, and C, the maximum level occurred in stage 2, in which the maximum peak discharge existed. However, for cross-sections in stream E, the maximum level occurred in stage 3 of the flood. The bed elevation in stage 3 was much higher than that of stage 2, and the deposition lifted the water to a higher level with a relatively smaller discharge (Figure 8).

For a better understanding of the flow properties along the streams, we adopted the unit width flux to analyze the main flow variation. As the unit width flux, *q*, is equal to the product of depth and velocity magnitude, regions with larger *q* denote more flow passing through. Figure 11 shows the distribution of unit width flux (*q*) at the four stages of the

flood. Since the milldam was closed, the main flow came from streams A and B, and went to the reservoir only through stream E.

**Figure 10.** Peak water level (**a**) and relative peak water level (**b**) (Case Qs-0.1).

Siltation thickness denotes the differences between bed elevation at any time and that of the initial time. The larger the thickness, the more sediments deposited there. Figure 12 presents the distribution of siltation thickness at the four stages of the flood. Sediments deposited mainly at the confluence region of streams A and B in stages 1 and 2, while the main siltation region shifted to stream E in stages 3 and 4. Meanwhile, the siltation thickness in streams A and B became thinner in stages 3 and 4, which explains the finding in the field investigation whereby the flood marks were much higher but the elevation of sediments was much lower in stream A. In stages 1 and 2, the water level in streams A and B was higher due to the larger discharge and higher bed elevation, while in stages 3 and 4, the water level in streams A and B was lower due to the smaller discharge and lower bed elevation, due to erosion. After the flood, only the final bed elevation and flood marks could be measured.

**Figure 11.** Distribution of unit width flux (*q*) at the four stages of the flood ((**a**) stage 1, (**b**) stage 2, (**c**) stage 3, (**d**) stage 4).

**Figure 12.** Distribution of siltation thickness at the four stages of the flood ((**a**) stage 1, (**b**) stage 2, (**c**) stage 3, (**d**) stage 4).

Figure 13 shows the distribution of the Froude number (Fr) at the four stages of the flood. The Froude number is defined as the velocity magnitude divided by water depth (sqrt(uˆ2+uˆ 2)/h). If Fr > 1, the flow is subcritical, if Fr > 1, the flow is supercritical, and if Fr = 1, the flow is critical. For streams A, B, and C, the Froude number of most regions was larger than 3 in the four flood stages, while for stream E, the flow was subcritical in stages 1 and 2, and supercritical in stages 3 and 4. There is little siltation at places with a larger Froude number. On the contrary, for stream E, the regions with a larger Froude number were nearly coincident with those of higher sediment thickness (Figure 12). The highly raised bed decreased the water depth and increased the Froude number.

**Figure 13.** Distribution of the Froude number at the four stages of the flood ((**a**) stage 1, (**b**) stage 2, (**c**) stage 3, (**d**) stage 4).

#### *4.3. Flash Flood in Streams with Bifurcation (Case Qs-0.1-Bifurcation)*

As the gates of the milldam were open in this case (Qs-0.1-bifurcation), all of the flow coming from streams A, B, and C travelled to the reservoir through both stream D and stream E. Figure 14 presents the peak discharge along the streams. The peak discharge increased rapidly at CS11 and CS18 due to the combination of flows, decreased sharply after CS11, and showed little change at other parts of the stream.

**Figure 14.** Peak discharge along the streams (Qs-0.1-bifurcation, flood (**a**) stage 1, stage 2. (**b**) stage 3, stage 4).

Figure 15 shows the peak level at each cross-section. The water level varied in a similar way to that in Figure 10. The maximum water level of streams D and E occurred in flood stage 3 (the second largest discharge stage), while the maximum water level of streams A, B, and C occurred in flood stage 2 (the largest discharge stage).

**Figure 15.** Peak water level (**a**) and relative peak water level (**b**) (Qs-0.1-bifurcation).

Figure 16 shows the distribution of unit width flux (q) at the four stages of the flood. The main flow came from streams A and B, and travelled to the reservoir through stream E in stages 1 and 2, while more and more flow travelled to the reservoir via stream D in stages 3 and 4. Figure 17 presents the flow rates and flow rate ratios of streams D and E. In stages 1 and 2 (before 30 h), the flow rate ratio of stream D was about 0.1, while in stages 3 and 4, the flow rate ratio of stream D increased to around 0.3. Figure 18 presents the distribution of siltation thickness at the four stages of the flood as the milldam was open. Sediments deposited mainly at the confluence region of streams A, B, and C in stages 1 and 2. The siltation region shifted to streams D and E in stages 3 and 4, and siltation in stream E increased much quicker than that in stream D. The uplifted bed of stream E pushed more water to stream D. For the distribution of Froude numbers in Figure 19, supercritical regions were still coincident with places of thick sediment deposition in streams D and E, and with little sediment deposits at places with a higher Froude number in streams A, B, and C.

**Figure 16.** Distribution of unit width flux (q) at the four stages of the flood ((**a**) stage 1, (**b**) stage 2, (**c**) stage 3, (**d**) stage 4).

**Figure 17.** Discharge (**a**) and discharge ratios (**b**) between streams D and E.

**Figure 18.** *Cont*.

**Figure 18.** Distribution of siltation thickness at the four stages of the flood ((**a**)stage 1, (**b**) stage 2, (**c**) stage 3, (**d**) stage 4).

**Figure 19.** *Cont*.

**Figure 19.** Distribution of the Froude number at the four stages of the flood ((**a**)stage 1, (**b**) stage 2, (**c**) stage 3, (**d**) stage 4).

Figure 20 shows the flooded area of cases Qs-0.1 and Qs-0.1-bifurcation. Among all the four stages of the flood, the flooded area of case Qs-0.1-bifurcation was larger than that of case Qs-0.1. The flooded area ratio between case Qs-0.1-bifurcation and case Qs-0.1 is shown in Figure 20b, and the ratio had a maximum of 1.1. From the point of view of reducing the flooded area, it is better to keep the milldam closed during the flood.

**Figure 20.** Flooded area of cases Qs-0.1 (**a**) and Qs-0.1-bifurcation (**b**).

#### **5. Discussion**

#### *5.1. Boundary Conditions for the Numerical Inversion of Flash Flood*

In this study, we adopted the two-dimensional shallow-water models coupled with sediment transport models to analyze the flash flood that occurred in Sanjiang Town in 2019. On most occasions, flash floods occur in a region with limited monitoring facilities. It is difficult to accurately set the boundaries for flash flood simulations. Usually, the elevation of the bed form, the grade curve of sediment, and flood marks can be measured after the flood disaster. However, the time-discharge curve, the inlet sediment concentration, and the outlet level are rather difficult to estimate. Luckily, we obtained the monitored discharge curves from the hydrometric station near Sanjiang Town, and deduced the discharge for streams A, B, and C by their upstream catchment areas. For the inlet boundary, we made assumptions with a series of sediment concentrations, and used the sediment concentration that best fit the field investigation data. Theoretically, the method is more accurate than those using a constant inlet sediment concentration, although it is still not certain that the sediment concentration is accurate due to the complexity of sediment transport. The method of evaluating the boundaries in this paper is helpful for the numerical analysis of other flash flood disasters, and hydrological models coupled with hydrodynamic models may be taken into account if the discharge curve cannot be directly obtained [2].

Furthermore, the quantitative simulation results are still helpful for taking measures against the flood. For example, keeping the milldam closed can minimize the flooding area.

#### *5.2. Impact of Sediment on Flash Flood in Stream Confluence and Bifurcation*

Flash flood wave enlargement is common in flash flood propagation, and it is usually caused by the sediment intrusion and bed form obstruction [18]. For sediment intrusion, the enlargement comes from the bulk effect, and the volume and density of fluid are both enlarged. For the bed form obstruction, the enlargement is rooted in the sharp change of the water level and velocity of the bed form. For example, the landslide dams, the change in the bed slope, etc., can increase or decrease the peak level and discharge.

The confluence and bifurcation are prototypes of streams. The flash flood enlargement in confluence streams is impacted not only by the bed form, such as the junction angle, elevation differences, etc., but also by the flow rate ratio of the streams [19]. In this study, as the flow rate ratio among streams A, B, and C was constant, the flash flood enlargement in confluence streams was mainly induced by the merging of inflows. Meanwhile, the flash flood enlargement in bifurcation streams is largely affected by the sediment deposition. The flow rate ratio between the two branches was changed by the uneven sedimentation of streams D and E. The adjustment in the flow rate ratio is similar to that in plain rivers, except that the variation in mountainous streams is much quicker [8].

The distinctive part of this study is that the study region lies in the transition region from a supercritical slope to a subcritical slope. The hydraulic slopes of streams D and E were flattened by the Sangjiang reservoir. Huge sediments deposited in streams D and E as the velocity was slowed down. From the point of view of reducing the flooded area, the position of Sangjiang reservoir was not well-planned considering the sediment transport mode, and it should be moved downstream, far away from Sanjiang Town.

#### **6. Conclusions**

(1) In this study, the flash flood that occurred on 20 August 2019 in Sanjiang Town was analyzed by a depth-averaged two-dimensional model. The simulation showed the flash flood process in detail, and introduced a method of evaluating the boundary conditions for the numerical inversion of flash floods.

(2) Mountainous streams with many people living nearby are usually reformed by manmade facilities, which will affect the flash flood propagation properties. In this study, the milldam in stream D could largely reduce the flooded area. The manmade facilities should be taken into account in the study of flash floods, and be properly utilized during the

flood events. Furthermore, planning of reservoirs should consider the sediment transport properties, as the hydraulic slope will be flattened by the reservoir.

(3) For streams in the transition region from supercritical slope to subcritical slope, the flash flood enlargement in confluence streams is mainly induced by the combination of inflows, and the flash flood enlargement in streams with bifurcation is largely affected by the sediment deposition. Flow in the supercritical slope runs at a very fast velocity, and seldom deposits sediment in the steep channel. In the meantime, the sediment is transported to the streams with a flat hydraulic slope and deposited there. The transition region has a much larger probability of being submerged during the flood.

**Author Contributions:** Conceptualization, Q.Y. and X.W.; investigation, X.W.; resources, W.D.; data curation, Y.S. and S.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key R&D Program of China (Grant No. 2019YFC1510703-03), the National Natural Science Foundation of China (Grant Nos. 51639007, 51979007, 51609014, and 51809107), and the Fundamental Research Funds for Central Public Welfare Research Institutes (Grant No. CKSF2021482/SL).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Written informed consent has been obtained from the patient(s) to publish this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Glossary of Terms**


#### **References**

