**Impacts of Filled Check Dams with Di**ff**erent Deployment Strategies on the Flood and Sediment Transport Processes in a Loess Plateau Catchment**

#### **Honglei Tang, Hailong Pan and Qihua Ran \***

Institute of Hydrology and Water Resources, Zhejiang University; Hangzhou 310058, China; tanghonglei@zju.edu.cn (H.T.); panhailong@zju.edu.cn (H.P.)

**\*** Correspondence: ranqihua@zju.edu.cn

Received: 2 April 2020; Accepted: 3 May 2020; Published: 7 May 2020

**Abstract:** As one of the most widespread engineering structures for conserving water and soil, check dams have significantly modified the local landform and hydrologic responses. However, the influences of sedimentary lands caused by filled up check dams on the runoff and sediment transport processes were seldom studied. Employing an integrated hydrologic-response and sediment transport model, this study investigated the influences of filled check dams with different deployment strategies in a Loess Plateau catchment. Six hypothetical deployment strategies of check dams were compared with no-dam scenario and the reality scenario. Results showed that filled check dams were still able to reduce Flood peak (Qp) by 31% to 93% under different deployment strategies. Considerable delays of peak time and decreases were also found in scenarios, which were characterized as having larger and more connective sedimentary lands on the main channel. Reduction rates of Sediment yield (SY) and the total mass of Eroded sediment (ES) ranged from 4% to 52% and 2% to 16%, respectively, indicating that proper distributions of check dams can promote sediment deposition in the channel and reduce soil erosion. The results of this study indicate that (1) check dam systems could still be useful in flood attenuation and sediment control even when they were filled, and (2) optimizing the deployment strategies of check dams can help reduce erosion.

**Keywords:** check dam system; sedimentary land; flood control; sediment transport; Loess Plateau

#### **1. Introduction**

Building check dams is an effective engineering measure for flood attenuation in mountainous streams that experience torrents [1,2] and soil erosion control in arid and semiarid regions with large sediment yields [3–5]. Because of their high efficiency in sediment reduction and relatively low building investment, check dams are constructed in many countries, such as Mexico [6], Spain [7], Italy [8], Iran [9], and China [10]. The Loess Plateau in Northwest China has suffered severe soil erosion due to its combination of highly-erodible soil layers and considerable human activities [11]. The erosion rate in the region ranged from 2000 to 20,000 t km−<sup>2</sup> yr−<sup>1</sup> in the 1950s, resulting in the riverbed uplift of the downstream Yellow River [12]. To reduce the large catchment sediment yield caused by soil erosion, more than 110,000 check dams have been constructed in the Loess Plateau catchments over the past 50 years (especially in the 1970s) to prevent sediment from being transported to the Yellow River [13]. Significant decreases in annual runoff and sediment yields of different catchments due to the construction of check dams have been widely observed and investigated in the Loess Plateau [14–16]. However, most of the check dams have been filled quickly due to the large amount of sediment discharge during floods, and a large number of check dams were destroyed due to low design criteria and unscientific site location [17].

As an engineering structure chronically changing the landforms of a natural channel, the influences of filled check dams on the flood and sediment transport processes have seldom been quantitatively studied. In the early stage of a check dam (i.e., the design storage capacity has not been consumed), channel flow and sediment can be easily blocked by the dam body. When a check dam is filled, the block effect of dam body no longer exists, but a large sedimentary land is formed behind the check dam. The sedimentary land is normally a flat platform with relatively higher soil water content and much lower slope gradient, which is widely used as farmlands in the Loess Plateau. The hydrological impacts of sedimentary land, such as infiltration promotion [18], evaporation increase [19], and enhancement of groundwater recharge [20], have been reported. According to the low slope gradient of sedimentary land and its impacts on hydrologic response, it may influence the flood and sediment transport processes in the long term [21], even as the check dam is filled. In fact, approximately 320,000 ha of sedimentary lands behind check dams were formed during the 1950s to 2000s in the Loess Plateau [13], and the area of sedimentary lands are expected to continue increasing as more and more check dams are about to be fully filled. Considering that all check dams are going to be filled one day and the sedimentary lands formed behind check dams will exist in a long term, it is necessary to quantitatively study their influences on the flood and sediment transport processes.

On the other hand, improving the performance of check dams in flood attenuation and sediment control by optimizing the deployment strategies for treated catchments is still a challenge, given the fact that many of the check dams built in the past 50 years were prone to dam break [17]. A recent study has paid attention to the optimization of check dams in a Loess Plateau catchment [22]. Employing a distributed erosion model (WaTEM/SEDEM) based on the Revised Universal Soil Loss Equation, the life expectancy of check dams under different deployment strategies and their impacts on sediment delivery ratio of the Shejiagou catchment were compared in [22]. Several insights on extending the life expectancy of check dams and simultaneously keeping the system effective in sediment reduction were obtained in this simulation-based study. However, limited by the deficiency of WaTEM/SEDEM, which does not simulate the flow and sediment movement on the sedimentary land, the hydrologic responses of sedimentary land such as infiltration improvement and flow-velocity decrease were not considered in this study. Besides, including the studies introduced above, most of the studies on the optimization of check dam systems focused on improving the sediment trapping efficiency of the system during its lifespan [22–24]. Few studies considered how to optimize the system in advance so as to achieve a better performance of sedimentary lands in flood and sediment control when all the check dams are filled.

Driven by these twofold challenges, the study reported herein investigated the influences of filled check dams with different deployment strategies on the flood and sediment transport processes. Following the work in [22], the Integrated Hydrology Model (InHM), which is fully physics-based, was employed to the 4.26 km<sup>2</sup> Shejiagou catchment. The objectives were (1) to quantitatively understand whether and how much filled check dams with their sedimentary lands could attenuate flood processes and reduce sediment output, and (2) to compare the influences caused by filled check dams under different deployment strategies on the flood and sediment transport processes.

#### **2. Data and Methods**

#### *2.1. Study Site and Data Sources*

The Shejiagou catchment (109◦58 E, 37◦42 N) is a small catchment belonging to Wudinghe Basin in the Loess Plateau, China (Figure 1a). The catchment is 4.26 km2 with an elevation ranging from 920 m to 1143 m (22.8◦ mean slope, Figure 1b). The catchment is covered by a thick loess layer with an average depth of over 100 m. Being highly eroded, the hillslopes of the catchment are deeply incised by gullies. The land cover of the catchment was mainly sparse grassland, bare land, and cultivated land in the 1960s (Figure 1c) [25,26]. The Shejiagou catchment has a semiarid continental climate with a mean annual precipitation of 500 mm and a mean annual potential evaporation of 1570 mm [27]. Seventy percent of the precipitation falls in the wet seasons from June to September. Infiltration-excess surface flow dominates runoff generation during floods in the catchment, and most eroded sediment is carried downstream by surface flow in heavy storms. There is normally no flow in the stream channels and gullies during inter-storm periods. The mean annual specified sediment yield was 12,800 t km−<sup>2</sup> yr−<sup>1</sup> based on records from 1959 to 1969 (2100 t km−<sup>2</sup> yr−<sup>1</sup> to 24,400 t km−<sup>2</sup> yr<sup>−</sup>1). According to the Strahler stream order, the river network is relatively simple and mainly consists of one third-order stream (i.e., the main channel) and several first-order streams (i.e., gullies), which symmetrically located in the two sides of the main channel (Figure 1d). Because the left-side gullies were more incised, the left side of the catchment has a steeper slope distribution and contributes more sediment than the right side. According to the model simulation in [22], major erosion areas distributed on hillslopes with slope gradients greater than 30◦ (Figure 1b), and the headwater areas of the catchment and the fifth gully on the right side of the main channel (Figure 1d) contribute the largest amount of sediment (see the Supplementary Materials).

**Figure 1.** Information of the Shejiagou catchment: (**a**) location and Digital Elevation Model (DEM); (**b**) slope distribution; (**c**) land use types; (**d**) stream network classification and check dam distribution.

Table 1 summarizes the data set used in this study. The 20-m resolution Digital Elevation Model (DEM) of Shejiagou catchment was obtained from the State Bureau of Surveying and Mapping. Infiltration experiments and soil sampling were conducted by the Zizhou Experimental Station in Tuanshangou, a sub-catchment of Shejiagou, to obtain data of soil water characteristics during 1960–1962. By using these soil data, the parameters of van Genuchten Model [28], which was used to describe the soil-water retention and permeability functions, were determined. Runoff, sediment, and precipitation data during floods were obtained from the Experimental Station of Soil and Water Conservation of Zizhou County, which belongs to the Yellow River Conservancy Commission of the Ministry of Water Resources (YRCC). The Shejiagou hydrologic station was set at the outlet of Shejiagou catchment to monitor the runoff and sediment data including water level (m), discharge (m3 s<sup>−</sup>1), sediment concentration (kg m−3), and sediment discharge (kg s−1) during 1959–1969. During this monitoring period, the hydrographs and sedigraphs of 31 rainfall-runoff events were recorded. The time intervals of the flow and sediment discharge data ranged from 5 min around flow peaks to 1 h for low-flow periods. According to three field surveys of check dams, four check dams were constructed sequentially from 1957 to 1973, holding a total storage capacity of 465,100 m3 (Figure 1d). The Shejiagou catchment was chosen in this study mainly because (1) the data and information of the catchment were sufficient for model construction and (2) the land cover of the Shejiagou catchment was relatively simple in the 1960s, which was beneficial to model parameterization.


**Table 1.** Data set used in this study.

#### *2.2. The Integrated Hydrology Model (InHM)*

Based on the blueprint of the distributed physics-based hydrological model proposed by Freeze and Harlan [31], InHM was developed to quantitatively simulate, via a fully coupled approach, 3-D variably saturated flow in soil and 2-D flow and sediment transport across the land surface [32–34]. The 3-D Richards' equation was implemented to describe variably saturated flow in soil, while the 2-D diffusion-wave equation coupled with depth-integrated multiple-species sediment transport was applied to describe surface flow movement and sediment transportation. They were discretized in space using the control-volume finite element method and coupled in one coherent framework using a physics-based first-order flux relationship driven by pressure head gradients [35]. Newton iteration was used to implicitly solve each coupled system of nonlinear equations. More details of InHM can be found in the Appendix A.

As a research tool with robust, general and efficient solution methods [36], InHM can run event-based or continuous simulations in various catchments with different types of topography and climate [37]. Compared to many process-based models [38–41], InHM has the following advantages in this study. (1) The time step is adaptive, and the output information is spatially distributed. This feature enables the model to represent vital variables (e.g., flow velocity, water depth, and sediment concentration) in a certain area (e.g., the sedimentary land) and at a specified time step. (2) The triangle spatial elements used in InHM are flexible and accurate to construct different shapes of sedimentary lands under various check dam deployment strategies. Hydrologically active areas such as gully mouths and tail-end areas can be refined to capture small variations in flow and sediment. The model is more dynamic compared to the methods in which the check dam is generalized as a sediment-trapping function [22] and more physics-based compared to the approaches that calculate flow and sediment output of a check dam by a simple water depth-discharge relationship [42].

The model has been utilized to examine the environmental impacts of many engineering structures, such as dam removal [43] and forest roads [44]. In the Loess Plateau catchment of Shejiagou, InHM is capable of simulating rainfall-runoff processes dominated by the infiltration-excess overland flow mechanism and rainsplash erosion as well as hydraulic erosion processes in flood events [26]. These successful applications indicated that InHM is a competent tool to study the impacts of filled check dams under different deployment strategies in the Shejiagou catchment.

#### *2.3. Modelling Scenario*

Three variables, namely, number, size, and site location of check dams, were usually considered to compare the runoff-sediment reduction efficiencies of different deployment strategies [23]. Apart from increasing number and/or size of check dams, which inevitably increases the deployment costs, appropriate choice of site locations may also improve the sediment reduction rate of a check dam system [22]. To study the impacts on runoff and sediment processes of filled check dam systems under different deployment strategies, three steps were conducted to set up modelling scenarios: (1) add different check dam systems into the model, (2) fill the check dams via rising the elevations of the channels influenced by each check dam to represent a filled check dam, and (3) conduct rainfall-runoff events in each modelling scenarios.

In this study, eight scenarios were designed (Figure 2, Table 2). Scenario-0 and Scenario-7 were cases which represent the no-check dam situation and the current situation of Shejiagou catchment, respectively. Scenarios 1-6 were cases considering different numbers, sizes and site locations of check dams. Although in all scenarios check dams were all filled (i.e., the design storage capacity was used up), the sizes still mattered because the sedimentary land formed by the check dam was different between small check dams and large check dams. For example, the total areas of sedimentary lands in scenario-2 and scenario-3 were 11.21 <sup>×</sup> 104 <sup>m</sup><sup>2</sup> and 5.93 <sup>×</sup> 104 m2, respectively. As for the determination of site locations, whether the check dams were constructed on gullies (the first or second order streams) or on the main channel (i.e., the third order stream) was considered in this study. For example, the second check dam was constructed on a first-order stream, which receives a high sediment load (2000–3600 t yr−1, see the Supplementary Materials) in scenario-5, while the second check dam was on the main channel which receives a much higher cumulative sediment load (9000–37,521 t yr<sup>−</sup>1) in scenario-6. The six scenarios (i.e., scenarios 1–6) were conceived to compare the effects of deployment strategies with different number (e.g., scenarios 1, 2, 4), size (e.g., scenarios 5, 6), and site location (e.g., scenarios 2, 3) of check dams. The six scenarios were also compared with scenario-0 and scenario-7, which provides indications of how the performance of a filled check dam system could be improved.

**Figure 2.** Check dam systems under different deployment strategies. Note that scenario-7 is the current deployment strategy of Shejiagou catchment.


**Table 2.** Description of the eight scenarios in the study.

#### *2.4. Model Settings and Parameters*

#### 2.4.1. Finite-Element Mesh

Eight 3-D meshes for the scenarios were constructed by adding subsurface layers to 2-D triangular surface meshes. The major differences of the eight 3-D meshes were the different combination of check dams and the sedimentary lands formed behind check dams. Figure 3a shows the 3-D mesh that contained 10,062 nodes and 19,918 triangular elements in scenario-1, in which only one filled check dam was considered. The cell dimensions of all the surface meshes varied from 180 m near the headwater boundary to 20 m along the gullies, and the meshes of the check dam were refined to 3 to 10 m. Sixteen subsurface meshes were added under the surface mesh with varying thickness from 0.10 m near the surface to 39.30 m near the bedrock with a constant basal elevation of 900 m (see Figure 3b). This setting of mesh size ensured a fine resolution in the hydrologically active areas (i.e., the channel, sedimentary land, and the near-surface soil) and simultaneously saved computation resources in the relatively inactive areas (i.e., the headwater regions and deep spaces). Figure 3c shows the check dam (i.e., S1-1) and its sedimentary land in the finite-element mesh. To model a filled check dam, the sedimentary land was elevated to the same elevation of the spillway, indicating that the dam body can no longer block surface flow. Based on the fact that the sedimentary lands normally have low slope gradients (usually less than 2%) when compared to the original channel [45,46], it was assumed that, for each filled check dam, the nodes on the sedimentary land have the same elevations to simplify the model construction.

#### 2.4.2. Boundary Conditions and Initial Conditions

Two subsurface boundary conditions were assigned to each 3D mesh: (1) impermeable for each lateral face and the bottom face (A-E-D-C-B and F-J-I-H-G in Figure 3a); and (2) a local sink (i.e., head-dependent flux) at the down-gradient face [47], subsurface flowing out through A-B-G-F in Figure 3b. The surface boundary conditions were driven by the given precipitation to the whole catchment and constrained by a critical water depth boundary condition (no dam at the outlet, so the critical water depth = 0 m) at the catchment outlet (A–B in Figure 3b). A one-hour uniform rainfall with an intensity of 90 mm hr−<sup>1</sup> was chosen to conduct the rainfall-runoff simulations in this study because this rainfall intensity was common in the Shejiagou catchment to induce floods, according to the 31 recorded rainfall-runoff events during 1959–1969. The evapotranspiration process was not considered in these storm event simulations.

**Figure 3.** 3D mesh used in scenario-1: (**a**) the whole catchment and boundary conditions; (**b**) the outlet area and the outflow boundary; (**c**) the surface mesh representing the filled check dam (S1-1) and its sedimentary land. Note that the elevation of sedimentary land was the same as that of the spillway on the right side of the check dam, indicating the check dam was filled.

Before all the scenario runs, a drainage run with no rainfall input was started to drive the whole catchment from a fully saturated state to a realistic state, and it stopped when the simulated discharge at the catchment outlet became stable and matched the average no-rainfall discharge of the catchment at the beginning of wet seasons (i.e., 0.006 m s<sup>−</sup>1). This drainage simulation was performed to obtain a continuous initial distribution of the hydraulic head that was close to the reality of the catchment. Using this initial condition, a one-week warming-up stage with 0.75 mm hr−<sup>1</sup> rainfall was conducted at the beginning of each scenario run, and then the one-hour 90 mm hr−<sup>1</sup> rainfall for the scenario run was input to the model. Each simulation ended 42 h after the rainfall, lasting for 2 days total.

#### 2.4.3. Model Parameters

A near-surface loess layer (0–0.5 m) and a deep loess layer (0.5 m to the bedrock) constituted the soil profile in the model. The near-surface soil layer was classified into three types, i.e., the normal surface loess soil (labeled as "SJG-soil (0–0.5 m)"), the silted soil of the sedimentary land (labeled as "Sedimentary land"), and the compacted soil of the check dam (labeled as "Check dam").

Table 3 summarizes the values of parameters that represent the soil hydraulic characteristics. The saturated hydraulic conductivities (*Ksat*) of soils were calibrated (see Section 2.4.4), while the porosity values and van Genuchten parameters of the soils were derived from the experiment data by the Zizhou Experimental Station. Table 4 summarizes the land surface parameters representing the hydraulic characteristics and erosion/deposition capabilities. For the hydrologic-response module in InHM, Manning's roughness coefficients (*n*) were obtained by calibration (see Section 2.4.4). For the sediment-transport module, the rainsplash coefficient (*cf*) and surface erodibility coefficient (ϕ), the two important parameters in InHM to describe the rainsplash erosion and hydraulic erosion processes, were determined by model calibration (see Section 2.4.4). Derived from previous studies based on soil texture [48,49], the damping coefficient (*cd*), rainfall intensity exponent (*b*), and raindrop turbulence factor (ξ) were set to 600 m−1, 1.6, and 0.25, respectively. According to the soil sample data in the Shejiagou catchment, the median diameter of the soil was set to 0.05 mm to represent the uppermost homogeneous loess soil for the sediment transport simulation. The soil cohesion coefficient was 0.30 [43] for all surface soils except that of the check dam body, which was assigned a larger cohesion coefficient (i.e., 0.60) due to compaction.


**Table 3.** Soil parameters used in InHM for the eight scenarios.

<sup>a</sup> Saturated hydraulic conductivity.; <sup>b</sup> Parameter related to the inverse of the air-entry pressure.; <sup>c</sup> Parameter related to the pore-size distribution.; <sup>d</sup> Residual soil-water content [28].


**Table 4.** Surface parameters used in InHM for the four scenarios.

<sup>a</sup> Manning's roughness coefficient. <sup>b</sup> Immobile water depth, identical for all surface zones. <sup>c</sup> Surface residual saturation, identical for all surface zones. <sup>d</sup> Average height of non-discretized micro-topography, identical for all surface zones. <sup>e</sup> Rainsplash depth dampening factor, identical for all surface zones. <sup>f</sup> Rain intensity exponent, identical for all surface zones. <sup>g</sup> Rain-induced turbulence coefficient, identical for all surface zones. <sup>h</sup> Rainsplash coefficient. <sup>i</sup> Surface erodibility coefficient.

#### 2.4.4. Model Calibration and Validation

Two rainfall-runoff events were selected for the model calibration (event of 1964/7/14) and validation (event of 1964/8/2). The two events were chosen among the 31 recorded events because they were two of the largest events with the most detailed (i.e., complete and fine in temporal resolution) precipitation data. The saturated hydraulic conductivity and Manning's roughness coefficient were calibrated to fit the observed discharge at the Shejiagou hydrologic station. After that, the rainsplash coefficient and the surface erodibility coefficient were calibrated to fit the observed sediment discharge. The Nash-Sutcliffe efficiency [50] was used as the criterion to evaluate the model performance:

$$MSE = \left[ \sum\_{i=1}^{n} \left( O\_i - \overline{O} \right)^2 - \sum\_{i=1}^{n} \left( S\_i - O\_i \right)^2 \right] / \sum\_{i=1}^{n} \left( O\_i - \overline{O} \right)^2 \tag{1}$$

where *Oi* was the observed value (water discharge or sediment discharge), *O* was the average value of the observed values, and *Si* was the simulated value and *n* was the number of samples. The Nash-Sutcliffe coefficients of the calibration were 0.763 and 0.738 for river flow and sediment discharge, and those of the validation were 0.781 and 0.732, respectively (Figure 4). These results showed good performance of the InHM with the optimized parameter set in the Shejiagou catchment. This parameter set (Tables 3 and 4) was then used in the scenario modelling runs.

**Figure 4.** Hydrographs and sedigraphs of the model calibration ((**a**) and (**b**), event of 14 July 1964) and validation ((**c**) and (**d**), event of 2 August 1964). Note that the blue bars represent the rainfall intensity.

The Manning's roughness coefficient values were calibrated to be 0.010 s m−1/<sup>3</sup> for the hillslopes and 0.008 s m−1/<sup>3</sup> for the gullies (including the main channel and sedimentary land) and 0.004 s m−1/<sup>3</sup> for the dam body in the catchment (Table 4). The calibrated Manning's roughness coefficient values were small but reasonable due to the low vegetation coverage and the large area of bare land and cultivated land in the Shejiagou catchment during the simulation period (i.e., 1964) [25,27]. Similar values have been measured on bare land and cultivated land through field experiments in the Loess Plateau in a few studies e.g., [51,52].

#### 2.4.5. Evaluation Criteria

To evaluate the influences on flood processes under different deployment strategies, the hydrograph and sedigraph at the catchment outlet were compared among the eight scenarios. Flood peak discharge (Qp), peak time (Tp), and flood volume (Qv) of each scenario were derived from the hydrographs. The runoff ratio (RR), which reflected the influences of check dams on the runoff-precipitation relationship in this study, was also calculated through dividing the runoff depth by the rainfall depth. To compare the efficiencies of sediment reduction under different deployment strategies, indicators such as Qs,p (sediment discharge peak), SY (sediment yield observed at the catchment outlet), and ES (the total mass of eroded sediment in the catchment) were derived from the modelling results. Besides, the sediment delivery ratio (SDR) was calculated as follows:

$$\text{SDR} = \frac{\text{SY}}{\text{ES}} \times 100\% \tag{2}$$

SDR has multi-fold environmental implications both in evaluating the soil erosion processes of a catchment and the effectiveness of soil and water conservation measures conducted in a catchment [53]. Variations of SDR of a specified catchment can be induced either by changes in SY or in ES. A comprehensively analysis of the changes in SDR, SY, and ES under different scenarios led to a more complete understanding of the changes in sediment transport processes caused by filled check dams. Moreover, the final distribution of eroded or deposited sediment in the catchment was also analyzed to furtherly understand the influences of filled check dams on sediment transport processes.

#### **3. Results**

#### *3.1. Variation Characteristics of Flood Processes under Di*ff*erent Deployment Strategies*

The flood processes were significantly changed by filled check dams under different deployment strategies (Table 5, Figure 5). Compared to scenario-0 (i.e., no check dams), the existence of check dams and their sedimentary lands reduced and delayed the flood peak. The efficiencies of flood-peak reduction and delay under the 7 deployment strategies were scenario-6 > scenario-4 > scenario-2 > scenario-7 > scenario-5 > scenario-3 > scenario-1. The hydrographs at the catchment outlet showed two different patterns: (1) the flood discharge rose and dropped sharply with a large peak discharge in scenario-0, scenario-1, scenario-3 and scenario-5; (2) the rising limb of hydrographs were delayed and more gentle, and the recession limb of hydrographs were longer in scenario-2, scenario-4, scenario-6, and scenario-7. The major difference between scenarios 2,4,6 and scenarios 1,3,5 was the existence of the large check dam (e.g., S2-2) on the middle reach of the main channel (Figure 2), which indicated that two large check dams constructed in series on the main channel could be more beneficial in flood peak attenuation than scenarios with only one large check dam. For scenarios with only one large check dam (i.e., scenarios 1,3,5), increasing the number of small check dams on the gullies led to significant reduction of flood peak and insignificant delay of peak time. For scenarios with two large check dams (e.g., scenarios 2,4,6), the delay of flood peak was longer when adding more small check dams in the catchment. Besides, proper site selection of large check dams on the main channel also made a difference in flood peak attenuation (e.g., scenario-6 versus scenario-7).


**Table 5.** Flood characteristics under different scenarios.

Different from unfilled check dam systems, which not only reduce flood peak but also reduce flood volume significantly [24], the reductions of flood volume by filled check dam systems were smaller (Table 5). For unfilled check dams, surface flow is intercepted by the dam body until the water depth is higher than the inlet of spillway, which forces a large amount of water to retain behind the check dam. For filled check dams, the elevation superiority of the spillway over the upstream channel no longer existed, and the flood volumes were mainly reduced by the promoted infiltration processes and the slowdown of flood velocity on the wider channel (i.e., sedimentary lands). Even so, the flood reduction rate can still be as much as 34.10% in scenario-6, in which the size combination was two large check dams and two small check dams. The efficiencies of flood-volume reduction under the 7 deployment strategies were scenario-6 > scenario-7 > scenario-4 > scenario-2 > scenario-5 > scenario-3 > scenario-1. Different from the flood peak, the flood-volume reduction rate of scenario-7 was larger than scenario-3 and scenario-4. A possible explanation for this result may be that there were more check dams and sedimentary lands in scenario-7 in the headwater areas, where more runoff was generated compared to other areas [22]. In scenario-7, the sedimentary lands of the three upstream check dams were connected and became a larger infiltration zone for the upstream runoff than those in scenario-2 and scenario-4. Under the same rainfall condition (i.e., 90 mm), changes of the runoff coefficient under different deployment strategies were consistent with those of the flood volume observed at the catchment outlet.

**Figure 5.** Hydrographs at the catchment outlet under different scenarios (Q: discharge).

#### *3.2. Variation Characteristics of Sediment Transport Processes under Di*ff*erent Deployment Strategies*

The changes of sediment transport processes at the catchment outlet were generally consistent with the variation characteristics of flood processes under different deployment strategies (see Figures 5 and 6). Along with the reduction and delay in flood peaks, the sediment discharge peaks (i.e., Qs,p) were also reduced and delayed. The efficiencies of Qs,p reduction under the 7 deployment strategies were scenario-6 > scenario-4 > scenario-2 > scenario-5 > scenario-7 > scenario-3 > scenario-1 (Table 6). This result highlighted the importance of adding a large check dam on the middle reach because there was a large check dam on the middle reach of the catchment in scenarios 6,4,2 and no large check dam on the same position in scenarios 5,7,3,1 (see Figure 2). Differently, the reduction rate of Qs,p in scenario-5 (i.e., 67.43%) was larger than that in scenario-7 (i.e., 59.02%), while the reduction rate of Qp in scenario-5 (i.e., 69.12%) was slightly smaller than that in scenario-7 (i.e., 71.45%). This was because that the gully controlled by check dam S5-2 contributed a large amount of eroded sediment into the main channel (see Figure 6 in [22]), and there was no check dam on the heavily eroded gully in scenario-7, indicating that targeted treatment to heavily eroded gullies could lead to a better performance of the check dam system in sediment reduction.

The reduction rates of sediment yield at the catchment outlet were basically higher than those of flood volume except for scenario-2 and scenario-7 (Tables 5 and 6), indicating that the deposition-promotion performance of sedimentary lands was normally better than their infiltration-promotion performance under most of the deployment strategies. The efficiencies of sediment-yield reduction under the 7 deployment strategies were scenario-6 > scenario-4 > scenario-5 > scenario-7 > scenario-2 > scenario-3 > scenario-1, implying that an increase in the number of check dams and targeted treatment to heavily eroded gullies could lead to a better performance in reducing sediment output of the catchment. Table 6 also showed that the total amount of eroded sediment in the catchment after a flood decreased due to the existence of filled check dams, and the reduction efficiencies were scenario-6 > scenario-4 > scenario-7 > scenario-2 > scenario-5 > scenario-3 > scenario-1. The reduction of total eroded sediment might be attributed to the lifting of erosion base of channels and hillslopes around the sedimentary lands, which led to decrease of runoff erosion power [24,54]. Large decreases of the sediment delivery ratio (SDR) were found in scenario-6 and scenario-4, which were consistent with the results of sediment yield at the catchment outlet. However, for scenarios with both small decreases of sediment yield and eroded sediment mass (e.g., scenario-1, scenario-2, scenario-7), their benefits for sediment reduction were hard to distinguish from SDR.

**Figure 6.** Sedigraphs at the catchment outlet under different scenarios (Qs: sediment discharge).


**Table 6.** Sediment transport characteristics under different scenarios.

#### *3.3. Erosion*/*Deposition Distribution in the Catchment under Di*ff*erent Deployment Strategies*

The final distributions of erosion and deposition areas under the eight scenarios were shown in Figure 7. When there was no check dam in the catchment (i.e., scenario-0), erosion processes occurred both on hillslopes and in channels, while deposition areas were sparsely scattered in gully-head areas (Figure 7a). However, the existence of check dams, even when they were filled, could transform a large area of channels and gullies into deposition sinks (Figure 7b–h). By comparing the deposition depths behind each check dam (i.e., in the sedimentary land), we found that more sediment was deposited on the upstream sedimentary lands (i.e., the tail-end areas of the sedimentary land) than on the areas that were close to the check dams (e.g., Figure 7b). This deposition pattern was different from that of an unfilled check dam. For an unfilled check dam, flood processes were stopped at the near-dam area, and consequently, deposition processes occurred there. Therefore, the thickness of the deposition layer was thicker at the near-dam area than at the upstream area [55]. For a filled check dam with a large and flat sedimentary land, which was considered in this study, flood processes were weakened as soon as the surface flow entered the tail-end areas of sedimentary lands due to the decrease of slope gradient [21], consequently promoting sediment to deposit on these areas.

**Figure 7.** The final distributions of erosion and deposition areas in scenario-0 (**a**); scenario-1 (**b**); scenario-2 (**c**); scenario-3 (**d**); scenario-4 (**e**); scenario-5 (**f**); scenario-6 (**g**); and scenario-7 (**h**). Note that the inset pictures in each figure show the deposition depths on the sedimentary lands, and the areas with deposition depth below zero (i.e., erosion) were blanked for the consideration of readability.

Perusal of the deposition distributions under different deployment strategies led to the following implications:


Deducing from those implications, in the situation where all check dams were filled, the best deployment strategies of check dams were those with large and connective sedimentary lands on the main channel and targeted treatments to heavily eroded gullies (e.g., scenario-6 in this study).

#### **4. Discussion**

#### *4.1. Potential Benefits of Filled Check Dams in Flood and Erosion Control*

Unlike check dams with large remaining storage capacity (i.e., low fill rate), filled check dams influence the flood processes mainly by the sedimentary land formed behind the check dams. Characterized as having a relatively high permeability [20] and a low slope gradient [46], the sedimentary land is able to attenuate flood velocity, and locally promote water infiltration and sediment deposition, consequently reducing and delaying the flood processes (Figure 5). Compared to the original channel, the sedimentary land is less susceptible to hydraulic erosion due to the attenuation of flood and the relatively large size of surface particles on the sedimentary land [56]. Therefore, the upstream expansion of sedimentary land will protect the soils of some channel areas and slope areas from being eroded, indirectly reducing the erodible area of the catchment. For example, the sedimentary land of S1-1 (i.e., the large check dam in scenario-1) protected approximate 950 m of the original channel reach and covered around 2.33 <sup>×</sup> <sup>10</sup><sup>4</sup> m2 of hillslopes around the channel, indirectly reducing by 900 t (2.18%) eroded sediment. This indirect benefit of a filled check dam, which has not been quantitatively analyzed by previous studies, was termed as the lifting of erosion base by [54].

The increase of sedimentary lands resulted in considerable reductions in flood peak discharge (Qp) and sediment peak discharge (Qs,p), and the reduction rates responded linearly to the area of sedimentary land (Figure 8), which implies the long-term potential benefits of the large amount of check dams in the Loess Plateau. Although the remaining storage capacity of them will be totally consumed in the future, the sedimentary lands behind check dams will still be effective in flood control. As the area of sedimentary land increases, the reduction rates of flood volume (Qv), sediment yield (SY), and eroded sediment (ES) were smaller than those of Qp and Qs,p (Figure 8), because these indicators were also sensitive to the size combination and site selection of check dams (Tables 5 and 6).

Check dams, especially small check dams, are filled quickly in the catchment with severe erosion problems (e.g., in the Loess Plateau). For example, 444 check dams (including 39 large check dams) were constructed during 1953–1977 in the 200 km2 Chabagou watershed and more than 60% of them were filled according to the check dam survey in 1978 (Figure 9). Several studies have reported that the flood processes were obviously attenuated and the runoff and sediment yield out of the Chabagou watershed were largely decreased due to check dam construction and hillslope treatments (e.g., terrace farming, afforestation) [57–59]. Based on the results of this study, it can be inferred that these filled check dams and their sedimentary lands also had a nonnegligible contribution to the reduction of runoff and sediment yield. Quantitatively modelling the cumulative effect of flood attenuation by the

large number of filled check dams will help to obtain a more comprehensive understanding of the potential benefits of check dams in the long-term.

**Figure 8.** Relationship of the total area of sedimentary lands and reduction rate of flood processes and sediment transport processes.

**Figure 9.** Fill rate of check dams in the Chabagou watershed. Note that when the check dam is filled, its fill rate equals to 1.

#### *4.2. Implications of Future Check Dam Deployment*

The quick-fill of check dams due to the large sediment yields during wet seasons and frequent dam breaks caused by flood overtopping often shorten the lifespan of check dams in the Loess Plateau [17]. For example, the two upstream check dams in the Shejiagou catchment (i.e., S7-3 and S7-4 in scenario-7) were totally filled after ten year's operation according to the check dam survey in 1978 [22]. There were technically two possible solutions to this problem: (1) dam maintenance such as regular sediment removal behind check dams can create new sediment storage capacities for

filled check dams; (2) construction of open check dams, which have a permanent outlet and can let more sediment be transported downstream [60]. However, the first solution seems to be unfeasible in most of the Loess Plateau catchments due to the large amount of check dams on the Loess Plateau and the high expense for large excavators to transport and work in the gullies. On the other hand, the sediment reduction rate of open check dams is much lower than the check dam studied in this study, which cannot meet the urgent needs of sediment reduction for the heavily eroded Loess Plateau catchment. Besides, the construction materials of the dam body for open check dams are usually stone or concrete, which increase the construction expense in the hilly-gullied Loess Plateau catchment. Instead, check dams in the Loess Plateau catchment are usually gravity dams constructed by the local earth using materials from the nearby hillslopes, which largely reduces the construction expenses (see the Supplementary Materials).

Therefore, a scientific deployment strategy of check dams is necessary in order to obtain the largest benefits of flood attenuation and sediment interception when the storage capacities of check dams are consumed in the future. According to Table 7, scenarios with two large check dams (e.g., scenarios 2,4,6,7) performed better than scenarios with only one check dam (e.g., scenarios 1,3,5) in the reduction of Qp and Qv. Except for scenario-5, scenarios with two large check dams (i.e., scenarios 2,4,6,7) also performed better than scenarios with only one large check dam (i.e., scenarios 1,3) in the reduction of Qs,p, SY, and ES. Apart from the number of large check dams in the system, the locations and number of small check dams also made a difference in the reduction performance. For example, scenario-5 performed better than scenario-7 in the reduction of Qs,p and SY because there was targeted treatment (i.e., check dam S5-2) to a heavily eroded gully in scenario-5. In the aspect of flood control, constructing large check dams on the main channel and making the sedimentary lands of the large check dams as connective as possible will obtain the best performance in the reduction and delay of flood peak (Table 7). Similar results were reported in [20], which showed that the combination of more than one large check dams on the main channel led to a considerable increase of groundwater recharge. In the aspect of sediment control, targeted treatment to specified gullies combined with large and connective sediment sinks in the main channel will largely promote sediment deposition in the channel and reduce soil erosion of hillslopes and channels (Table 7). Compared the current deployment strategy (i.e., scenario-7), the deployment strategy in the hypothetical scenario-6 which was more in line with the description above indeed performed better both in flood and sediment control.


**Table 7.** Performance ranking in flood and sediment control among different scenarios.

<sup>a</sup> Scenarios with different check dam systems (i.e., scenarios 1–7) were compared; <sup>b</sup> **(1)** represents the best reduction performance among all scenarios, **(7)** represents the worst reduction performance among all scenarios. For example, the best performance in the reduction of Qp was scenario-6.

#### *4.3. Limitations and Future Work*

There were two simplifications in the parameterization of sedimentary land to simplify the problem in this study: (1) the elevations of surface nodes on the same sedimentary land were the same (i.e., the slope gradient of sedimentary land was zero) and (2) no croplands were planted on the sedimentary land. In reality, the slope of a sedimentary land in a check dam system is normally less than 2%. The first simplification might lead to a slight overestimation of the flood attenuation effects of sedimentary lands. However, most sedimentary lands in the Loess Plateau are planted with crops (e.g., corn), which increase surface roughness and local infiltration. This led to a small underestimation of

the decrease in flow velocity and the sediment deposition. The two simplifications only had small influences on the flood processes, and their influences might offset each other [61].

While the implications on deployment strategies of check dams were specified to Shejiagou catchment, we believe that similar findings may be obtained for larger catchments with more complex check dam systems (e.g., other sub-catchments of Chabagou watershed in Figure 9). For this reason, future research will focus on larger sites with more complicated stream networks and check dam systems, to validate the implications and to provide a more comprehensive guidelines for the optimization of check dam deployment strategies. Besides, the influence of sedimentary land on catchment hydrologic responses and sediment transport processes in a larger time scale instead of event-scale need to be studied in order to aid the prediction of future changes in runoff and sediment yield in the Loess Plateau.

#### **5. Conclusions**

Employing the InHM, this study investigated the influences of filled check dams with different deployment strategies on the flood and sediment transport processes in a Loess Plateau catchment, which is a 4.26 km2 erosion-prone area and currently has four check dams. Six hypothetical scenarios representing different deployment strategies were compared with no-dam scenario and the scenario that is currently applied in this catchment. Event-scale rainfall-runoff simulations were conducted with a 90 mm hr−<sup>1</sup> rainfall intensity for the eight scenarios. Indicators such as flood peak discharge (Qp), peak time (Tp), flood volume (Qv), sediment peak discharge (Qs,p), sediment yield (SY), and eroded sediment (ES) were used to compare the performances of different deployment strategies in flood and sediment control.

Filled check dams were still able to reduce flood peak by 31% to 93% under different deployment strategies. Considerable delays of peak time and decreases were also found in scenarios 2,4,6, which were characterized as having larger and more connective sedimentary lands on the main channel. Reduction rates of sediment yield and eroded sediment ranged from 4% to 52% and 2% to 16%, respectively, indicating that a proper distribution of check dams was beneficial to promote sediment deposition in the channel and reduce soil erosion of hillslopes and channels.

Apart from showing the potential benefits of filled check dams, this modelling exercise provided the following insights on the deployment strategy of check dams, which led to a better performance in flood and sediment control when check dams are filled:


Compared to the current deployment strategy (i.e., scenario-7), a more scientific deployment strategy for the selected catchment, which is more in line with the insights above, was proposed (i.e., scenario-6).

The simulations reported herein took a first step to explore the potential benefits of filled check dams and demonstrated that check dams could have a long-term influence on the flood and sediment transport processes via forming sedimentary lands that changes the catchment landforms permanently. *Water* **2020**, *12*, 1319

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/5/1319/s1, Figure S1: (a) Sediment erosion (negative value) and deposition (positive value) rates for Shejiagou catchment; (b) Sediment inflow to stream segments; (c) Cumulative flow without check dams (Debasish et al., 2018), Figure S2. Open check dams. Beam dam (left), silt dam (right). Derived from Armanini and Larcher (2001), Figure S3. Characteristic components of a sediment trap with an open check dam: (a) inlet structure: solid body dam; (b) scour protection; (c) basin; (d) lateral dikes; (e) maintenance access; (f) open check dam; (g) counter dam. Derived from Piton and Recking (2015), Figure S4. Examples of check dams in the Loess Plateau. (a) check dam LDG-1 in Liudaogou catchment; (b) check dam MDZ-1 in Wangmaogou catchment; (c) check dam NYG-2 in Nianyangou catchment; (d) check dams system in Shejiagou catchment.

**Author Contributions:** Conceptualization, H.T.; methodology, H.T.; software, Q.R.; validation, H.T.; formal analysis, H.T.; investigation, H.P.; writing—original draft preparation, H.T.; writing—review and editing, Q.R.; supervision, Q.R.; funding acquisition, Q.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by the China Major Science and Technology Program for Water Pollution Control and Treatment (2017ZX07206-001) and the National Natural Science Foundation of PR China (Grant no. 51379184 and no. 51328901).

**Acknowledgments:** We thank the Zizhou Experimental Station of YRCC and the National Earth System Science Data Center, National Science & Technology Infrastructure of China (http://www.geodata.cn) for providing essential observation data. We also thank the two anonymous reviewers for providing constructive suggestions. Special thanks to Amplin for helping improve the language of the manuscript.

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

#### **Appendix A. The Integrated Hydrologic Model (InHM)**

#### *Appendix A.1. Hydrologic-Response Module*

The 3D subsurface flow in variably saturated porous medium is estimated in InHM by:

$$\nabla \cdot f^a \stackrel{\rightharpoonup}{q} \pm q^b \pm q\_{ps}^c = f^v \frac{\partial \phi S\_w}{\partial t} \tag{A1}$$

$$
\overrightarrow{q} = -k\_{rw} \frac{\rho\_w \mathcal{g}}{\mu\_w} \overrightarrow{k} \nabla(\psi\_p + z) \tag{A2}
$$

where *<sup>f</sup> <sup>a</sup>*[-] is the area fraction related to each continuum, *q* [*LT*<sup>−</sup>1] is the Darcy flux, *qb*[*T*<sup>−</sup>1] is a specified rate source/sink, *qe ps*[*T*<sup>−</sup>1] (equal to <sup>−</sup>*qe sp*) is the rate of water exchange between the porous medium and surface continua, *f <sup>v</sup>*[-] is the volume fraction associated with each continuum, φ[*L*3*L*−3] is porosity, *Sw*[*L*3*L*<sup>−</sup>3] is the water saturation, *t*[*T*] is time, *krw*[-] is the relative permeability, ρ*w*[*ML*<sup>−</sup>3] is the density of water, *<sup>g</sup>*[*LT*−2] is the gravitational acceleration, <sup>μ</sup>*<sup>w</sup>* [*ML*−1*T*−1] is the dynamic viscosity of water, 

*k* [*L*2] is the intrinsic permeability vector, ψ*p*[*L*] is the pressure head, and *z*[*L*] is the elevation head.

The diffusion wave approximation to the depth-integrated shallow wave equations is employed in InHM to estimate the 2D transient surface flow on the land surface (both overland and open channel), with the conservation of water on the land surface expressed by:

$$\nabla \cdot \psi\_s^{\text{multi}} \stackrel{\text{combile}}{q}\_s \pm a\_s q^b \pm a\_s q\_{sp}^\varepsilon = \frac{\partial (S\_{w\_s} H\_l + \psi\_{im})}{\partial t} \tag{A.3}$$

$$\overrightarrow{q}\_s = -\frac{\left(\psi\_s^{\text{mobile}}\right)^{2/3}}{\overrightarrow{n}\,\Phi^{1/2}}\nabla(\psi\_s + z)\tag{A4}$$

where <sup>ψ</sup>*s*[*L*] is the surface water depth, *<sup>q</sup> <sup>s</sup>*[*LT*−1] is the surface water velocity calculated by a two-dimensional form of the Manning's equation, *as*[*L*] is the surface coupling length scale, *qb*[*T*<sup>−</sup>1] is the source/sink rate (i.e., rainfall/evaporation), *qe sp*[*T*<sup>−</sup>1] is the rate of water exchange between the surface continua and porous medium, *Ht*[*L*] is the average height of non-discretized surface microtopography, *n*[*TL*−1/3] is the Manning's surface roughness tensor, and Φ[-] is the friction (or energy) slope; ψ*mobile <sup>s</sup>* [*L*] and ψ*im*[*L*] refer to the water depth exceeding and held in depression storage, respectively.

The calculated daily potential evapotranspiration was then incorporated into InHM as actual evapotranspiration (ET) using a set of sink functions:

$$Q\_b^E = \alpha(\psi) q\_{\text{max}}^E A\_b \tag{A5}$$

where *QE <sup>b</sup>* [*L*3*T*<sup>−</sup>1] represents the volumetric evapotranspiration rate, <sup>α</sup>(ψ)[-] is a response function of the saturation of the porous medium and the degree of ponding at the land surface, *qE* max[*LT*<sup>−</sup>1] is the potential evapotranspiration rate per unit area estimated by the revised Penman–Monteith method, ψ[*L*] is the pressure head of the subsurface nodes or water depth of the surface nodes, and *Ab*[*L*2] is the area associated with the surface water equation

The first-order coupling between the surface and subsurface continua, driven by pressure head gradients, occurs via a thin soil layer of thickness, *as* in Equation (A3). The control volume finite-element method is employed to discretize the equations in space. Each coupled system of nonlinear equations is solved implicitly using Newton iteration.

#### *Appendix A.2. Sediment-Transport Module*

Depth-integrated multiple-species sediment transport, restricted to the surface continuum, is calculated for each sediment species by:

$$\frac{\partial \mathbb{C}\_{\rm{scd}}}{\partial t} = -\nabla \cdot \overrightarrow{q}\_{s} \mathbb{C}\_{\rm{scd}} + \frac{1}{V\_{\rm{uv}}} (\mathbf{e}\_{\rm{scd}} + \sum\_{j=1}^{\rm{BC}} q\_{s\_{j}}^{b} \mathbb{C}\_{\rm{scd}\_{j}}^{\*}) \tag{A6}$$

$$
\varepsilon\_{\rm sel} = \varepsilon\_{\rm 6} + \varepsilon\_{\rm h} \tag{A7}
$$

where *Csed*[*L*3*L*<sup>−</sup>3] is volumetric sediment concentration, *<sup>q</sup> <sup>s</sup>*[*LT*<sup>−</sup>1] is the depth-averaged surface water velocity, *Vw*[L3] is the volume of water at the node, *esed*[*L*3*T*−1] is the volumetric rate of soil erosion and/or deposition, *qb sj* [*L*3*T*−1] is the rate of water added/removed via the *<sup>j</sup>*th boundary condition, *C*∗ *sedj* [*L*3*L*<sup>−</sup>3] is the sediment concentration of the water added/removed via the *j*th boundary condition, *BC* is the total number of boundary conditions. The net erosion rate is the sum of the rainsplash erosion rate *es*[*L*3*T*<sup>−</sup>1] and the hydraulic erosion rate *eh*[*L*3*T*<sup>−</sup>1]. The rainsplash erosion rate of each species in Equation (A7) is calculated as:

$$\mathfrak{e}\_{s\_i} = \begin{cases} \sigma\_i c\_f F(\psi\_s) (\cos(\theta) \cdot r)^b A\_{3D\_\prime} & q > 0 \\ 0 & \text{,} \qquad q < 0 \end{cases} \tag{A8}$$

$$F(\psi\_{\mathfrak{s}}) = \exp(-c\_d \psi\_{\mathfrak{s}}) \tag{A9}$$

where σ*i*[-] is the source fraction of species *i*, *c <sup>f</sup>*[(*TL*−1) *<sup>b</sup>*−1] is the rainsplash coefficient, *b* [-] is the rainfall intensity exponent, θ[-] is the angle of the element from horizontal, *r*[*LT*−1] is the rainfall intensity, *A*3*D*[*L*2] is the three-dimensional area associated with the node, and *q*[*LT*−1] is the sum of rainfall intensity and infiltration intensity; *F*(ψ*s*)[-] is a damping function, related to surface water depth ψ*s*[*L*] and the surface water damping-effectiveness coefficient *cd*[*L*<sup>−</sup>1], to represent the reduction in splash erosion with increasing surface water depth. The hydraulic erosion rate in Equation (A7) is estimated as:

$$
\sigma\_{h\_i} = \alpha\_{\text{scd}\_i} (\mathbb{C}\_{\text{scd}\_{\text{mu} \times\_i}} \sigma\_i - \mathbb{C}\_{\text{scd}\_i}) \tag{A10}
$$

$$\mathcal{C}\_{\text{sol}\_{\text{max}\_i}} = 0.05 \frac{\overrightarrow{q}\_s q\_\*^3}{g^2 d\_{\text{sol}\_i} \psi\_s \left(\nu\_{\text{scd}\_i} - 1\right)^2} \tag{A11}$$

$$\alpha\_{\rm scd} = A \begin{cases} 2v\_{\rm scd\_i} \xi & , \qquad \mathsf{C}\_{\rm scd} > \mathsf{C}\_{\rm scd\_{\rm max\_i}} \ (erosion) \\ \varphi \, \eta\_{\, s} \psi\_{\rm s} \chi\_{i} \, \, \vert \, & \qquad \mathsf{C}\_{\rm scd} < \mathsf{C}\_{\rm scd\_{\rm max\_i}} \ (depposition) \end{cases} \tag{A12}$$

where α*sedi* [*L*3*T*−1] is the hydraulic erosion transfer coefficient for species *<sup>i</sup>*, *Csed*max*<sup>i</sup>* [*L*3*L*−3] is the concentration at equilibrium transport capacity for species *<sup>i</sup>*, *<sup>q</sup>*∗[*LT*<sup>−</sup>1] is the local shear velocity, *dsedi* [*L*] and γ*sedi* [-] are the particle diameter and specific gravity, respectively; *A*[*L*2] is the area associated with the node in Equation (A12), *vsedi* [*LT*−1] is the particle settling velocity, ξ[-] is a coefficient related to turbulence in the surface water due to raindrop impact, ϕ[*L*<sup>−</sup>1] is an erodibility coefficient related to surface properties and texture, and χ*i*[-] is the particle erodibility factor (ranging from zero to one).

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Internal Stability Evaluation of Soils**

#### **Qingfeng Feng 1, Hao-Che Ho 2, Teng Man 3, Jiaming Wen 1, Yuxin Jie <sup>1</sup> and Xudong Fu 1,\***


Received: 18 May 2019; Accepted: 9 July 2019; Published: 12 July 2019

**Abstract:** Suffusion constitutes a major threat to the foundation of a dam, and the likelihood of suffusion is always determined by the internal stability of soils. It has been verified that internal stability is closely related to the grain size distribution (GSD) of soils. In this study, a numerical model is developed to simulate the suffusion process. The model takes the combined effects of GSD and porosity (*n*) into account, as well as Wilcock and Crowe's theory, which is also adopted to quantify the inception and transport of soils. This proposed model is validated with the experimental data and shows satisfactory performance in simulating the process of suffusion. By analyzing the simulation results of the model, the mechanism is disclosed on how soils with specific GSD behaving internally unstable. Moreover, the internal stability of soils can be evaluated through the model. Results show that it is able to distinguish the internal stability of 30 runs out of 36, indicating a 83.33% of accuracy, which is higher than the traditional GSD-based approaches.

**Keywords:** suffusion; internal stability; grain size distribution (GSD)

#### **1. Introduction**

Suffusion is the mass movement of fine particles through the pore space of a coarser matrix driven by seepage forces [1]. It is one of the main mechanisms initiating internal erosion within levees, earth dams and foundations [2–4] as well as watershed hillslopes [5,6].

A lot of experimental [7–14] and numerical (CFD) studies [15,16] were launched and expected to reveal the underlying mechanisms corresponding to suffusion. Most of these studies were macroscopic and under the important assumption that the movement of soil is a continuous process. This assumption has never been verified through or connected directly to a particle-scale study. However, some of following studies in granular system start to show their potential to support this assumption. For instance, as DEM models [17–19] and CFD+DEM models [20] are developed, the suspension, collision or spin of particles have been described in detail, but still seem to be irrelevant to the internal structure fracture (based on regime transition research [21]) or massive movement (based on rheology research [22–24]). Since then, macroscopic studies on suffusion are still mainstream, especially when focusing on engineering practices.

Suffusion can take on many forms, but it is particularly insidious when it occurs under a relatively low hydraulic gradient. Hydraulic engineers determine the internal stability by examining the critical gradient of a soil which directly affects the likelihood of diffusion. Generally, the soil is internally stable if *i*<sup>c</sup> < 0.7 [7].

The internal stability of soils is closely related to its grain size distribution (GSD). As pointed out by Wan and Fell [1], internally unstable soils are usually broadly graded soils with particles from silt (or clay) to gravel size, whose GSD curves are concave upward (or known as gap-graded soils). GSD-based assessment approaches, based on experimental evidence, have been used extensively in practice [8–14]. Most of them distinguish coarse and fine particles within the soil and incorporate characteristic sizes into an assessment index or indices.

Early studies accounted for two particles sizes with an additional division size. For example, Kézdi [25] divided soils into coarse and fine components. The maximum value of *D*15/*d*85, which depends on the division size, is used as the index of internal stability. *D*<sup>15</sup> (*d*85) is the size for which 85% (15%) by weight is finer in the coarse (fine) components, respectively. The soil is evaluated as unstable if *D*15/*d*<sup>85</sup> > 4. Similarly, Kenney and Lau [26] introduced a "shape curve" to identify an internally unstable soil. They defined *F* as the percentage of particles that are finer than an adjustable size *d* and *H* + *F* as the percentage of particles finer than the size of 4*d*. The material is assessed as internally unstable when the minimum value of *H*/*F* is less than 1.3 originally, [26] and a modified lower value of 1.0 later [27]. Different from the usage of an adjustable particle size, Burenkova [28] considered three characteristic sizes and evaluated the soils unstable when *d*90/*d*<sup>60</sup> is greater than 1.86log(*d*90/*d*15) + 1 or less than 0.76log(*d*90/*d*15), where *d*<sup>x</sup> (*x* = 15, 60 or 90) hereafter is the diameter for which *x*% by weight of the soil particles are finer.

In recent decades, the idea of a more precise representation of the GSD curve has been further developed. Wan and Fell [1] extracted four sizes from the GSD curve. Two are for coarse particles and the other two for the fine counterpart of the soils. They defined two variables, *s*<sup>1</sup> = 15/log(*d*20/*d*05) and *s*<sup>2</sup> = 30/log(*d*90/*d*60), for each soil sample. The well-known assessment criterion is that a soil sample is internally unstable when *s*<sup>1</sup> < 22 and *s*<sup>2</sup> > 80. That is, an internally unstable soil is featured as *d*20/*d*<sup>05</sup> > 4.8 and *d*90/*d*<sup>60</sup> < 2.4. Also, its accuracy is higher than any other GSD-based approaches mentioned above.

The accuracy of the GSD-based approaches has been improving since more and more GSD information has been added in. In terms of seepage flow dynamics, however, these approaches are still empirical and with ambiguous meaning. Early research suggests that the gap-graded soils are apt to being internally unstable since fine particles easily pass through, but do not fill in the voids between the coarse particles if subjected to seepage [9]. The filtration of detached fine particles has therefore been well accepted as the mechanism responsible for the internal stability. Overall, suffusion appears as the combination of detachment, transport, and filtration of fine particles. They may depend on the excessive shear stress related to particle inception [15,29], sediment transport capacity if definable, and the size of the voids, respectively. In addition to the filtration effect, the other mechanisms may dominate in some cases, but have received less attention so far. Such incomplete understanding hampers the explanation of other important factors such as porosity [30,31] and loading history [32].

This paper presents a numerical simulation of fine particle dynamics subjected to seepage flow, within an extended numerical framework of Fujisawa's study [15]. We take into account the combined effects of GSD and porosity. Wilcock and Crowe's theory [33] is adopted to quantify the inception and transport of each size group of soil particles. We reproduce the phenomena observed in seepage tests and explain why the soils with specific GSD behave more internally unstable, and we predict the internal stability of soil through the model.

#### **2. Numerical Model**

Numerical models have been successfully applied to analyze seepage flow through embankments, dikes, and soil foundations. With a reasonable model for particle movement within soil matrix, numerical models may be able to evaluate internal stability of soils with a certain GSD. Our new model is an alternative approach for exploring internal stability.

#### *2.1. Governing Equations*

In terms of representing the dynamics of fine particles subjected to seepage forces, numerical simulation at the microscopic scale provides a potential approach. In 2010, Fujisawa et al. described the process of suffusion (or the process of porosity *n* change with time *t*) within a unit volume of seepage field by using ∂*n*/∂*t* = *EA* [15]; The involved parameters *E* and *A* are the erosion rate defined as the volume of eroded particles per erodible area per time, and the erodible area per unit volume.

Taking the GSD into consideration, as shown in Figure 1, we assume that the particles from the soil skeleton would be eroded into the pore flow and increase its concentration, considering the two-dimensional seepage flow in the directions 1 and 2. During the time period of *dt*, the volume of eroded soil particles is -*EkAkdVdt*, where *Ek* denotes the erosion rate of *k*th size group section defined as the volume of eroded particles per erodible area per time, *Ak* denotes the erodible area per unit volume of *k*th size group, *dV* = *dx*1*dx*2*dL*, and *dx*1, *dx*2, and *dL* are the dimensions of the unit volume, respectively.

**Figure 1.** Unit volume of seepage field as the modeling framework.

The governing equations of mass and momentum conservation are as follows:

$$\frac{\partial \mathfrak{n}}{\partial \mathfrak{t}} = \sum E\_k A\_k \tag{1}$$

$$\frac{\partial n\mathbf{C}}{\partial t} + \frac{\partial v\_{\text{i}}\mathbf{C}}{\partial \mathbf{x}\_{\text{i}}} = \sum E\_{k}A\_{k} \tag{2}$$

$$m\frac{\partial \mathbf{S}\_r}{\partial t} + \frac{\partial v\_i}{\partial \mathbf{x}\_i} = 0 \tag{3}$$

$$
v\_i = -k\_s \frac{\partial \mathbf{h}}{\partial \mathbf{x}\_i} \tag{4}$$

where *n* denotes the porosity of soil mass, *C* denotes the solid concentration of seepage flow, *vi* denotes the seepage velocity in the two directions (*i* = 1, 2), *Sr* denotes the saturation of soils, *h* denotes the local water head, *ks* denotes the hydraulic conductivity, *xi* denotes the distance along the *i*th direction (*i* = 1, 2), and *t* denotes the time passing by. Details about the closure equations will be introduced in Section 2.2.

Item -*EkAk* differs from *EA* in the study of Fujisawa et al. [15] since the dependence of *Ek* on the size of *k*th group is taken into account. If the dependence is neglected, i.e., *Ek* = *E* holds for each size group, the formulation by Fujisawa et al. [15] will be recovered.

#### *2.2. Closures*

Closures of items erosion rate *Ek*, erodible area *Ak*. and hydraulic conductivity *ks* are listed below:

#### 2.2.1. The Erosion Rate *Ek*

For suffusion, a linear relationship has been used between the erosion rate *E* and the excessive shear stress τ acting on eroded particles. However, this assumption does not distinguish the difference in the inception of particles with different sizes. Here, we adopt Wilcock and Crowe 's theory [33] for mixed sediments which found wide applications in geomorphic science and river engineering. That is, the erosion rate *Ek* is evaluated taking into account GSD, i.e.,

$$E\_k = \begin{cases} \alpha(\tau - \tau\_{ck}) & \text{for } \tau \ge \tau\_{ck} \\ 0 & \text{for } \tau < \tau\_{ck} \end{cases} \tag{5}$$

where τ is shear stress exerted on the soil particles, *Ek* and τ*ck* are the erosion rate and critical shear stress for the *k*th size group, respectively, provided that the soil particles can be divided into several size groups (group number ≥*k*) and the particles within the same size group have the same critical shear stress.

The selection of α is based on calibration. More details will be introduced in Section 3.3.

The shear stress τ that is exerted on soil particles, is generated by the seepage flow with the soil fabric. In 1980, Hillel [34] assumed that the seepage field within a unit volume (see Figure 2a) can be treated as equivalent tubes in a length of *L* and radius of *R* = *D*/2 (see Figure 2b).

**Figure 2.** (**a**) Seepage field and (**b**) its equivalent tubes within a unit volume.

As shown in Figure 2, the equivalence guarantees the same discharge due to seepage flow under the same hydraulic gradient as well as the same porosity in average. The shear stress τ is then defined at the wall of tube. Using the Poiseuille law, it can be obtained by the following:

$$
\pi = I \sqrt{2\rho g \frac{k\_s \mu}{n}} \tag{6}
$$

where *I* = *h*/*L* denotes the local hydraulic gradient; *g* denotes the acceleration of gravity; ρ = *C*ρ*<sup>s</sup>* + (1 − *C*)ρ*<sup>w</sup>* is the density of seepage flow; ρ*<sup>s</sup>* denotes the density of soil particles; ρ*<sup>w</sup>* denotes the density of pure water; μ = μ*w*(1 + 2.5*C*) is the dynamic viscosity of seepage flow with suspended sediments, which is corrected by the suspension concentration; and μ*<sup>w</sup>* denotes the viscosity of pure water. It is noted that Equation (6) quantifies an equivalent shear stress within the unit volume, which balances the pressure drop driving the seepage flow.

The equivalent shear stress τ in the model is a macroscopic measure of the microscopic seepage shear stress exerted on the soil particles in the unit volume, and the potential difference between these two items is compensated through the calibration of erodible coefficient α in Equation (5). Meanwhile, the erosion of the fine particles is assumed to be equivalent to the process of sediment transportation, such an idea is also introduced when calculating the critical shear stress τ*c*.

In previous studies [15,35], the critical shear stress τ*<sup>c</sup>* was assumed to be a constant. However, since the GSD cannot be ignored, coarse particles will protect the fine ones from being eroded (known as the "Hiding Function"). Thus, the critical shear stress varies when particle size differs (see Figure 3). By integrating Wilcock and Crowe's theory [33] for the inception and transport of non-uniform sediments, for the *k*th group containing the soil particles of size between *dk* and *dk*+1, the critical shear stress τ*ck* is evaluated based on Equations (7) and (8):

$$
\pi\_{ck} = \pi\_{csm} \left(\frac{d\_k}{d\_{sm}}\right)^b \tag{7}
$$

$$\frac{\tau\_{\rm csm}}{(\rho\_{\rm s} - \rho\_{\rm w})gd\_{\rm sm}} = A + B\exp(-20F\_{\rm s})\tag{8}$$

where *dsm* denotes the medium size of soil particles; τ*csm* is the corresponding critical shear stress; *b* denotes an exponent that may differ from the erosion in open channel flows and is calibrated in the simulation; *Fs* denotes the content ratio of fine particles (diameter *d* < 2 mm); and *A* and *B* are the two empirical parameters, and can also be calibrated during the simulation.

**Figure 3.** Effect of GSD in suffusion process.

Parameter *A* in Equation (8) is usually related to the cohesion of the soil. When there are cohesive particles in the soil, the value of *A* will be relatively higher [33]; parameter *B* in Equation (8) is usually related to fluid properties [33]. Due to the nature differences between the Poiseuille flow in our model and the open channel flow in the sediment transport process in Wilcock and Crowe's study [33], the value of *B* may also be different. In our simulation, these two parameters are obtained through calibration. More details will be introduced in Section 3.3.

#### 2.2.2. The Erodible Area Ak

Within a volume of *dV*, the *Ak* for the *k*th size group can be quantified with the total projective area of the given size *dk*, that is, the product of the projective area of a particle, π*dk <sup>2</sup>*, and the number of such particles per unit volume, 6(1 − *n*)(*Pk*<sup>+</sup><sup>1</sup> − *Pk*)/π*dk <sup>3</sup>*. Therefore, the *Ak* reads,

$$A\_k = 6(1 - n)(P\_{k+1} - P\_k) / d\_k \tag{9}$$

where *Pk*<sup>+</sup><sup>1</sup> and *Pk* denote the percentage passing by weight of soil particles for the size finer than *dk*<sup>+</sup><sup>1</sup> and *dk* respectively.

#### 2.2.3. The Hydraulic Conductivity *ks*

The hydraulic conductivity *ks* is related to the saturation, porosity, and GSD of the soil. A relationship for the hydraulic conductivity of the saturated soil is proposed, and it reads

$$k\_s = \Gamma d\_{10}^{1.565} \frac{n^{2.3475}}{\left(1 - n\right)^{1.565}} \tag{10}$$

where *d10* denotes the diameter of 10% for which soil particles by weight are finer, and Γ is a constant (calibrated in Section 3.3). Such a relationship is based on Chapuis' study [10], as well as the selection of exponents of each item.

#### 2.2.4. Update of *d*<sup>10</sup> in Equation (10)

In the case of suffusion driven by a seepage flow, the moving out of eroded particles will change the local porosity and GSD. As a result, the hydraulic conductivity *ks* as well as the critical shear stress τ*ck* has a feature of dynamic adjustment during the erosion process. In order to capture the dynamic features, it is necessary to formulate the temporal change of *Pk* corresponding to the particle size *dk*, and then to specify *d*<sup>10</sup> in Equation (11). For this purpose, one can quantify the *Pk*(*t* + *dt*) after a time period *dt* as

$$P\_k(t+dt) = \frac{P\_k(t)(1-n) - \sum\_{i=1}^k E\_i A\_i dt}{(1-n) - EA dt} \tag{11}$$

A rearrangement of Equation (11) without considering the second-order term results in

$$\frac{dP\_k(t)}{dt} = \frac{P\_k(t)EA - \sum\_{i=1}^k E\_i A\_i}{1 - n} \tag{12}$$

Accordingly, the size of *d*<sup>10</sup> corresponding to 10% for which soil particles by weight are fine after the time period *dt* can be quantified approximately as

$$d\_{10} = (d\_x + d\_{x+1})/2 \tag{13}$$

where *dx* and *dx*+<sup>1</sup> can be obtained through exhaustive method among all soil size groups, under the condition that item (*Px* − 10%) + (*Px*+<sup>1</sup> − 10%) reaches the smallest absolute value.

#### *2.3. Numerical Method and Solution Procedure*

The numerical model is solved by using a central scheme of the finite volume method (FVM) built on the OpenFOAM platform (version 2.2.2, https://openfoam.org/release/2-2-2/). More details when assigning basic variables (seepage velocity *vi* and porosity *n*) and governing equations are shown in Appendix A.

The flow chart indicating the solution procedure is illustrated in Figure 4.

**Figure 4.** Flow chart of simulation program (Eq(s): Equation(s)).

Variables in closures are updated in every cell of simulation domain, and calculations above are repeated at following time steps.

#### **3. Model Verification**

In this section, the proposed model is examined by previous experiments (studies by Skempton and Brogan in 1994 [9] and Wan and Fell in 2008 [1]). The experimental works have detailed records of hydraulic gradient (or eroded soil mass) at different stage of seepage flow and are therefore used for model calibration and validation.

Furthermore, in the next section, more experiments (studies by Lau in 1984 [36]; Lafleur et al. in 1989 [8]; and Chapuis in 1996 [10]) are also included when examining the model's capability in predicting internal stability.

#### *3.1. Overview of Experiments*

#### 3.1.1. Devices

The designs in experimental devices between these studies were similar: their experiments were conducted using acrylic cylinders filled up with fully mixed soils (see Figure 5). The upper end of the cylinder was the outlet of the seepage field, which allowed free outflow (also known as being downstream of the seepage field), and the discharge of the seepage flow was measured when running the experiment. The lower end was the inlet of the seepage field (also known as being upstream of seepage field), which connected the filter layer and the water pipe. The filter layer was designed to make the inflow distribution more stable, and the water pipe was used to connect the soil device with the water tank.

**Figure 5.** Designs of experimental devices.

#### 3.1.2. Observed Stages of Suffusion Process

With the increase of hydraulic gradient, Skempton and Brogan [9] observed the seepage flow and identified three stages of suffusion: (1) the "dancing-like movement of particles" stage, during which initial signal demonstrating that the suffusion had begun; (2) the "slight but general movement of particles" stage, during which suffusion kept developing and formed the shape of the "pipe" started stretching towards the upstream, and (3) the "strong general piping of particles" stage, where the soil particles were strongly moving out of the soil fabric.

Wan and Fell [1] gave similar characterizations. They identified the flow at the outlet as "clear", "cloudy", and "very cloudy" by manually observing and recording, because the degree of turbidity indicated the degree of suffusion, just like the judgment shown in Skempton and Brogan [9]. The "clear" flow indicated no suffusion and other conditions meant erosion at different stages. The soil samples were evaluated as internally unstable when the stage of "slight but general movement of particles" or "cloudy" flow was reached below the hydraulic gradient of a manually fixed value which was set as the critical one (e.g., 0.7 mentioned in Section 1).

#### 3.1.3. Soil Samples

Skempton and Brogan [9] used three different samples of non-cohesive soils and reported three runs of data, namely Runs A, B, C, respectively. In comparison, Wan and Fell [1] presented 18 runs of experimental data with 14 soil samples (upward flow tests only), their soils contained a portion of cohesive particles (diameter *d* < 0.02 mm). The maximum soil diameter was *d* = 30 mm. Both gap-graded soils and continuous-graded soils were used in the experiments. (see Figure 6).

The definition of gap-graded soils varies in different studies. However, among them, most gap-graded soils have irregular shapes. For instance, a horizontal part of the curve could be graded from 0.02 mm to 2 mm. Otherwise the soils are defined as continuous-graded.

For the runs that lacked complete curves (e.g., Runs 2 and 7), we assume that the diameter of particles in the missing part is the same size as the finest soil particles shown in curves, so that we can perform simulations for these runs.

**Figure 6.** GSD curves of soils in (**a**) internal stable runs; (**b**) internal unstable runs in the studies of Skempton and Brogan [9] (dashed lines) and Wan and Fell [1] (solid lines).

#### *3.2. Simulation Domain and Boundary Settings*

The simulation domain is proposed considering the circular symmetry of the seepage field, and assuming a two-dimensional flow at the vertical cross-section plane. As shown in Figure 7a, the domain covers a rectangle cross-sectional area. The size of the area is corresponding to the actual size of experimental device. Therefore, for the runs in Skempton and Brogan [9] the height of the domain is set as 155 mm and the width is set as 139 mm; and for the runs in Wan and Fell [1] the height is set as 300 mm and the width is set as 300 mm. The inlet boundary is imposed with a stable, consistent water head at each stage of hydraulic gradient. Taking Run B in Skempton and Brogan [9] as an example, the loading history is shown in Figure 7b.

In 2017, Rochim et al. [32] noted that the loading history might affect the experimental results of critical hydraulic gradient as well as the eroded mass. They found that if they set a longer duration time under a specific hydraulic gradient, they would record a lower critical hydraulic gradient and more eroded mass.

Such a phenomenon can be explained: when the duration time is not long enough under a specific hydraulic gradient, the erosion process may not reach a steady/stable status, which will lead to the unreliable judgment of critical hydraulic gradient.

To eliminate the effect of different loading histories, in our simulation, a long enough duration time and corresponding stable results are used to validate our model.

**Figure 7.** (**a**) Simulation domain and (**b**) loading history of hydraulic gradient of Run B by Skempton and Brogan [9].

#### *3.3. Parameter Calibration*

During the simulation, the pure water density ρ*<sup>w</sup>* is set as 1000 kg/m3, and the viscosity coefficient <sup>μ</sup>*<sup>w</sup>* is set as 10−<sup>6</sup> Pa·s. Besides GSD, other initial conditions of the soils (including initial porosity *n*0, initial permeability coefficient *ks*0, and density of soil particles ρ*s*) are also reproduced according to the experimental records.

There are two ways to obtain the initial porosity *n*<sup>0</sup> of the soil: (1) directly from the experimental records, and (2) by calculating *n*<sup>0</sup> = 1 − *V*s0/*V*, where *V*s0 denotes the initial soil volume and *V* denotes the initial volume of seepage field.

Parameters need to be calibrated include α in Equation (5), *b* in Equation (7), *A* and *B* in Equation (8), and Γ in Equation (10). The calibrated parameters are presented in Table 1. These parameters are calibrated using a part of experimental data, i.e., Runs A and B from Skempton and Brogan [9] (for non-cohesive soils) and Runs 2 and 3 from Wan and Fell [1] (for cohesive cases).


**Table 1.** Details of parameter calibration.

Based on these runs, *A*, *B*, *b* and α are calibrated by trial and error. *A*<sup>0</sup> = 0.021, *B*<sup>0</sup> = 0.015, and *b*<sup>0</sup> = 0.067, as suggested by Wilcock and Crowe [33] are used as initial values for *A*, *B*, and *b*, and α<sup>0</sup> = 0.0000021, as suggested by Fujisawa et al. [15], is used as the initial value of α.

It should be noted again that Wan and Fell [1] did not record the amount of soil erosion, and α ranges from 0.0001 to 0.01 would not affect the calculation results in Section 3.4.1. However, Skempton and Brogan [9] did record the amount of soil erosion during the experiment, the model can reasonably reproduce the *q-i* relationship as well as the eroded mass when α = 0.005 is subtracted. Therefore, α is fixed as 0.005 in our model. Value of Γ is calculated based on Equation (14):

$$\Gamma = \frac{k\_{\! 0}}{\left(d\_{10}^{1.565} \frac{n\_0^{2.3475}}{\left(1 - n\_0\right)^{1.565}}\right)}\tag{14}$$

#### *3.4. Calibration Results*

#### 3.4.1. *q-i* Relationship

Figure 8 presents the calibration results for runs in Section 3.3:

**Figure 8.** Calibration: dashed lines for computation results of runs in Skempton and Brogan [9]; solid lines for computation results of runs in Wan and Fell [1]; Solid points for test results of runs in Skempton and Brogan [9]; and hollow points for runs in Wan and Fell [1].

It shows that the computed *q*-*i* relationships agrees well with measured ones. It should be noted that ±10% change of these parameters will not affect the calibration results.

#### 3.4.2. Eroded Mass

The computed eroded mass in Run A (duration time, 80 min; last hydraulic gradient, *i* = 0.28) and Run B (duration time, 90 min; last hydraulic gradient, *i* = 0.34) in Skempton and Brogan [9] are 220.40 g and 77.16 g, which are close to the measured mass 225 g and 75 g, with relative errors less than 3%.

#### *3.5. Model Verification*

#### 3.5.1. *q-i* Relationship in Other Runs

Using the calibrated parameters, we conducted the simulation with the data from other runs, the results of representative runs are shown in Figure 9.

**Figure 9.** Verification: dashed lines for computation results of runs in Skempton and Brogan [9]; solid lines for computation results of runs in Wan and Fell [1]; solid points for test results of runs in Skempton and Brogan [9]; and hollow points for runs in Wan and Fell [1].

The internally unstable runs (e.g., Runs A and B) and internally stable runs (e.g., Run C) in Skempton and Brogan [9] show obvious difference in shape of *q*-*i* curves. For example, when *i* = 0.28 and 0.34, curves of Runs A and B show a sudden rise, respectively, whereas the curve of Run C shows a stable trend all over the experiment. The explanation is as follows: for Runs A and B, as particles are massively eroded by the seepage flow under a specific hydraulic gradient *i*, the porosity *n* and the hydraulic conductivity *ks* increases, as well as the discharge *q*; while for Run C, less soil particles are eroded under such a hydraulic gradient which triggers massive suffusion in Runs A and B, or an even higher hydraulic gradient. The simulation result of the data by Wan and Fell [1] show a similar phenomenon.

#### 3.5.2. Eroded Mass for Run C in Skempton and Brogan

The computed eroded mass for Run C in Skempton and Brogan [9] is 57.4 g, close to the measurement of 55 g. The relative error is less than 5%.

#### *3.6. Determination of Critical Hydraulic Gradient*

When *i* reaches critical hydraulic gradient *ic*, it can be inferred that the soil has undergone a sudden, irreversible, and largely destructive adjustment due to suffusion. Run 10 in Wan and Fell [1] is used as an example to reveal the reason. In Figure 10, the suffusion process can be divided into three different stages according to the different slopes of the *q*-*i* curve:

**Figure 10.** Three stages of suffusion process for Run 10 in Wan and Fell [1].

In stage I (*i* < 0.55), discharge of seepage flow *q* grows linearly with the hydraulic gradient *i*, during this stage the out flow is "clear" according to the corresponding experimental record; and in stage II (0.55 < *i* < 0.76), *q* grows in a nonlinear, but controlled, trend with *i*, during this stage the outflow is "cloudy" according to the corresponding experimental record; and in stage III (*i* > 0.76 = *ic*), one may notice an abrupt change of *q*-*i* curve, and this abrupt change may indicate the failure of soil fabrics.

When *i* < *ic*, the suffusion behavior has already taken place in the formation of losses of fine particles in limited scale, and such losses will not affect the stability of the whole soil fabrics. When *i* = *ic* or *i* > *ic*, more particles are lost until the failure of whole structure occurs.

Following the justification, *ic* is evaluated by observing the boundary between stage II and stage III for all the runs in the three studies. Figure 11 compares the computed critical gradients with the measured gradients. 18 of the total 21 runs are reproduced satisfactorily (relative error <10%). The three exceptions are Runs 5, 11, 13 in Wan and Fell [1].

**Figure 11.** Comparison between computed critical hydraulic gradient and experimental critical hydraulic gradient.

#### *3.7. Porosity Distribution in Di*ff*erent Su*ff*usion Stages*

In this study, taking Run A in the study by Skempton and Brogan [9] as an example in order to describe the suffusion process more directly and concretely, we present the spatial distribution of porosity corresponding to stage I, II, and III in 3.6 (see Figures 12 and 13), respectively:

Figure 12a denotes the initial porosity distribution (corresponding to stage I), Figure 12b,c denotes the porosity distribution when it reaches a stable one under hydraulic gradient *i* = 0.1 and *i* = 0.2 respectively (corresponding to stage II since *i* < *ic* = 0.28). When *i* = 0.1, the high-porosity zone only appears limitedly in the local region near the outlet, which means that there is no significant erosion in the seepage field. When *i* = 0.2, the high-porosity zone expands for a while, however, it is still limited to the local region near the outlet.

**Figure 12.** Porosity distribution for Run A in Skempton and Brogan [9] when *i* < *ic* = 0.28.

Figure 13a–c denotes porosity distribution when *i* = *ic* = 0.28 (corresponding to stage III), during the loading time *t* = 3000 s, 4000 s and 5000 s. The high-porosity zone grows so rapidly and intensely that it connects the upstream to the downstream. Such a phenomenon indicates an unsteady status under the massive suffusion behavior.

**Figure 13.** Porosity distribution for Run A in Skempton and Brogan [9] when *i* = *ic* = 0.28.

#### **4. Predicting Internal Stability**

In Section 3, we describe the suffusion process numerically and evaluate the critical hydraulic gradient *ic*, and in this section we attempt to give a further investigation in predicting internal stability.

#### *4.1. Shear Stress Comparison Approach*

Averagely, suffusion will be triggered when τ > τ*cr*, where τ*cr* denotes the critical shear stress of the finest soil particles (particle size *dr*) remaining in soil fabrics. According to Equation (7), τ*cr* = min(τ*ck*) as it has been already known that *dr* = min(*dk*), and it will increase when suffusion proceeds. Meanwhile, as the particles are continuously eroded out of the seepage field, the porosity ρ and hydraulic conductivity *ks* of the soils become higher, and based on Equation (10), τ will also increase. When *i* = *ic*, τ will keep exceeding τ*cr* as the suffusion continues. The suffusion will develop inevitably until the whole fabrics are destroyed.

For each of the soils, there must exist a hydraulic gradient *ik* < *ic*, when *i* = *ik*, if the loading time is long enough, and all the particles with sizes *d* < *dk* will be eroded out of the soil fabrics. Under this assumed situation, *dr* = *dk*, and τ*cr* = τ*ck*. At the moment, the volume of the soils *Vs* will be decreased by *k*% (ρ*<sup>s</sup>* is a constant), and the porosity of the soil will be *n* = *n*<sup>0</sup> + (1 − *n*0)*k*%.

Therefore, when *n* = *n*<sup>0</sup> + (1 − *n*0)*k*% is given, τ*cr* = τ*ck*; where τ*ck* is determined by Equations (7) and (8). Since then, a quantitative τ*cr*-*n* relationship can be computed. A τ-*n* relationship can also be obtained through Equation (6).

Here, an approach to predict *S* (stable) and *U* (unstable) based on the comparison of the relationship between τ*cr*-*n* and τ-*n* is illustrated: when *i* is given, by comparing the relative locations of the τ*cr*-*n* curve and the τ-*n* curve in the same coordinate system, we can judge whether suffusion starts, stops or proceeds irreversibly.

#### *4.2. Evaluation Results*

Taking Run A in Skempton and Brogan [9] as an example (see Figure 14), under the hydraulic gradient *i* = 0.28, when *n* < 0.405, τ > τ*c*, which means suffusion is triggered and developed; when *n* > 0.405, the τ-*n* curve almost coincides with the τ*c*-*n* curve, which means that the suffusion process is irreversible. When *i* = 0.30, the τ-*n* curve is always above the τ*c*-*n* curve and there is no intersection between the two curves. This indicates that under this hydraulic gradient, the suffusion process of the soil is still irreversible.

**Figure 14.** Change of shear stress and critical shear stress with porosity (τ*-n* and τ*c-n* curves).

According to the results shown in Figure 14, the critical hydraulic gradient is *ic* = 0.28, which is in accordance to the experimental data of Skempton and Brogan [9], and since *ic* < 0.7, according to the study by Hillel [34], the internal stability of soils in this run is *U*. It should be noted that α does not involve with the analysis process, which indicates that it does not affect the results of internal stability prediction.

To verify the approach, experimental data from Lau [36], Lafleur et al. [8], and Chapuis et al. [10] are also studied. The GSD curves are presented in Appendix B. This approach is able to distinguish the internal stability of 30 runs out of 36, indicating 83.33% of accuracy, while the most accepted traditional GSD-based approach based on the study by Wan and Fell [1] is 75%.

#### **5. Conclusions**

In this study, a new numerical model taking combined effects of GSD and porosity into account is developed to describe the suffusion process. The theory by Wilcock and Crowe [33] in field of sediment transportation is adopted to quantitively describe the inception and transport of finer soil particles. During the simulation, all of the related parameters are calibrated in a solid way. The numerical results showed satisfactory performance when validating with the experimental data of Skempton and Brogan [9] and Wan and Fell [1].

The model can reproduce the *q*-*i* relationships of experiments as well as critical hydraulic gradients *ic* of runs, and describe the suffusion process by showing porosity spatial distribution at different stages.

Inspired by the modeling work, an approach to predicting internal stability is also proposed. Disclosing the relationships among the factors, we found that the soils with specific GSD and porosity are prone to induce the suffusion. The model demonstrates better performance in the judgments of the internal stability (83.33% of correctness) of 36 experimental runs. The accuracy is higher than most of traditional GSD-based approaches. The results provide useful guidelines in macroscopic engineering practice.

**Author Contributions:** Conceptualization, Q.F. and X.F.; methodology, Q.F. and H.-C.H.; software, Q.F. and T.M.; validation, Q.F., Y.J. and J.W.; numerical runs, Q.F. and X.F.; investigation, Q.F.; data curation, H.-C.H. and J.W.; writing—original draft preparation, Q.F.; writing—review and editing, Q.F., Y.J. and X.F.; visualization, Q.F. and T.M.; supervision, X.F.; project administration, X.F.; funding acquisition, X.F. and Y.J.

**Funding:** This research is funded by The National Key Research and Development Program of China (Grant No. 2017YFC0804602), and The National Natural Science Foundation of China (Grant No. 51525901). The second author acknowledges the support of open fund of State Key Laboratory of Hydroscience and Engineering, Tsinghua University.

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

#### **Appendix A**

The numerical method for assigning basic variables (seepage velocity *vi* and porosity *n*) within each cell are listed as:

$$v\_i^{l,t} = \frac{\int\_{\Omega^l} v\_i^t dV^l}{V^l} \tag{A1}$$

$$m^{l\downarrow} = \frac{\int\_{\Omega^l} n^t dV^l}{V^l} \tag{A2}$$

where Ω*<sup>l</sup>* and *V<sup>l</sup>* denote the region and volume of the *l*th cell, respectively. *vi l,t* and *nl,t* denotes the seepage velocity and porosity within *l*th during *t*th time step, respectively.

Based on Equations (A1) and (A2), the unknown variable porosity *nl,t*+<sup>1</sup> in governing Equation (1) can be solved numerically as:

$$n^{l,t+1} = n^{l,t} + dt \sum E\_k^{l} A\_k^t \tag{A3}$$

where *Ek <sup>t</sup>* denotes the erosion rate of *k*th size group section at *t*th time step, and *Ak <sup>t</sup>* denotes the erodible area per unit volume of *k*th size group at *t*th time step.

Adopting the derivation by Fujisawa et al. [21], the unknown variable concentration *Cl,t*<sup>+</sup><sup>1</sup> in governing Equation (2) can be solved numerically as:

$$V^l \frac{n^{l,t+1}\mathcal{C}^{l,t+1} - n^{l,t}\mathcal{C}^{l,t}}{dt} + q\_i^t t^l {}\_{i\bar{j}}\Delta s\_{\bar{j}}^l = V^l \frac{n^{l,t+1} - n^{l,t}}{dt} \tag{A4}$$

As shown in Figure A1, *tij <sup>l</sup>* and Δ*sj <sup>l</sup>* denotes the unit normal vector and length of the *j*th (*j* = 1, 2, 3 and 4) boundary segment of the *l*th cell. Also, item *qi <sup>t</sup> tij <sup>l</sup>* can be calculated through Equation (A5) using the FVS scheme:

$$\mathbf{q}\_{i}^{t}\mathbf{t}\_{ij}^{l} = \mathbf{C}^{l,t}\mathbf{v}^{l,t}\mathbf{t}\_{ij}^{l} + \mathbf{C}^{r,t}\mathbf{v}^{r,t}\mathbf{t}\_{ij}^{l} \tag{A5}$$

where superscript *r* = *r*1, *r*2, *r*3, and *r*<sup>4</sup> denotes the neighboring cells (shearing the boundary) of *l*th cell.

**Figure A1.** Geometry of the cells.

In the simulation, item ∂*Sr*/∂*t* is always 0 since the soil are fully saturated all the time. Similar to Equations (A4) and (A5), the unknown variable water head *hl,t* in Equations (3) and (4) can be solved numerically based on the following Equation:

$$k\_s^{l,t} \frac{\text{grad} h\_1^{r\_1,t} - 2 \text{grad} h\_1^{l,t} + \text{grad} h\_1^{r\_2,t}}{\Delta s\_1^l} + k\_s^{l,t} \frac{\text{grad} h\_2^{r\_3,t} - 2 \text{grad} h\_2^{l,t} + \text{grad} h\_2^{r\_4,t}}{\Delta s\_1^l} = 0 \tag{A6}$$

where,

$$\text{grad}h\_1^l = \frac{h^{r\_1,t} - 2h^{l,t} + h^{r\_2,t}}{\Delta s\_1^l} \tag{A7}$$

$$\text{grad}h\_2^l = \frac{h^{r\_3, t} - 2h^{l, t} + h^{r\_4, t}}{\Delta s\_2^l} \tag{A8}$$

#### **Appendix B**

The GSD curve of soils in Lau [36], Lafleur et al. [8], and Chapuis et al. [10] are shown in Figure A2:

**Figure A2.** GSD curves of soils in (**a**) internal stable runs; (**b**) internal unstable runs in studies by Lau [36], Lafleur et al. [8], and Chapuis et al. [10].

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

### *Review* **SPH Modeling of Water-Related Natural Hazards**

#### **Sauro Manenti 1, Dong Wang 2, José M. Domínguez 3, Shaowu Li 4, Andrea Amicarelli <sup>5</sup> and Ra**ff**aele Albano 6,\***


Received: 4 July 2019; Accepted: 6 September 2019; Published: 9 September 2019

**Abstract:** This paper collects some recent smoothed particle hydrodynamic (SPH) applications in the field of natural hazards connected to rapidly varied flows of both water and dense granular mixtures including sediment erosion and bed load transport. The paper gathers together and outlines the basic aspects of some relevant works dealing with flooding on complex topography, sediment scouring, fast landslide dynamics, and induced surge wave. Additionally, the preliminary results of a new study regarding the post-failure dynamics of rainfall-induced shallow landslide are presented. The paper also shows the latest advances in the use of high performance computing (HPC) techniques to accelerate computational fluid dynamic (CFD) codes through the efficient use of current computational resources. This aspect is extremely important when simulating complex three-dimensional problems that require a high computational cost and are generally involved in the modeling of water-related natural hazards of practical interest. The paper provides an overview of some widespread SPH free open source software (FOSS) codes applied to multiphase problems of theoretical and practical interest in the field of hydraulic engineering. The paper aims to provide insight into the SPH modeling of some relevant physical aspects involved in water-related natural hazards (e.g., sediment erosion and non-Newtonian rheology). The future perspectives of SPH in this application field are finally pointed out.

**Keywords:** SPH (Smoothed Particle Hydrodynamics); water-related natural hazards; sediment scouring; dense granular flow; fast landslide; surge wave; flooding on complex topography; HPC (High Performance Computing); FOSS (Free Open Source Software)

#### **1. Introduction**

Thanks to the availability of high-performance computers, in the last few years, computational fluid dynamics (CFD) has been widely applied to simulate natural hazards in the field of hydraulic engineering. Due to the fast and large deformations characterizing the problems in this research field, meshless techniques allow for some intrinsic limitations of traditional grid-based methods (e.g., mesh deformation and cracking; free-surface, and interface treatment) to be overcome. Among the meshless techniques, the smoothed particle hydrodynamics (SPH) method has proven to have advantages over other methods. Following a Lagrangian approach, each continuum is discretized through a discrete set of material particles that lack connective mesh and follow the deformation undergone by the material. The dynamics of material particles obeys Newton's laws of motion and the discretized form of the governing equations (i.e., momentum and mass balance) is obtained through particle approximation

using an interpolant kernel function (i.e., a central function of the particles' relative distance). Based on the solution algorithm of the discretized governing equations, two different approaches can usually be defined: weakly compressible SPH (WCSPH), if the continuum is assumed to be slightly compressible (governing equations can be decoupled), and incompressible SPH (ISPH).

Even if SPH was originally developed for astrophysical problems, it has been subsequently extended to the analysis of free-surface flows and several kinds of multiphase flows including non-Newtonian fluid and post-failure dynamics of rainfall induced landslide. Many works have shown that the multiphase coupled dynamics of water and non-cohesive sediment mixture may be simulated following a fluidic approach [1–3]. Through the proper definition of appropriate yielding criteria and non-Newtonian rheological models, the triggering and propagation of multiphase flow as well as dense rapid granular flow can be simulated through the SPH approach with a suitable degree of accuracy in reproducing the experimental results, showing that this approach can be conveniently applied to real case studies. In this way, it is possible to simulate scouring around complex structures that may cause hazards related to both structural instability and threats to the safe management of hydro-power plants [4–7]. Another important branch of non-cohesive sediment dynamics is related to two-phase flows on complex 3D topography such as the sediment eroded by a dam-break wave propagating on a mountain valley [8].

In the field of dense granular flows, fast landslide dynamics can be properly simulated by accounting for the interactions with water (both pore water and stored water) that influences landslide deformability and run-out characteristics. Thus, complex multiphase flows can be simulated for investigating natural hazard related to a tsunami wave generated by a quick landslide at the slope of a reservoir [9–13]. Additionally, a fast shallow landslide triggered by intense rainfall at the slope of a hill represents the most common natural hazards in some areas of the world [14,15] and can be simulated with proper numerical modeling. Furthermore, SPH can be successfully applied to the analysis of slope stability, with the detection of failure initiation and propagation leading to delineation of the potential slip surface [16].

Computational methods are also successful in simulating complex nonlinear flows and transport processes involved in natural hazards related to the flooding risk, as in the case of dam and dyke failures. Water wave propagation in complex 3D topography can be properly predicted [17–20] including urban areas where the complexity of the 3D flow pattern is increased by sudden changes in the propagation direction due to buildings as well as different types of floating bodies such as cars, trucks, etc., as in the works in [21,22].

Even if mathematical models have been widely established to simulate the relevant hydraulic and hydro-geologic processes involved in the simulation of water-related natural hazards [23–25], the relevant complexity (especially from a geometric point of view) of practical applications in the field of applied engineering poses some difficulties in the applications of the computer codes that implement the above-mentioned mathematical models. The numerical modeling of a complex 3D problem frequently requires the discretization of large domains with a high resolution, which implies simulating a very large number of elements (usually at least some millions), thus increasing the required computational time and resources in a non-linear way. Fortunately, the parallel computational power of the computer hardware has increased significantly in recent years, especially in the case of graphics processing units (GPUs) which are now several times more powerful than traditional central processing units (CPUs). However, it is important to note that the GPU is more powerful only for properly posed problems that are parallelizable on GPUs, and not necessarily for general computing. Furthermore, it is imperative that the models are implemented using high performance computing (HPC) techniques to take advantage of the power of current hardware (i.e., multi- cored CPUs, GPUs, or other hardware accelerators).

Another limitation that computational modeling suffers is the inherent uncertainties of some modeling parameters that influence (in some cases even significantly) the model response and lower the reliability of the results. Despite this intrinsic limitation, such deterministic models may be useful anyhow to support risk analysis and mitigation through the development of fast-running numerical models that can help to create probabilistic maps of risk variables. In this way, a multi-disciplinary decision support system for natural hazard risk reduction and management could be developed that allows for the exploration of several scenarios including potential risk-reduction options [26].

In this context, the aim of the paper was to provide an overview of some recent applications of the SPH method to the modeling of water-related natural hazards. The selected works illustrate some of the relevant problems of practical and theoretical interest in the field of hydraulic engineering that could be useful to provide an introduction to the SPH method as applied to the analysis and mitigation of water-related natural hazards.

Additionally, some widespread free and open-source software (FOSS) CFD codes for multiphase engineering applications are reviewed as well as their capabilities in modeling some of the problems described in this paper.

#### **2. Two-Phase Coupled Dynamics**

This section illustrates some engineering applications of the SPH method in the simulation of two-phase flows involved in hydraulic engineering problems of practical interest in the field of water-related natural hazards.

When dealing with the numerical simulation of two-phase flows of immiscible fluids, a central aspect to be considered is the interface treatment. Concerning incompressible SPH (ISPH) simulations of gravity currents induced by the inflow of a certain fluid into another fluid with a density difference, two different approaches were proposed in [27] and classified as a coupled model and a decoupled model, respectively. The coupled model is based on the assumption that standard ISPH governing equations can be applied to every particle in the computational domain regardless of the phase it belongs to. Therefore, when an interface particle is considered, the kernel interpolation is extended over all the heterogeneous neighbors in its interaction domain, which is the kernel compact support (circular in 2D problems or spherical in 3D problems) centered on the considered interface particle. In contrast, the decoupled model assumes that standard ISPH governing equations (including the Poisson equation) should be solved separately for each phase by identifying the free-surface. Therefore, if an interface particle is considered to be the kernel, interpolation is restricted to its homogeneous neighbors solely within the kernel support, thus kernel truncation occurs. Proper interface boundary conditions should be satisfied to couple both phases: these conditions can be represented by the continuity of both normal and tangential stresses between heterogeneous particles at the interface. According to the results in [27], the decoupled model also provides higher accuracy than the coupled model in the case of multi-fluid flow with a high-density ratio.

A modified formulation of the standard weakly compressible SPH (WCSPH) governing equations were presented in [28] for modeling rapid free-surface multiphase flows with a high-density ratio involving violent impact against a rigid wall. This formulation, which is based on the coupled approach, allows for the numerical instability induced by the discontinuity of the fluid properties across the interface to be eliminated without excluding the heterogeneous particle during kernel interpolation. Fairly accurate results were achieved with this coupled model simulating some flows with continuous pressure across the interface. Other good examples of multi-phase SPH formulations (both WCSPH and ISPH) that deal with the interface between fluids under certain constraints are described in the works in [29–32].

Among the variety of problems involving the numerical simulation of multi-phase flows, two topics of concern in the field of water-related natural hazards have been examined in the following: (i) scour with non-cohesive sediment transport induced by rapidly varied free-surface flow and erosion around complex structures; and (ii) fast dynamics of dense granular flows and landslides. Some relevant works on these topics are discussed, pointing out the peculiar aspects of numerical modeling.

#### *2.1. Scouring and Sediment Transport*

An important aspect of the design and maintenance of the bottom-founded submerged structure is non-cohesive sediment scouring. Bottom sediment erosion around the structure is caused by a complex flow pattern induced by the presence of the structure that strongly modifies the upstream undisturbed flow field [33]. The scouring evolution over time must be properly analyzed and mitigated in order to avoid worsening of structural stability due to foundation exposure.

This topic is widely investigated in fluvial hydraulics, in coastal areas, and in the marine environment. When dealing with a safety assessment of a hydraulic structure placed in the riverbed (e.g., bridges, barrages, etc.), the erosive action of the river stream must be reliably evaluated [34]. Several empirical formulas to predict final scouring are available, but the phenomenon is time-dependent and is affected by several uncertainties (related to both sediment and flow characteristics) as well as stochastic factors influencing the flow evolution such as transport and deposition of wood debris [5,35,36].

Additionally, in the design of highly demanding marine structures, the scour around the foundations should be dependably predicted [37,38]. For instance, a continuous one-way current due to the overtopping flow of an extreme tsunami wave may cause long-term erosion of the foundation on the rear side of the coastal protection structures [39,40]. Scour pit evolution behind the seawall induces the formation of eddies with increasing size to adapt the growing dimension of the scour cavity. Excessive sediment erosion directly results in the instability, and even destruction, of these hydraulic structures [41–44]. Most empirical relations have been established to predict the bulk and time-averaged sediment transport rate, but their application requires experimental data for calibration. Furthermore, these empirical relations were obtained from small scale laboratory experiments and their extension to a real scale problem may lead to significant erosion estimate errors [45].

An assessment of the mechanics of detailed temporal erosion processes as well as for the reliable design and assessment of these structures requires both physical investigations and advanced numerical modeling that accounts for the effects of the stochastic variables on the scouring process. Of course, numerical modeling is generally less expensive than physical experimentation, but, in some cases, experimental data may be available from technical literature. As with all engineering problems, any numerical tool must be properly validated against experimental data. After the model has been validated, its parameters require tuning based on the model/experiment results.

Many sediment modeling techniques use mesh-based approaches such as finite difference method and finite volume method to simulate the erosion and fluid-sediment coupled dynamics. However, these mesh-based methods suffer from some intrinsic limitations due to the fixed grid system, which leads to some difficulties in effectively simulating the bed-grain interactions, fluid-sediment momentum transmissions, and the dynamics within the deposits because of their fixed grid system. Furthermore, accurately tracking the free surface and fluid-sediment interface is also a big challenge for these approaches. Lagrangian meshfree particle methods overcome these limitations by intrinsically capturing the free-surface and tracking the particles. These methods have been widely used in recent years for the analysis of erosion processes.

The SPH technique has been proven to be effective in complex multiphase applications [46] such as water–gas flows or bubbly flow simulations [47,48]. In addition, when simulating non-cohesive bottom sediment scouring by rapid water flow, both weakly compressible (WCSPH) and incompressible (ISPH) formulations allow for a reliable description to be obtained. The basic idea is to model the sediment dynamics, likewise a pseudo-fluid, once the onset of sediment particle motion is attained. According to this approach, a criterion should be defined that represents the critical condition for the incipient motion of sediment; this can be done in terms of either critical velocity [49,50], or critical shear stress. The first approach may lead to some problems if critical velocity is evaluated as the depth-averaged value that in some cases may not be representative of the scouring, as for the continuous tsunami overflow behind a seawall that induces rapid water depth variation inside the scour pit [2]. The second approach based on the critical shear stress has been widely adopted to analyze non-cohesive sediment scouring by rapid water flows.

*Water* **2019**, *11*, 1875

The authors in [7] implemented two different criteria in a WCSPH model to simulate the erosion process of non-cohesive bed sediment due to constant bottom water discharge from 2D tank (flushing). These criteria are based on Mohr–Coulomb (M–C) yielding theory and on Shields (SH) theory, respectively. The M–C critical condition defining sediment failure and the onset of erosion requires the introduction of a maximum viscosity as the product between sediment viscosity, μ*s*, and a magnification factor, η, representing the numerical parameters to be tuned. When the apparent viscosity is lower than the maximum viscosity, the sediment is treated as a non-Newtonian fluid of Bingham type and solid particles are set in motion with constant viscosity μ*<sup>s</sup>* (green curve in Figure 1). The strategy of introducing an upper viscosity limit for the sediment was also followed in [51] in the WCSPH simulation of complex problems in the field of marine engineering; below this maximum limit, the work in [51] adopted a variable apparent viscosity calculated through the M–C theory for the soil phase. The SH critical condition does not require the introduction of a numerical threshold for the viscosity of the solid phase. However, in [7], both the M–C and SH approaches require tuning of the mechanical parameters of the bottom sediment such as the angle of internal friction, ϕ, and sediment viscosity, μ*s*, that became numerical parameters to fit the experimental eroded profile. This may be not be practical when calibration data are not available for the investigated problem.

**Figure 1.** Rheological models for non-cohesive sediment erosion with bed load transport and for dense granular flows. *I*2: second invariant of the rate of deformation tensor; μ*app*: apparent viscosity.

In order to overcome this limitation, [8] proposed a WCSPH formulation of a mixture model for the analysis of dense granular flows consistent with the kinetic theory of granular flow (KTGF). This mixture model, which avoids the use of an erosion criterion, has been integrated into the FOSS code SPHERA v.9.0.0 (RSE SpA). The relevant physical properties (i.e., density and velocity) of the mixture of pure fluid and granular material are expressed as a function of the volume fraction ε, occupied by each phase at a material point. The balance equations for the mixture are defined accordingly and discretized consistently with the WCSPH approach. The frictional regime of the mixture dynamics is represented under the packing limit of the KTGF, which holds for the volume fraction ε*<sup>s</sup>* of the solid (granular) phase and is characteristic of bed-load transport and fast landslides (see also Section 2.2). In the frictional regime, the mixture (or apparent) viscosity, μ, is calculated as a weighted sum of the pure fluid viscosity μ*f*, and the frictional viscosity μ*fr*, the latter being evaluated on the basis of the mean effective stress σ'*m*, angle of internal friction ϕ, and the second invariant *I*<sup>2</sup> of the rate of the deformation tensor of the sediment. The frictional viscosity increases as the shear rate tends to zero, in accordance with the pseudo-plastic rheological behavior (dashed blue curve in Figure 1). To avoid the unbounded growth of apparent viscosity of the mixture, a threshold (or maximum) viscosity μmax is introduced with a physical meaning. Threshold viscosity acts when approaching the zero shear rate; mixture particles with an apparent viscosity higher than the threshold viscosity are considered in the elastic–plastic regime of soil deformation where the kinetic energy of solid particles is relatively small and the frictional regime of the packing limit in the KTGF does not apply. For these reasons, the threshold viscosity is assigned to those particles that are excluded from the SPH computation (continuous red curve in Figure 1; below μmax, the red curve coincides with the dashed blue curve of the pseudoplastic model). The excluded particles represent a fixed boundary with suitable values of the relevant physical properties and are included in the neighbor list of the nearby moving particle. The value of the threshold viscosity does not require tuning or calibration, but it should be selected for the specific problem; the value assigned to the threshold viscosity is the higher value that does not influence the numerical results significantly. Further increase of the value assigned to the threshold viscosity only affects the computational time because it determines an increase in the maximum value that can be assumed by the apparent mixture viscosity during the computation. The relationship between the mixture viscosity and the time step value that assures the stability of the adopted explicit time-stepping scheme is given by Equation 2.33 in [8]. When the numerical stability of the time integration scheme is dominated by the viscous criterion, the threshold viscosity reduces the computational time.

In [8], some experiments of erosional dam breaks adopting the physical values for the mechanical parameters of the sediments including the angle of internal friction were simulated with good reliability. The 2D experimental tests, involving rapidly varied mixture flows with erosion and the transport of bed sediment, were simulated in [7] and [8]. Moreover, in [6], mixture flow computation was performed with the open-source DualSPHysics solver [52] accelerated with a graphic processing unit (GPU). The adopted WCSPH algorithm, representing the improvement of the model in [53], reproduces the dynamics of the bottom sediment phase combining two erosion criteria with the non-Newtonian rheological model. The adopted rheological model is based on the Bingham-type Herschel–Bulkley–Papanastasiou (HBP) model providing the apparent viscosity μ*HBP* of the sediment:

$$\mu\_{HBP} = \frac{|\tau\_c|}{\sqrt{l\_2}} \left( 1 - e^{-m \sqrt{l\_2}} \right) + 2\mu \left( 2\sqrt{l\_2} \right)^{n-1} \tag{1}$$

This model allows for the transition between un-yielded (fixed) and yielded (mobile) granular material to be described without introducing a maximum value of the viscosity for the solid phase. Proper selection of the HBP model's parameters, *m* and *n*, allows for the reproduction of the required rheological behavior as the Newtonian or the pseudoplastic model. The yield stress parameter τ*<sup>c</sup>* in the HBP model was evaluated in two different calibration procedures depending on the erosion criterion that holds for modeling the yielding mechanism of sediment. The qualitative representation of Equation (1) with *n* < 1 is denoted by the blue dashed curve in Figure 1.

If a solid particle is identified as being near the sediment–water interface by means of the criteria described in [6], then the Shields criterion determines the onset of bed erosion at the sediment–water interface and provides the critical bed shear stress τ*bcr*, replacing the yield stress parameter τ*c*, in the HBP model.

If the particle is near the water–sediment interface, but Shields erosion criterion is not satisfied or if the particle is far from the water–sediment interface, then the yielding mechanism is modeled through the Drucker–Prager criterion that allows the critical shear stress τ*y*, replacing the yield stress parameter, τ*c*, in the HBP model to be defined to determine the apparent viscosity for the sediment.

As a result, three distinct regions may be defined, starting from the water–sediment interface and going downward: (a) eroded sediment exceeding the Shields critical bed shear stress (bed load transport); (b) yielded sediment (plastic deformation with slow kinematics); and (c) un-yielded sediment (high viscosity, static condition).

The model proposed in [6] also accounts for suspended transport. Following [51], the identification of the suspension layer is obtained through a non-dimensional concentration, *cv*, computed for an interface particle as the ratio of the sediment particle volume to the total particle volume within the interaction domain. The onset of suspended transport is determined by *cv* falling below the threshold value of 0.3 and the suspended sediment viscosity is computed through [54] the experimental colloidal equation, which is more simple to implement with respect to the piece-wise function adopted in [51]. The density of suspended particles is computed by solving the mass balance equation. Even if, in some cases, the percentage gap between the experimental and numerically predicted maximum scouring depth is significant, it can be seen that the scour process is affected by several random factors and therefore reliable predictions of scouring effects are quite difficult to obtain, even with experimental modeling. The sediment dynamics models based on a synthetic rheological law (e.g., [6]) assume the same rheological behavior for the bed-load transport (frictional regime of KTGF), suspension for dense granular flows (kinetic-collisional regime of KTGF),a and suspension for diluted granular flows (kinetic regime of KTGF). This feature provides advantages and drawbacks with respect to KTGF-based sediment dynamics models (e.g., [8]). The model of [6] can reproduce several sediment transport regimes (not only bed-load transport), but is not coherent with KTGF, some parameters require tuning procedures, and non-transported granular material (e.g., landslides) is not treated.

The work in [2] investigated by means of ISPH the water-induced 2D sediment scouring where the erosion process is mainly related to loose sediment particles suspended in the water flow, as in the case of scouring behind a seawall produced by continuous tsunami overtopping. In contrast to the physics of sediment flushing, where sediment moves as bed load at very high density, in the case of overtopping erosion, the solid particles move more loosely and the density of turbid water is significantly lower than the sediment density. In such situations, the erosion dynamics is controlled by the balance between the suspension effect due to turbulent mixing and settling of suspended solid particles owing to gravity force. This process is strongly affected by the size of numerical particles, which is usually far bigger than the size of a real sediment grain in practical problems. For the reasons explained above, adoption of real sediment density in the computation may lead to an unreliable representation of the erosion. The model proposed in [2] introduces a simplification due to the size (and hence number) of the particles required to properly model the turbulent mixing. This model is based on the concepts of numerical turbid water particle (TWP) and clear water particle (CWP) to simulate the sediment-entrained flow in cases where the granular particles of sediment move more loosely and sediment is washed away, mainly in the form of a suspended load. Due to the size considerations discussed above, numerical particles should be treated as a combination of clear water and turbid water particles if a suspended load is simulated; therefore, eroded solid particles represent a water–sediment mixture whose density is calculated based on the integral interpolation theory, thus accounting for density reduction as solid particles are suspended and mixed with water. The value of 1250 kg/m<sup>3</sup> has been suggested for the reasonable initial density of the TWPs based on studies of bottom sediment movement. Figure 2 shows that the proposed ISPH model can simulate the real-time processes of the 2D overflow induced scouring. The detailed comparisons between numerical and experimental data can be found in [2].

**Figure 2.** Snapshots of incompressible smoothed particle hydrodynamics (ISPH) calculated equilibrium scouring pit with different overflow depth: (**a**) η = 3.3 cm; (**b**) η = 4.7 cm, and (**c**) η = 6.0 cm.

In the subsequent work in [1], the ISPH model of [2] was successfully extended to the simulation of 3D local bed scour induced by clear water stationary flow passing a submerged vertical cylinder of relatively large size. Additional formulations were introduced to account for transversal and longitudinal bed slope. The erosion model was based on the turbidity water particle concept and the sediment motion was initiated when the calculated shear stress on the interface particles exceeded the critical value. The 3D ISPH erosion model was used to simulate the scouring process around a large vertical cylinder with a diameter of 60 cm in Figure 3, where the vorticity and the shape of the scour pit are illustrated at time *t* = 2.5 s and *t* = 3.5 s. The numerical results show that the proposed ISPH model could simulate the relevant features of the flow and the scour process (Figure 3). The detailed comparisons between the numerical and experimental data can be found in [1] and the scour dynamics were validated with a suitable degree of accuracy for engineering purposes. Even if a suitable representation of the complex scouring process can be obtained, the vorticity field shows some numerical noise (Figure 3a). Improvement of the calculated vorticity field could be obtained through the approach proposed in [46].

**Figure 3.** ISPH computed sediment bed scouring process: (**a**) snapshots of ISPH calculated vorticity in the scouring pit; (**b**) Snapshots of the ISPH calculated scouring pit.

#### *2.2. Fast Landslides and Dense Granular Flows Interacting with Water*

Numerical modeling of dense granular flows and landslides is still a challenging topic, especially when considering the interaction between the sediment and the water that may be both an internal interaction, related to pore water in landslide-prone saturated soil, and an external interaction with stored water in a basin with unstable slopes.

The interaction between pore water and soil matrix is a fundamental aspect that influences the triggering and propagation of shallow landslides induced by intense rainfall events that represent the most common natural hazards in some areas of the world [14]. Intense rainfall events induce water infiltration at slopes, leading to an increase of the volumetric water content and pore water pressure that worsen the slope stability of the soil layer close to the surface and may cause its failure. Reliable assessment of landslide susceptibility also requires a proper definition of the rainfall characteristics considering recent climate trends affecting rainfall [55,56].

If landslide triggering occurs, in the post-failure phase, the pore water content, combined with geo-mechanical properties of the soil, influence the sediment dynamics, and may induce in some cases flow-like fast earth movements classified as complex landslides because their run-out start as shallow rotational-translational failure, but changes into dense granular flows due to the large water content [15]. In this case, the fast landslide dynamics is more difficult to predict using traditional analytical models and a numerical approach could be helpful for landslide susceptibility assessment and the creation of debris–flow inundation maps.

The work in [57] proposed a combined triggering–propagation modeling approach for the evaluation of rainfall induced debris flow susceptibility. They adopted the transient rainfall infiltration and grid-based regional slope-stability model (TRIGRS) [58], which is based on the infinite slope stability approach, to obtain a map of potentially unstable cells within a study catchment under an intense rainfall event. Since not all unstable cells, in general, evolve to a debris flow, an empirical instability-to-debris-flow triggering threshold is calibrated on the basis of observed events and used to identify, among unstable cells, the so-called triggering cells that could most likely contribute to debris flow. The triggering threshold is, therefore, the key element that allows coupling between the TRIGRS slope instability model and the debris flow propagation model FLO-2D [59], a finite volume model that numerically solves the depth-integrated flow equations. The work in [57] assumed a zero excess rainfall intensity and a total friction slope depending on the Bingham-type rheological parameters as a function of the sediment concentration. The calibration of the triggering threshold with geo-morphological data of the catchment area represents a crucial step for obtaining reliable susceptibility maps in the nearby areas. Back-analysis of a catastrophic event that occurred on 1 October 2009 in the Peloritani mountain area (Italy) provided fairly good results.

In order to quantify the level of risk addressing the uncertainty inherent in landslide hazard, the susceptibility evaluation represents only one of the relevant issues involved in risk analysis. Additionally, run-out dynamics should be properly assessed in order to provide a quantitative estimate of the landslide hazard and select appropriate protective measures for risk mitigation. In order to reach such a goal, reliable predictive models should be used to obtain quantitative information on the destructive potential of the landslide. This is mainly related to the following characteristics: run-out distance; width of damage corridor; travel velocity; characteristic depth of the moving mass; and characteristic depth of the deposits [60].

The work in [61] proposed a 2D depth-integrated, coupled, SPH model for predicting the path, velocity, and depth of flow-like landslides. While the post-failure flow model of [57] assumes the heterogeneous moving mass as a single-phase continuum, [61] modeled the dense granular flow as a two-phase mixture composed of a solid skeleton with the voids filled by a liquid phase. Assuming that the shear strength of the liquid phase can be disregarded, the stress tensor within the mixture is composed of pore pressure and effective stress. Then, the mixture dynamics was described by quasi-Lagrangian depth integrated governing equations of mass and momentum balance, and the pore pressure dissipation equation that were discretized according to the standard SPH approach. Depth-integrated equations do not provide information on the vertical flow structure that is needed for evaluating shear stress on the bottom and depth integrated stress tensor. For this reason, [61] assumed that this vertical structure would be the same as the uniform steady-state flow according to the so-called model of the infinite landslide having constant depth and moving at a constant velocity on a constant slope. The work in [61] adopted the simple method proposed in [62] to obtain the bottom shear stress in a non-Newtonian fluid of Bingham-type. The model in [61] was used to simulate some catastrophic event that occurred on May 1998 in the Campania region (Italy), showing the relevant role of geotechnical parameters (especially the fluid phase and angle of internal friction) for the reliable prediction of the run-out distance, velocity, and height of the landslides. Proper selection of the value assigned to these parameters assured the best agreement with the field observations.

Some flow-like landslides are characterized by a relatively small average depth if compared with the horizontal linear dimensions, therefore the assumption of a depth-integrated model is strongly consistent with the physics of the phenomenon. Furthermore, there could be rainfall induced landslides where the initial average depth is comparable with the horizontal length and width. In these cases, significant variations of the vertical thickness and the vertical velocity profile may occur along the landslide body in the flow direction. This is illustrated in Figure 4, showing the preliminary results of an on-going study. In this new study, a narrow landslide that occurred during an intense rainfall event on April 2009 in a hilly area of the Oltrepò Pavese, named the Recoaro Valley, Northern Italy, was reproduced by means of the 2D WCSPH simulation. It can be seen that in the early phase, just after the failure, the mass portion close to the landslide front moved faster than the rear portion. Additionally, just after the impact against the downstream vertical wall of a damaged building, the landslide front began decelerating while the rear mass portion on the steep slope still maintained a relatively high average speed.

**Figure 4.** Frames of the 2D smoothed particle hydrodynamics (SPH) simulation of the rainfall induced shallow landslide occurred on April 2009 at Recoaro, Oltrepò Pavese (Northern Italy). (**a**) Velocity field; (**b**) Density field.

In the peculiar case described above, the modeling approach based on the infinite landslide with constant depth and moving at a constant velocity on a constant slope may be less appropriate. Therefore, the solution of the governing equation in the three-dimensional form seems more appropriate than the depth-integrated model. The simulation shown in Figure 4 was carried out with the code SPHERA v.9.0.0 adopting the mixture model for dense granular flow discussed in [8]. Even if the code has a 3D formulation, a 2D approach may be conveniently adopted in this case because the landslide is relatively narrow and the flow may be assumed to be identical on the vertical planes along the flow direction.

The work in [63] used a finite volume approach to simulate mudflows and hyper-concentrated flows characterized by suspended fine material by adopting a Bingham rheological model. Similar to [64], the Cross model was adopted for modeling the non-Newtonian flow, assuming that the constant parameter *m* was equal to 1, resulting in the following formulation of the apparent viscosity:

$$\begin{cases} \mu\_{\varepsilon f\bar{f}} = \frac{\mu\_0 + \mu\_{\infty} K \bar{\cdot} \bar{\cdot}}{1 + K \bar{\cdot}}\\ K = \frac{\mu\_0}{\tau\_B} \mu\_{\infty} = \mu\_B \end{cases} \tag{2}$$

In Equation (2) . γ denotes the shear rate defined through the second invariant *I*<sup>2</sup> of the rate of deformation tensor; *K*, μ0, and μ<sup>∞</sup> are three constant parameters that can be conveniently related to common Bingham rheological parameters, namely the yield stress τ*<sup>B</sup>* and viscosity μ*B*. From a physical point of view, μ<sup>0</sup> and μ<sup>∞</sup> denote the viscosity at very low and very high shear rate, respectively. In order to avoid numerical divergence caused by the unbounded growth of effective viscosity as the shear rate approaches zero, μ*e*ff is limited to a suitably high threshold value, which is set to μ<sup>0</sup> = 103μ*<sup>B</sup>* to assure convergence. The test cases simulated in work in [63] considered the flow on inclined surfaces and analyzed the role of the Froude number [65] during the propagation phase, which may be helpful in designing the control works.

Landslides occurring at the slopes of confined water bodies (e.g., artificial basin or river-valley reservoir) or at coastal regions involve complex interactions between the solid and the fluid phase. In the post-failure phase of rapid landslides generated at the slope of a water body, an impulse wave is generally induced, which is usually referred to as tsunami [66,67]. The characteristics of the generated water wave are related to the velocity and the shape (e.g., thickness and slope angle) of the landslide front, whose dynamic deformation is in turn affected by the water induced stresses at the interface while the landslide is penetrating the water body. Accurate prediction of the landslide induced wave hazard depends on reliable numerical models that can simulate the coupled dynamics.

The work in [11] proposed a WCSPH model for the analysis of impulsive wave generated in a basin by a deformable landslide. In order to properly reproduce landslide deformation during the post-failure phase, they adopted an elastic–plastic constitutive model for the soil combined with the Drucker–Prager criterion. To account for the interaction with stored water, occurring at a very small timescale due to the fast dynamics, [11] introduced a bilateral coupling model consisting of two sequential steps: the interface soil particles were initially considered as the moving deformable boundary whose velocity and position was used for solving the governing equations of the water phase; then the soil constitutive equation was solved with the corrected stress taking into account the water-induced surface force. This coupling model would require, in theory, an iterative procedure to assure the consistency of both stress and deformation at the interface at each time step. However, the sequential approach is quite suitable because the interface deformation is mainly caused by the landslide dynamics and small displacements occur within a time step. The mechanical parameters in the model of [11] have a physical meaning and the corresponding values can be deduced from conventional soil mechanic experiments. Two experiments were simulated concerning the wave generated by a slow landslide [12] and a fast landslide [68], respectively, obtaining in both cases good agreement with the experimental data.

The work in [10] proposed a hybrid model for simulating the coupled dynamics between the landslide and stored water as well as the propagation of generated surge waves. The discrete element method (DEM) was used for simulating the landslide dynamics, while WCSPH was used to solve the governing equations for water. Spurious numerical noise affecting the water pressure field was removed by applying the δ-SPH technique [69]. The interaction mechanism between the solid and fluid phase in the hybrid DEM-SPH model was based on the drag force and buoyancy. At each time step, these forces were first calculated considering (i) the initial position and velocity of both fluid and solid particles, (ii) the pressure of the fluid, and (iii) the local soil porosity, ε, evaluated at the landslide–water interface with a kernel interpolation of the solid particle volume. Based on the calculated drag force and buoyancy, the corresponding new position and velocity of both fluid and solid particles were then calculated to solve the corresponding discretized governing equations. The updated position and the velocity field were subsequently used to recalculate drag force and buoyancy. Therefore, this interaction mechanism requires an iterative process to assure convergence of the calculated forces and consistency of the displacements and velocity field of the involved phases. The hybrid DEM-SPH model was used to simulate the sliding along a 45◦ sloping plane of a rigid body that mimics a submarine landslide, obtaining a suitable prediction of the experimental time evolution of water surface elevation [70]. A modified simulation was additionally performed by considering the deformability of the sliding body, showing that smaller and less violent surge waves were generated owing to the landslide deformation.

A similar approach was proposed in [71] for the analysis of underwater granular collapse by coupling WCSPH (for the fluid phase) with DEM (for the solid incoherent phase). A coupling module was developed for fluid–grain interaction: the force exerted by the fluid on a solid particle was obtained by integrating the contributions on its surface; in turn, the effect of the solid particles on the fluid motion was calculated by including the neighboring DEM particles in the SPH interpolation of the governing equations for the considered liquid particle. Another attempt of coupling SPH, in particular, the DualSPHysics model, with a distributed-contact discrete-element method (DCDEM) was proposed

in [72] to explicitly solve the fluid and solid phases to model a real case of an experimental debris flow. An experimental setup for stony debris flows in a slit check dam was reproduced numerically, where solid material was introduced through a hopper, assuring a constant solid discharge for the considered time interval.

The reference mixture model of [8] for the dynamics of dense granular flow was modified by introducing a numerical parameter, μ0, referred to as limiting viscosity. The effect of limiting viscosity arises in the frictional regime at low deformation rates near the transition zone to the elastic–plastic regime: in this shear rate interval, a constant value μ<sup>0</sup> (lower than μ*max*) is assigned to the mixture viscosity (see red dot-dashed curve in Figure 1), thus reducing the computational time in the case where the viscous criterion dominates the numerical stability of the time integration scheme. There are other alternative approaches to keep control of computational time in the simulation of high-viscosity flows. In the work of [73], a semi-implicit integration scheme was proposed to overcome the severe time-stepping restrictions caused by the WCSPH explicit integration scheme when simulating highly viscous fluids, as in the case of lava flow with thermal-dependent rheology. According to this approach, only the viscous part of the momentum equation is solved implicitly, thus saving computational time and obtaining an improved quality of the results with respect to the fully explicit scheme.

The mixture model of [8], modified with limiting viscosity, was successfully applied in [9] to the analysis of a fast massive landslide at the slope of an artificial reservoir. The simulated case reproduced a two-dimensional scale laboratory test carried out in 1968 at the University of Padua reproducing in Froude similitude a characteristic cross-section of the Vajont artificial basin. In 1963, a catastrophic landslide, estimated to be about 270 million cubic meters, fell into the Vajont reservoir, generating a tsunami that caused about 2000 casualties [74].

The experimental campaign of Padua aimed to evaluate the effects of both the material type and the landslide falling time on the maximum run-up of the generated wave over the opposite side of the valley. The landslide velocity was imposed through a rigid plate pulled by an engine through a steel cable. Several tests were carried out considering two types of rounded gravel (3–4 mm and 6–8 mm), crushed stone, and squared tile. For each of these material types, different values of the run-out velocity were tested by varying the plate stroke (in the range 0.5–0.8 m) and its velocity, thus resulting in a landslide falling time ranging approximately between 15 s and 500 s at full scale.

The frames in Figure 5 show the simulation results obtained with SPHERA v.9.0.0 [75]; some representative instants were selected during the acceleration phase where the vertical rigid plate pushes the landslide toward the basin, and the subsequent run-out phase of the generated wave climbs the opposite slope. Fairly good results were obtained in terms of the maximum height reached by the rising wave front. The left-hand panels show the velocity field evolution. The right-hand panels show the density field; it can be noticed that the lower landslide layer (light-brown) has a higher density because of the saturation. The effect of pore water allows for the landslide dynamics to be obtained more close to the experimental results because at *t* =1.35 s, its front comes into contact with the opposite side of the basin. No tuning of the physical parameters of the model was necessary. The required accuracy was controlled by the proper adoption of the limiting viscosity value, μ0, reaching a reasonable compromise with consumed computational time. This test case is provided as tutorial number 35 in the documentation of SPHERA v.9.0.0 that is freely available in [76].

**Figure 5.** Frames of the simulation of the Padua experiment. (**a**) Velocity field; (**b**) Density field.

Calculation of the flow-impact induced forces on the submerged structure may be required in order to estimate landslide and water-related hazards. Using the standard WCSPH approach usually leads to high-frequency numerical noise in the pressure field [29] that may taint the calculation of the load time history acting on the solid structure. To overcome this problem and obtain the correct prediction of the hydrodynamic load for structural stability assessment, several strategies have been successfully introduced [77].

#### **3. Flooding in Complex Topography with the Transport of Sediments**

This section provides synthetic discussions on the following topics: state-of-the-art of CFD mesh-based codes for flood propagation on complex topography with sediment transport; advantages of SPH modeling for flood propagation on complex topography; state-of-the-art of SPH codes for sediment transport; example of a SPH code for flood propagation on complex topography with sediment transport; and pre-processing and post-processing tools for floods on complex topographies (i.e., real topographic surfaces or their scale models). The most popular codes used to represent erosional floods (i.e., flood propagation over granular beds) rely on mesh-based numerical methods and the shallow water equations (SWE). These are briefly recalled in the following. The work in [78] presents a finite element code for bed-load transport. A major novelty is the 3D mathematical formulation. However, the code validations only refer to 1D configurations in 2D domains where analytical solutions are available. The work in [79] introduced a Godunov-type 1D SWE model for 2D erosional dam-break floods. The work in [80] presented a 2D SWE-FVM (finite volume method) model for 3D erosional dam-break floods, which uses Exner's equation for the bed top evolution, a 1D heuristic formula for the bed-load transport rate, and a 1D spatial reconstruction scheme based on Riemann solvers suitable for multi-phase flows. The work in [81] presented a 2D SWE model for 3D erosional dam-break floods; validations refer to 1D and 2D bed-load transport configurations, whereas a 3D demonstrative configuration (still represented with a 2D code) was reported on a quake-induced erosional dam-break flood (Tangjiashan Quake Lake). The work in [82] applied a 2D SWE-FVM code to simulate an erosional

flood with bed-load transport in the Yellow River. The last two examples might represent typical applications of 2D SWE mesh-based models to 3D erosional floods on complex topography with 1D schemes for sediment transport.

The above mesh-based models share the following drawbacks: 2D modeling; no consistency with the KTGF; ad-hoc tuning procedure for the viscosity of the 2-phase mixture (water and granular material); and shallow-water approximation (e.g., hydrostatic pressure profiles, velocity is uniform along the vertical).

With respect to the above state-of-the-art mesh-based models, the meshless SPH numerical method introduces several advantages, which are synthesized in the following. SPH allows the detailed 3D fluid dynamics fields within the urban canopy (urban fabric) to be simulated including fluid-structure interactions, whilst the mesh-based 2D porous models do not provide a direct modeling of the flood-building interactions. The SPH method can simulate the 3D transport of solid structures and fluid interactions with other mobile and fixed structures and can also represent 3D bed-load transport and its impact on mobile and fixed solid structures. The SPH method provides a direct estimation of the position of the free surface and the fluid and phase interfaces due to its Lagrangian nature. The direct representation of Lagrangian derivatives is also responsible for the absence of the advective non-linear terms arising in the Eulerian formulation of the balance equations. No computational mesh generation is requested, thus saving person-months, software licenses, and computational resources. Nonetheless, several drawbacks have been reported. SPH is slightly more time consuming then mesh-based methods (fixed with the same spatial resolution and accuracy) and has a narrower, but peculiar range of application fields, whose number is nonetheless elevated and relevantly involve floods. Furthermore, a SPH code can usually cover a wide range of spatial resolutions at different accuracy levels, but maintain stable algorithms. This allows the same code to be used for both preliminary analyses at a coarse spatial resolution and accurate simulations (at fine spatial resolution).

Although several SPH studies are available for 2D and 3D applications for granular flows, only a few have been dedicated to floods with sediment transport (over a simplified topography). Among these studies, [7] introduced a 2D erosion criterion to represent the sediment removal from water bodies by means of discharge channels (i.e., flushing procedures). Regardless of the application, the SPH models for granular flows were mostly restricted to 2D codes, featured by either 2-phase models or ad-hoc tuning for mixture viscosity.

Recently, a numerical mixture model for dense granular flows was presented in [8] to simulate the sediment dynamics phenomena, which typically involve the failure of earth-filled dams and dykes, bed-load transport, and fast landslides over complex topography. This model is discussed in the following. This was integrated into the FOSS SPH code SPHERA (RSE SpA) [76]. This mixture model permits the simulation of the above phenomena by solving a system of balance equations, coherent with the theoretical state-of-the-art frame represented by the KTG [83], under the conditions of the "packing limit". The model is based on mixture parameters (velocity, density, viscosity, and pressure) and phase variables (e.g., mean effective stress, frictional viscosity, and liquid phase pressure). The viscous parameters (mixture viscosity, frictional viscosity, and liquid viscosity) do not need any calibration/tuning (Section 2.1) [8]. Filtration is partly and implicitly represented; despite the absence of an explicit filtration scheme, a Lagrangian sub-scheme for saturation conditions was based on the hypotheses of stratified flows and local 1D filtration flows parallel to the local seepage [8]. A separated treatment involves the mixture particles under the elastic–plastic strain regime: they are held fixed as their velocities are negligible for applications such as bed-load transport and fast landslides [8]. The mixture model allows a high number of fluids to be simultaneously represented in the same domain, provided that they are either liquids or granular materials (both fully saturated or dry).

The model above was validated through laboratory experiments [8] and applied to a 3D erosional dam-break flood on complex topography (Figure 6). This demonstrative test case showed the applicability of the SPH method in simulating a flood on complex topography with bed-load transport. The 3D erosional dam break was triggered by an instantaneous and almost complete failure of a

gravity dam, whose structure was not simulated. The water flow impacts, thrusts, or erodes, and then transports a portion of the downstream mobile bed, which is composed of a bed of granular material. Its original sedimentation was related to the presence of a weir, whose structure was then removed before the dam building. The 3D fields of the velocity vector and pressure were computed and the 2D fields of the maximum (over time) water depth and specific flow rate (i.e., flow rate per width unit) were elaborated. These quantities, requested for risk analyses, were derived from the particle and the topography heights, and the magnitude of the depth-averaged velocity vectors. The time series of the water depth, the fluid volumes (cumulated in selected sub-domains), and the flow rate were assessed at specific monitoring sections. This demonstrative test, reported in [8] and available as tutorial no. 18 of SPHERA v.9.0.0 [76], shows the potential of the SPH method in simulating a full-scale 3D flood on complex topography with sediment transport. Although validations were reported in terms of comparisons with laboratory datasets [8], these only refer to simplified topographies. A full scale validation still has to be investigated due to the lack of available measures for complex topographies with sediment transport.

**Figure 6.** Demonstrative test case for the meshless SPH method applied to a 3D erosional dam-break flood on complex topography [8]. (**a**) 3D view. (**b**) Top view. Digital elevation model (grey, with a black vertex triangulation), liquid SPH particles (blue), mixture SPH particles (brown).

CFD codes for flood propagation on complex topography need a suitable numerical chain. A non-exhaustive list of pre-processing and post-processing free tools for floods on complex topography is discussed hereafter.

The dataset SRTM3 (USGS) represents the most accurate open data DEM ("digital elevation model", not to be confused with the "discrete element method" from Section 2.2) archive with a spatial resolution of 1" (spatial resolution length scale approximately equal to 31 m) and an almost global cover. The SRTM3 files are available in the ".tif" format. The numerical tool GDAL [84] can be used to convert the DEM ".tif" file format of SRTM3 in the alternative format ".dem". DEM2xyz [85] can read the DEM file (".dem" format), convert the geographic coordinates in Cartesian coordinates over a regular grid, and write the resulting DEM on an output file (".xyz" format), possibly coarsening the spatial resolution and reconstructing bathymetry where it is not available. Paraview (Kitware) [86] can read the ".xyz" output file of DEM2xyz and elaborate a 2D Delaunay grid starting from the DEM vertices. Paraview also allows cutting the numerical domain, drawing particular regions of interest (e.g., water bodies), engineering works (e.g., dams), and monitoring elements (e.g., points, lines, surfaces). The above information, derived from Paraview, can be provided to DEM2xyz, which can be executed again to add the new elements designed with Paraview.

The 3D fluid dynamics output files of a CFD code for floods can be further elaborated by means of Grid Interpolator [87], which reads a 3D field of values from an input grid and interpolates them on an output grid with a different spatial resolution. At this point, Paraview can be used to visualize the 2D and 3D fluid dynamics fields and return the associated image files.

All the pre-processing and post-processing items above are FOSS, but SRTM3, which is "open-data" (dataset available upon public and free access) [88] can be used within the same numerical chain and can be replaced by other free or proprietary items depending on the user resources and targets.

#### **4. High Performance Computing Solutions for Complex Hydraulic Engineering Problems**

Numerical modeling is becoming more useful and practical thanks to the capability of the current computer hardware. Years ago, significant simplifications of the problem needed to be undertaken in order to make the numerical simulation feasible due to past hardware limitations. Nowadays, thanks to the continuous hardware improvement and the use of HPC techniques that allow for the advantages of the enormous calculation power of the current hardware to be taken, it is possible to now simulate complex problems with fewer simplifications, and problems that were impossible to be simulated due to their spatial and/or temporal scale can now be carried out. Therefore, several codes have been developed that use these HPC techniques to simulate complex hydraulic engineering problems of water-related natural hazards in reasonable execution times. Some examples are SPHERA for sediment transport [8] and landslides [9]; DualSPHysics to model the scouring of two-phase liquid-sediments flows [6,53]; and GPUSPH for simulating lava flows [74]. As explained in [89], due to the limitation of sequential computing and the high cost of increasing performance in sequential architectures, the era of sequential computing was replaced by the era of parallel computing. Currently, hardware performance is fundamentally increased by expanding the number of processing units, instead of increasing the power of a single processing unit. Therefore, it is mandatory to use parallel programming techniques to distribute the workload among the available processing units and synchronize their execution.

OpenMP (Open Multi-Processing) is the most common parallel programming technique for multi-core systems with shared memory because it does not involve major code changes. MPI (message passing interface) is a standard for distributed memory systems that allows for the work to be divided among the multiple calculation nodes connected to each other. On the other hand, and more recently, GPGPU (general-purpose computing on graphics processing units) has emerged as an alternative to the use of traditional CPUs (central processing units). This technique allows the code to be executed on the GPU (graphics processing units). GPUs initially developed for rendering graphics have greatly increased their parallel computing power over the last decade due to the increasing computing demands of the gaming market, and more recently due to the significant growth of markets related to AI (artificial intelligence), DL (deep learning), and DNN (deep neural networks). They are in fact now several times faster than CPUs for problems that demand a high computational cost whilst offering a high degree of parallelism, to which SPH falls under [90]. At the same time, the emergence of programming languages such as CUDA (compute unified device architecture) and OpenCL makes the development of GPU codes easier. These two reasons have encouraged many researchers to implement their models for GPU, thus starting the era of GPU computing [91].

Particle methods such as smoothed particle hydrodynamics (SPH) [92] or moving particle semi-implicit (MPS) [93] have a very high computational cost. The high computational cost of particle methods is due to the fact that each particle has to interact with a large number of neighboring particles, which, in addition, may change every time step. The way that particles are organized in memory and the neighbor search (NS) algorithm is key to maximizing performance in parallel architectures. The most common algorithms for NS are the Verlet list and Cell-linked list. The work in [94] compared both algorithms applied to SPH and, in [95], variants of these algorithms were evaluated to conclude that the Verlet list may be faster than Cell-linked list, but the memory consumption is prohibitive for a high number of particles, so the Cell-linked list was found to be more efficient. The same conclusion was confirmed in [96] for GPU simulations. Quad-tree partitions combined with the Morton space filling curve (SFC) was recommended when using variable resolution in multi-core architectures [97] and GPU devices [98]. The work in [99] proposed another algorithm based on hierarchical cell decomposition

for variable resolution in distributed memory architectures using MPI. A very fast NS, where the speed is the main priority, was implemented in [100] for real-time animations, but this approach is invalid for real-physic simulations since some interactions of neighboring particles are missing.

The high computational cost of particle methods can be countered using the programming techniques discussed above. Thus, different implementations of the SPH method using OpenMP can be found [97,101,102], although the same algorithms can be applied to other particle methods. The work of [101] included a performance analysis of SPH in multi-core CPUs and studied how to avoid bottlenecks. The work in [102] implemented an optimized SPH for parallel machines with shared memory and compared the effective computing performance on multi-core CPUs and Xeon Phi MIC (many integrated core) coprocessors with 61 cores. The work in [97] evaluated the efficiency of a neighbor search algorithm based on a quad-tree partitioning on 32 processors. The SPH method has also been implemented for distributed memory systems using MPI for WCSPH [99,103–106], and MPI for ISPH [107,108]. The execution on distributed systems requires the division of the simulation domain into multiple subdomains. Each subdomain is executed by an independent process, but all processes must be synchronized for every calculation step. Therefore, load balancing between the different processes is essential to minimize the time that each process waits for the rest of the processes with a consequent loss of performance. In addition, this balancing needs to be dynamic where the workload can shift from one subdomain to another during simulation in Lagrangian methods. On the other hand, the decomposition of the domain into distributed memory systems requires the transfer of data between the nodes that process each domain, especially in the case of GPUs where transfers via PCI Express on the node itself are also required. The overhead of these memory transfers can be significant, so it is a major problem in the dynamic decomposition of the domain. In [103], the dynamic load balancing was based on the METIS package. The work in [104] implemented the code JOSEPHINE using the programming language Fortran 90 and OpenMPI, obtaining a good efficiency on 16 processors. The implementation of [99] used a hierarchical cell decomposition to allow variable resolution. The work in [105] obtained a high efficiency with 32,768 CPU cores using a decomposition algorithm based on orthogonal recursive bisection [109]. Most recently, [109] implemented dynamic domain decomposition between Voronoi subdomains and achieved a good speedup using 1024 processes. The work in [107] presented an ISPH solver where a simple spatial domain decomposition in cells was combined with the Hilbert space filling curve (SFC) and load balancing was only applied when particle imbalances were detected. This implementation reported an efficiency of 63%, simulating a still water case with 66.5 million particles on 512 processors compared to the execution time using 32 processors. The ISPH implementation of [108] achieved an efficiency close to 80% on 1536 MPI nodes, simulating a dam break with 100 million particles compared to the runtime using one node with 24 cores. Nowadays, powerful GPUs can be used instead of expensive parallel machines or clusters based on CPU processors to simulate several millions of particles. The first full GPU implementation of the SPH method was reported in [110] using shader programming. Other implementations on GPU are described in the works [52,111–115]. The work in [111] presented the SPH code, named DualSPHysics, implemented for GPU or CPU executions where a Cell-linked list [116] is used for a neighbor search. This work showed a comparison between several GPU models of different generations and several CPU models, achieving an improvement of almost two orders of magnitude. Hérault et al. [112] described another SPH solver for GPU named GPUSPH using a Verlet list. This solver allowed them to simulate lava flows on a complex topography with accuracy and efficiency, obtaining speedups of up to two orders of magnitude with respect to an equivalent CPU as claimed by the authors of [112]. The work in [52] presented an improved DualSPHysics code where the CPU and GPU optimization strategies described in [96] were applied to obtain a fair comparison between GPU and multi-core CPU executions. A multi-phase SPH model for liquid–gas simulations on GPUs based on DualSPHysics was shown in [113]. The use of GPUs is mandatory in this multi-phase approach since both fluids, water, and the surrounding air volume must be discretized in SPH particles. Therefore, the total number of particles and the execution time are usually increased several times. The work in [114] described a 2D SPH code optimized for long hydraulic simulations using the Hilbert space filling curve to improve the data locality and increase the performance. The work in [115] presented the AQUAgpusph code using the programming language OpenCL. The use of OpenCL allows the use of NVIDIA or AMD GPUs and the same code can be executed on multi-core CPUs. GPU implementations of ISPH can be found in [117–119]. The work in [117] presented a solver written in OpenCL that could be executed on a CPU or GPU and the performance of both architectures was compared. The work in [118] showed an ISPH code to perform realistic fluid simulations in real time, so some simplifications were applied to increase the speed of execution. The work in [119] implemented an optimized GPU solver starting from [111] and achieved a speedup close to 5 against CPU simulations with 16 threads. GPU implementations of MPS can also be found in [120] for 2D and in [121] for 3D free surface simulations.

Several examples of solvers based on shallow water equations on GPU are in [98,122–126]. These allow for large-scale flood simulations to be performed using only one GPU. The work in [98] presented a SPH solver for shallow water with adaptive resolution using a quad-tree algorithm for neighbor search. The work in [122] presented a 2D flood model that allowed for the simulation of long duration floods in a few minutes. The work in [123] simulated more than one hour of the Malpasset dam failure case in less than 30 s and proved that the use of single precision was enough for that problem. The work in [124] presented a finite volume solver with explicit discretization and an efficient algorithm to deactivate dry zones improving performance. The work in [125] introduced an optimized version of the IBER solver [127] for CPU and GPU simulations. That new version of IBER solved 24 h of an extreme flash flood in less than 10 min. The work in [126] presented a shallow water model based on the finite volume approach implemented for GPU with the OpenACC language [128]. OpenACC allows GPU parallelization to be implemented automatically using compiler directives like OpenMP for CPUs. Some examples of flood simulation, possibly thanks to the computing power of a single GPU, can be found in [18,129,130]. The work in [129] used more than 1.5 million particles to simulate a complicated city layout including underground spaces with the GPUSPH model [112]. The work in [18] simulated runoff on a real terrain generated from photogrammetry information obtained by an UAV (unmanned aerial vehicle). This simulation with 7.5 million particles was performed using the DualSPHysics code [111] and took 138 h for 15 physical minutes. The work in [130] showed a large-scale urban flood performed by the GPU model implemented in [98]. This simulation with 230,000 particles took 135 h to simulate five physical hours.

One GPU can perform simulations with a high particle number in affordable execution times as demonstrated by the works mentioned before. The same results on CPU systems would require expensive CPU clusters. However, a fair comparison between GPU and CPU performance is not straightforward, since optimized solvers for both architectures and hardware from the same period should be used to avoid unreal speedups [116,131]. The work in [132] disproved the speedups of 100x or 1000x as shown in some GPU–CPU comparisons.

Although today's GPUs provide high computation power, the simulation of real cases implies huge domains with a high resolution, which implies simulating tens or hundreds of million particles and these simulations are not viable in a single GPU due to memory limitations and prohibitive execution times. The solution is to use multi-GPU machines with shared memory or GPU clusters. Examples of multi-GPU SPH solvers can be found in [133–136]. The work in [133] explored the use of MPI in the DualSPHysics code to perform multi-GPU simulation on GPU clusters. This approach allowed for the simulation of 32 million particles on four GPUs, achieving an efficiency of 80% using weak scaling by simulating eight million particles per GPU. The work in [134] described the implementation of an optimized multi-GPU version of DualSPHysics using an MPI that was able to simulate 1024 million particles on 128 GPUs, achieving almost 100% efficiency using weak scaling by simulating eight million particles per GPU. The work in [135] extended the GPUSPH solver to multi-GPU machines with shared memory using threads and domain decomposition in two dimensions. This implementation is limited by the number of GPUs hosted by one machine or computation node, typically six or eight GPUs. The work in [136] presented a multi-GPU implementation with 3D domain decomposition that achieved 89% efficiency using strong scaling by simulating more than 30 million particles on eight GPUs using MPI.

#### **5. Conclusions and Future Perspectives**

This paper collected some recent works showing the application of CFD techniques for modeling problems of practical and theoretical interest involving complex multiphase flows relevant for the analysis and mitigation of water-related natural hazards. The paper focused on meshless techniques for the numerical modeling of fast landslides, tsunami wave, flooding in complex geometry and sediment scouring; few relevant examples have also been mentioned concerning traditional grid-based methods applied to the analysis of environmental risks related to flooding in complex topography. The peculiar features of the examined works in terms of mathematical modeling and numerical implementation have been illustrated and outline the principal differences.

Different approaches are commonly used for modeling the non-Newtonian dynamics of dense granular flows. Some of them consider the sediment as a single-phase moving mass; other approaches model the sediment as a two-phase mixture where the voids of the solid matrix are filled by some liquid phase in order to account for the pore pressure effects. In all of the works examined where the SPH method was applied to the solid, the motion of the granular phase was treated as a pseudo-fluid once the sediment particles had been mobilized. Furthermore, the criterion adopted for determining the onset of sediment motion is, in general, not uniform. In some of the considered works, a numerical threshold was used for this purpose and the involved mechanical parameters, despite their physical meaning, require proper tuning so that the model results can suitably fit the experimental data. This approach is used both in the analysis of bed sediment erosion and landslide run-out where the definition of a triggering threshold that establishes the critical condition for the onset of sediment motion is required. In these cases, the need of tuning parameters, which assume a numerical role rather than a physical meaning, may limit the applicability of the model to those situations of practical interest where calibration data are available. An alternative approach to overcome this limitation is provided by the kinetic theory of dense granular flows (KTGF) that is put into SPH formalism and has been proven to be able to simulate dense granular flows (the so-called "packing limit" of KTGF) and fast landslides with a suitable degree of accuracy.

Proper numerical modeling of the erosion processes that are related to the motion of loose sediment particles suspended in the water flow would require that the size of the numerical particles should be of the same order of the size of real sediment grains; but this task is not allowed in practical problems where the geometrical complexity may induce excessive computational cost. However, reliable results may be obtained by simulating numerical particles as a combination of clear water and turbid water particles that mimic a suspended load.

Concerning the interaction between the phases (i.e., water and sediment), different approaches have been followed. In some cases, the governing equations of motion are solved simultaneously for both phases, thus obtaining coupled dynamics. In other cases, the governing equations are solved separately for each phase and the interface conditions should be enforced to assure kinematic and dynamic continuity between phases through an iterative procedure that can, however, affect the computational time. While some models were examined in terms of the mathematical features, these were not the same methods that were examined in terms of numerical implementation.

The above-mentioned relevant aspects may exert significant influence both in the numerical representation of the multiphase flow and in the reliability of the model results. Therefore, they should be carefully evaluated in the analysis of water-related natural hazards in order to obtain a reliable representation of the investigated problem. Given the complexity (especially geometrical) that usually characterizes the problems of practical interest, these numerical models could be used to support risk analysis and mitigation if appropriate programming techniques and modern architectures for scientific computation are used to obtain fast-running computer codes. This goal is not simple to

accomplish, even with the adoption of HPC techniques and parallel computing. The simulation of real cases with large domains and high resolution will probably become more affordable with the use of GPU clusters. However, this task needs better implementation of the identified useful formulations of SPH for multi-GPU multi-node. This requires overcoming the problems associated with domain decomposition and memory transfer across nodes, which are particularly difficult for the fully-coupled two-phase formulations.

**Author Contributions:** Conceptualization, R.A. and S.M.; Design and enhanced structure of the manuscript, R.A. and S.M.; Writing Sections 1, 2 and 5–original draft preparation, S.M.; Writing Section 2.1–original draft preparation, D.W.; Writing Section 3–original draft preparation, A.A.; Writing Section 4–original draft preparation, J.M.D.; Writing—review and editing, S.M., D.W., J.M.D., S.L., A.A., and R.A.

**Funding:** This work was partially financed by Xunta de Galicia (Spain) under project ED431C 2017/64 "Programa de Consolidación e Estructuración de Unidades de Investigación Competitivas (Grupos de Referencia Competitiva)" co-funded by European Regional Development Fund (FEDER), and by the Xunta de Galicia postdoctoral grant ED481B-2018/020. The work was also funded by the Ministry of Economy and Competitiveness of the Government of Spain under project "WELCOME ENE2016-75074-C2-1-R" and by the EU under the ERDF (European Regional Development Fund) through the Interreg project "MarRISK" (0262\_MARRISK\_1\_E). The contribution of the RSE author (Section 3) was financed by the Research Fund for the Italian Electrical System (for "Ricerca di Sistema -RdS-") in compliance with the Decree of Minister of Economic Development 16 April 2018.

**Acknowledgments:** We acknowledge the CINECA award under the ISCRA initiative for the availability of high performance computing resources and support. HPC simulations on SPHERA refer to the following HPC research project: HPCNHLW1—High Performance Computing for the SPH Analysis of Natural Hazard related to Landslide and Water interaction (Italian National HPC Research Project); instrumental funding based on competitive calls (ISCRA-C project at CINECA, Italy); 2019; S. Manenti (Principal Investigator) et al.; 360,000 core-hours. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. HPC simulations on SPHERA (Section 3) refer to the following HPC research project: HPCEFM17—High Performance Computing for Environmental Fluid Mechanics 2017 (Italian National HPC Research Project); instrumental funding based on competitive calls (ISCRA-C project at CINECA, Italy); 2016–2017; 200,000 core-hours; Amicarelli A. (Principal Investigator) et al.

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

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **A New Parallel Framework of SPH-SWE for Dam Break Simulation Based on OpenMP**

#### **Yushuai Wu 1, Lirong Tian 1, Matteo Rubinato 2, Shenglong Gu 1,3,\*, Teng Yu 1, Zhongliang Xu 4, Peng Cao 5,6, Xuhao Wang 6,7 and Qinxia Zhao <sup>1</sup>**


Received: 7 April 2020; Accepted: 8 May 2020; Published: 14 May 2020

**Abstract:** Due to its Lagrangian nature, Smoothed Particle Hydrodynamics (SPH) has been used to solve a variety of fluid-dynamic processes with highly nonlinear deformation such as debris flows, wave breaking and impact, multi-phase mixing processes, jet impact, flooding and tsunami inundation, and fluid–structure interactions. In this study, the SPH method is applied to solve the two-dimensional Shallow Water Equations (SWEs), and the solution proposed was validated against two open-source case studies of a 2-D dry-bed dam break with particle splitting and a 2-D dam break with a rectangular obstacle downstream. In addition to the improvement and optimization of the existing algorithm, the CPU-OpenMP parallel computing was also implemented, and it was proven that the CPU-OpenMP parallel computing enhanced the performance for solving the SPH-SWE model, after testing it against three large sets of particles involved in the computational process. The free surface and velocities of the experimental flows were simulated accurately by the numerical model proposed, showing the ability of the SPH model to predict the behavior of debris flows induced by dam-breaks. This validation of the model is crucial to confirm its use in predicting landslides' behavior in field case studies so that it will be possible to reduce the damage that they cause. All the changes made in the SPH-SWEs method are made open-source in this paper so that more researchers can benefit from the results of this research and understand the characteristics and advantages of the solution proposed.

**Keywords:** dam break; SWE; SPH; openMP; numerical modelling; computational time

#### **1. Introduction**

The Smooth Particle Hydrodynamics (SPH) is a meshless method [1] very commonly used nowadays [2–12]. Gingold and Monaghan [13] were the first to propose this method to solve astrophysical simulations, using statistical techniques to recover analytical expressions for the physical variables from a known distribution of fluid elements. The SPH method is typically used for solving the equations of hydrodynamics in which Lagrangian discretized mass elements are followed [14]. Compared with the limitations of the Eulerian grid method [15–17], the SPH method has unique advantages in dealing with free water surface and moving boundary conditions [18–21]. In fact, the SPH method strictly runs in accordance with the law of conservation of mass and can deal with free surface and moving boundary flexibly; hence, it is very suitable for simulating dam break flows [22]. Flooding due to dam break has potentially disastrous consequences, and multiple studies were conducted to numerically replicate the hydrodynamics of this phenomenon [23–28]. In most cases, the dam break flow occurs in a wide area and lasts for a long time. Therefore, the Shallow Water wave Equations (SWEs) have become the main application formulae for this specific dam break problem.

In 1999, Wang and Shen [29] applied the SPH method to SWEs for the first time. Dam break flows are unsteady open channel flows that can be described by the St. Venant equations, which are equations that can be used for flows with strong shocks [29]. The study conducted by Wang and Shen [29] has demonstrated that the method developed based on the SPH is capable of providing accurate simulations for mixed flow regimes with strong shocks. In 2005, Ata and Soulaimani [30] tried to reduce the difficulties that were associated with the treatment of the solid boundary conditions, especially with irregular boundaries. Ata and Soulaimani [30] derived a new artificial viscosity term by using an analogy with an approximate Riemann solver, and the several numerical tests conducted have confirmed that the stabilization proposed provides more accurate results than the standard artificial viscosity introduced by Monaghan [22]. However, Ata and Soulaimani [30] found that it was difficult to implement the Dirichlet boundary conditions for bounded domains. Different techniques such as symmetrization and ghost particles were implemented; nevertheless, results obtained for irregular boundaries or in presence of shocks were not satisfactory, and there was a need to make SPH method more competitive with standard approaches [30].

A new method with good stability was then needed for the SPH numerical simulation of shallow water equations to deal with dam break flows, flood waters, debris flows, avalanches, and tidal waves.

De Leffe et al. [31] proposed an improved calculation method based on a two-dimensional SPH solid wall boundary condition. By introducing a periodic redistribution of the particles and using a kernel function with variable smoothing length, this modification was tested and validated against dam break flows on a flat dry bottom in 1D and 2D. Comparisons conducted against literature results [30,32,33] have confirmed how this new approach is robust and able to simulate complex hydrodynamic situations.

However, to date, only a few studies have investigated the efficiency of solving the SPH-SWEs model. Vacondio et al. [34–37] have developed the serial code to solve the SWEs by using the SPH method and have made an open source version called SWE-SPHysics, which has been optimized and adapted based on the hydrodynamics investigated by other researchers [38–52]. Despite continuous progress, there is still a limitation related to the computational efficiency when the number of particles to simulate is very large, and this aspect still needs to be improved.

Xia and Liang [53] explored the Graphic Processing Units (GPUs) to accelerate an SPH-SWE model for wider applications such as dam breaks. Xia and Liang [53] demonstrated that the performance of the new GPU accelerated SPH-SWE model can significantly improve the calculation efficiency and verified that the quadtree neighbor searching method may reduce redundant computation when searching neighbor particles [53]. Furthermore, Liang et al. [54] developed a shock-capturing hydrodynamic model to simulate the complex rainfall-runoff and the induced flooding process in a catchment in England, Haltwhistle Burn, of 42 km2, and implemented it on GPUs for high-performance parallel computing.

GPU has been enhanced with a Fortran programming language capability employing CUDA (Compute Unified Device Architecture), known as CUDA Fortran [55]. Although the GPU parallel computing performance is strong, the GPU price is relatively expensive and support for the CUDA FORTRAN language compiler is limited. Furthermore, if CUDA FORTRAN language is compared with other parallels (OpenMP and MPI), the program design is also more complex [56–58].

Over the years, with the development of parallel computing, the development of OpenMP and MPI in parallel methods has matured, and MPI is widely used in the field of engineering computing [59–63]. However, in this multi-machine cluster environment, memory is not shared. For global shared data operations, data must be transferred by the communication between machines [64,65].

OpenMP is based on the shared storage mode of multi-core processors, and it is commonly used in parallel processing of single workstations. Although it is limited by the processing capacity and memory capacity of a single node, it can simplify the past multi-core computing to the present multi-core computers. The program design is relatively simple, and it can secure the advantages of economic and programming optimizations [66].

Therefore, this paper adopts the parallel computing method of CPU-OpenMP that is applied to a single machine and a multi-core to calculate the new SPH-SWEs framework for the parallel computing of the test case of *2-D dry-bed dam break with particle splitting* [67,68] and a *2-D debris flow with a rectangular obstacle downstream the dam* [2]. The accuracy of the new SPH-SWEs framework was then verified by comparing the serial algorithm, and the advantages of CPU-OpenMP parallel computing were analyzed.

The paper is organized as follows: Section 2 describes the methodology adopted presenting the theoretical derivation of the numerical model applied and the governing equations. Section 3 explains the setup of the SPH-SWEs method used. Section 4 provides the results of the application tested with a discussion of the results obtained. Lastly, Section 5 produces a brief summary and concluding remarks of the whole study.

#### **2. Methodology**

Ata et al., [30] and Paz and Bonet [32] have initiated the idea of SPH-SWE model solution, which is described in Section 2.1 with all the governing equations.

#### *2.1. Governing Equations*

By ignoring the Coriolis effect and the fluid viscosity, SWEs can be written in the Lagrangian form as follows:

$$\begin{array}{ll}\frac{d\mathbf{d}}{dt} = -d\nabla \cdot \mathbf{v} \\ \frac{d\mathbf{v}}{dt} = -\mathbf{g}\nabla d + \mathbf{g}(\nabla b + \mathbf{S}\_{\mathbf{f}}) \end{array} \tag{1}$$

*d* represents the water depth, *g* the acceleration of gravity, *v* is velocity, *b* represents the riverbed elevation, and *Sf* represents the riverbed friction. In the SWEs, the area density is defined as:

$$
\rho = \rho\_w d \tag{2}
$$

ρ represents the density, and ρ<sup>w</sup> represents the density of water.

#### *2.2. Water Depth Solutions*

According to the SPH idea, the area density (i.e., water depth) of particles is solved as shown below in the implicit function:

$$\begin{array}{rcl}\rho\_i &=& \sum\_j m\_j \mathbf{W}\_i(\mathbf{x}\_i - \mathbf{x}\_j, h\_i) \\ h\_i &=& h\_0 \binom{\rho\_0}{\rho\_i}^{1/d\_m} \end{array} \tag{3}$$

*xi*/*xj* represents the particle coordinates; *mj* represents the particle mass of *j*; *hi* and ρ*<sup>i</sup>* represent the smooth length and the area density of the particle *i*; *h*<sup>0</sup> and ρ<sup>0</sup> represent the initial values of the smooth length and the area density, respectively; *d*<sup>m</sup> represents the latitude (1 represents one dimension, 2 represents two dimensions) and *W* represents the kernel function.

#### *2.3. Speed Solution*

According to the Lagrangian equation of motion [69], *ai* is the acceleration of the particle *i*, and the solution formula of each particle can be obtained as follows:

$$\mathbf{a}\_{l} = \frac{\mathbf{g} + \mathbf{v}\_{l} \cdot \mathbf{k}\_{l} \mathbf{v}\_{l} + \mathbf{t}\_{l} \cdot \nabla b\_{i}}{1 + \nabla b\_{i} \cdot \nabla b\_{i}} \nabla b\_{i} - \mathbf{t}\_{l} + \mathbf{S}\_{f,i} \tag{4}$$

*ki* = ∇(∇*bi*) represents *b*(*x*) of the curvature tensor [70]; *ti* represents the acceleration caused by the internal force; ∇*bi* represents the riverbed gradient of the particle *i*, and in order to deal with any complex terrain problem, the riverbed gradient can be modified as follows [71,72]:

$$
\nabla b\_{\bar{i}} = \sum\_{\bar{j}} b\_{\bar{j}} \widetilde{\nabla} \mathcal{W}\_{\bar{i}} (\mathbf{x}\_{\bar{i}} - \mathbf{x}\_{\bar{j}}, h\_{\bar{i}}) V\_{\bar{j}} \tag{5}
$$

∇*Wi* denotes the gradient of the modified kernel function, which is modified by the correction matrix *Li*, as shown below:

$$\widetilde{\nabla} \mathcal{W}\_{i} \{ \mathbf{x}\_{i} - \mathbf{x}\_{j}, h\_{i} \} = \ L\_{i} \nabla \mathcal{W}\_{i} \{ \mathbf{x}\_{i} - \mathbf{x}\_{j}, h\_{i} \} \mathcal{L}\_{i} = \left\lvert \sum\_{j} \nabla \mathcal{W}\_{i} \{ \mathbf{x}\_{i} - \mathbf{x}\_{j}, h\_{i} \} \times \{ \mathbf{x}\_{i} - \mathbf{x}\_{j} \} V\_{j} \right\rvert^{-1} \tag{6}$$

In order to reduce the numerical oscillation and ensure the stability of the calculation, one method is to increase the viscosity term as introduced below [73]:

$$\begin{cases} \mathbf{t}\_{i} = \sum\_{j} m\_{j} \frac{\mathcal{S}}{2\rho\_{w}} \Big[ (\frac{1}{\mathcal{P}\_{j}} + \pi\_{ij}) \nabla W\_{j} (\mathbf{x}\_{i} - \mathbf{x}\_{j}, \mathbf{h}\_{j}) - (\frac{1}{\mathcal{P}\_{i}} + \pi\_{ij}) \nabla W\_{i} (\mathbf{x}\_{j} - \mathbf{x}\_{i}, \mathbf{h}\_{i}) \Big] \\ \quad \boldsymbol{\beta}\_{i} = \begin{array}{c} -\frac{1}{\rho\_{i} d\_{m}} \sum\_{j} m\_{j} r\_{ij} \frac{d \mathcal{W}\_{ij}}{d r\_{ij}} \\ \end{array} \\ \boldsymbol{\pi}\_{ij} = \begin{array}{c} \frac{\mathsf{T}\_{i\_{i}} \mathbf{v}\_{ij} \mathbf{x}\_{ij}}{\rho\_{j} \sqrt{\left\| \mathbf{x}\_{j} \right\|^{2} + \zeta^{2}}} \end{array} \tag{7}$$

β represents the correction coefficient caused by variable smooth length; *rij* represents the particle spacing; π*ij* represents the numerical viscosity added to maintain stability. However, this method has the problem of numerical dissipation.

To reduce this issue, the interaction between two particles was treated as a Riemann problem [74], as follows:

$$\begin{cases} \mathbf{t}\_{i} = \sum\_{j} m\_{j} p^{\*} \left[ \frac{1}{\rho\_{j}^{2} \beta\_{j}} \nabla \mathcal{W}\_{j} (\mathbf{x}\_{i} - \mathbf{x}\_{j}, h\_{j}) - \frac{1}{\rho\_{i}^{2} \beta\_{i}} \nabla \mathcal{W}\_{i} (\mathbf{x}\_{j} - \mathbf{x}\_{i}, h\_{i}) \right] \\ p^{\*} = 0.5 \boldsymbol{\mathcal{G}} \boldsymbol{\rho}\_{w} (\mathbf{d}^{\*})^{2} \\ d^{\*} = \frac{\frac{\mathcal{G} \boldsymbol{d} + \mathbf{g}, \mathbf{d}\_{r} + \mathbf{v}\_{l,n} - \mathbf{v}\_{r,n}}{\mathcal{G}\_{l} + \mathcal{G}\_{r}}}{\sqrt{0.5 \frac{\mathcal{G} \left( \boldsymbol{d}\_{l} + \mathbf{d}\_{k} \right)}{\mathcal{G}\_{l} \left( \boldsymbol{d}\_{l} \right)}}} \\ d\_{0} = \frac{1}{\mathcal{S}} \Big[ 0.5 (\boldsymbol{c}\_{l} + \mathbf{c}\_{r}) + 0.25 (\boldsymbol{v}\_{l,n} - \mathbf{v}\_{r,n}) \Big]^{2} \end{cases} \tag{8}$$

*dl* and *dr* represent the water depth on the left and right sides, respectively; *k* = *l* and *k* = *r* represent the left and right states, respectively; *<sup>d</sup>*<sup>0</sup> represents the initial estimated water depth; *<sup>c</sup>* = *gh* represents the shallow water wave velocity.

#### *2.4. Time Integration and Boundary Processing*

In order to update the particle velocity and displacement, the leap-frog time integration scheme [75] was used. By using this method, both time and space are of second-order accuracy, and the storage demand is relatively low, however the calculation efficiency is relatively high, as shown below:

$$\begin{array}{lcl}\nu\_{i}^{n+1/2} &=\nu\_{i}^{n-1/2} + \Delta t \mathfrak{a}\_{i}^{n} \\ \mathfrak{x}\_{i}^{n+1} &=\mathfrak{x}\_{i}^{n} + \Delta t \mathfrak{v}\_{i}^{n+1/2} \\ \nu\_{i}^{n+1} &=\nu\_{i}^{n+1/2} + \frac{1}{2}\Delta t \mathfrak{a}\_{i}^{n} \end{array} \tag{9}$$

Δ*t*represents the time step; where the time step must meet the Courant number condition [76] displayed as follows:

$$
\Delta t = \text{CFL} \min\_{i=1}^{N} \left( \frac{h\_i}{c\_i + \|\nu\_i\|} \right) \tag{10}
$$

To solve the boundary problem, this study adopts the Modified Virtual Boundary Particle (MVBP) method [36]. MVBP method is an improvement of the virtual boundary particle (VBP) method [77]. The virtual particles on the boundary will neither move with the fluid particles nor interact with them but generate the virtual particles similar to the mirror image through point symmetry. This method is easier and simpler in dealing with the complex boundary.

Compared with the VBP method, the MVBP method has two improvements: (1) When a virtual boundary particle is within the range of the kernel function of a fluid particle, two layers of newly generated virtual particles can be added (*Xk,*<sup>1</sup> = 2*Xv* − *Xi and Xk,*<sup>2</sup> = 4*Xv* − *Xi*). Among them, *Xk,*<sup>1</sup> and *Xk,*<sup>2</sup> represent the coordinates of the newly generated virtual particles, and *Xv* represents the virtual boundary particles; (2) When the internal angle of the boundary is less than or equal to 180◦, two newly generated virtual particles are added outside the corner. Compared with the single point of the VBP method, this improvement reduces the kernel truncation error.

In order to verify the effect of different kernel functions to simulate the dam break, this paper used the SPH-SWE open source code [34–37,68] and applied different kernel functions (*B-spline, super Gauss, quadratic spline, Gauss, quartic spline, quintic and Bell*) to simulate case 2 open source scenario [68]. These were the initial conditions considered for the dam-break: (i) simulation area was 2000 m long; (ii) the river bed elevation was 0 m; (iii) the initial fluid particle area was 1000 m long; (iv) the particle spacing was 10 m; (v) the initial water depth was 10 m, (vi) and the simulation duration was 50 s. The position, the water depth, and speed of the fluid particles were obtained as an output every 10 s. After verification, the advantages and disadvantages of different kernel functions and different numerical oscillation processing methods (see Equations (7)–(8)) identified from the results are consistent. At t = 50 s, the water depth dissipated after 1500 m, and the results of different kernel functions and different numerical oscillation processing methods can be visually reflected through the graphs in Figure 1; at the same time, it also illustrates the continuity of the numerical oscillation issue. Therefore, the data of t = 50 s (Figure 1 results) was selected for analysis in this study.

It can be noticed from Figure 1 that in the numerical simulation of dam break based on SPH-SWEs approach, all kernel functions are characterized by numerical oscillation except the option where the Bell kernel function is considered. The three kernel functions that provide a more accurate estimation of the water depth and the velocity are B-spline, quadratric spline, and quartic spline (above 85.7%). Figure 2a–d shows the absolute error between the calculated water depths and velocities using these three kernel functions vs. the analytical solution. Results displayed confirm the optimal performance of the B-spline kernel function in dam break simulations.

**Figure 1.** Simulation results of different kernel functions. (**a**,**c**) Results of water depths calculations. (**b**,**d**) Results of velocity calculations.

**Figure 2.** Absolute errors of different kernel functions and processing methods. (**a**–**c**) The absolute errors of water depth and velocity calculated by B-spline kernel function, quartic spline, and quadratic spline kernel function, respectively. (**d**) The absolute errors of water depth and velocity calculated by using the artificial viscosity method.

If the kernel function adopted is B-spline, the numerical oscillation is treated as adding the viscosity term (numerical viscosity, lax Friedrichs flux [78]). The results are shown in Figures 2d and 3. It is shown in Figure 3 that the Riemann solvers method has better characteristics in dealing with numerical oscillation problems when using the same kernel function. This confirms that it is possible to increase the viscosity by increasing πij in Equation (7), which helps to reduce the numerical oscillation; however, this can cause a decrease in the accuracy of the calculation results and an possible increase on computational time.

**Figure 3.** Simulation results of different processing methods.

According to Figure 2a,d and Figure 3, it can be found that the two-shocks Riemann solver is more advantageous in dealing with numerical oscillations when considering dam break cases. The absolute error of the solutions considered is significantly small. Table 1 displays the errors of the three methods adopted, and it can be seen that the standard deviations of the velocity and water depth (WD) of the three methods are small and, basically, of similar magnitude. However, the average relative errors of water depth and speed solved by the two-shocks Riemann solver seem to be the smallest, 0.92% and 5.86%, respectively.



• *WD* = *Water Depth*.

In order to verify the effect of the selected B-spline kernel function and the numerical oscillation processing method, a wet case was simulated (simulation range of 2000 m, initial water depth in the range of 0–1000 m is 10 m, and water depth in the range of 1000–2000 m is 5 m. The simulation time was 50 s, and the calculation results are generated every 10 s). The results are summarized in Table 2 and are shown in Figure 4.

According to the numbers displayed in Table 2, using the SPH-SWE model adopting the B-spline kernel function and the two-shocks Riemann solver method to solve the wet case, the simulation results are better, and their mean absolute errors to simulate velocity and water depth are within the 6%.

Based on these results, it was decided to select the B-spline kernel function and the two-shocks Riemann solver method to perform the two-dimensional dam-break numerical simulation to verify the computational efficiency of the new SPH-SWE model solution framework proposed.


**Table 2.** Error analysis for the wet case.

**Figure 4.** Simulation results of the SPH-SWE module. (**a**,**b**) water depth and velocity diagram at t = 30 s; (**c**,**d**) water depth and velocity diagram at t = 50 s.

#### **3. SPH-SWE Model Solution Framework**

When simulating dam break cases with a large number of particles involved, there is a challenge to be faced associated with low calculation efficiency. The code runs in serial steps and before each variable calculation, the mesh is divided into small grids (mesh size is *2H*, smooth length, in order to calculate the corresponding parameters of each particle), making the process more repetitive and requiring a lot of calculation time. Moreover, the code framework is complex and it is demanding to complete any modification. As the open source code solves the model with a large number of particles, it has the problem of low calculation efficiency and cannot even be calculated (the reason is that the array overflows). This paper proposes a new SPH-SWE model solution framework, which can dynamically allocate the storage space of particle information, solve the problems of repeated particle search and unsuccessful memory allocation of the array of stored particle information, and can quickly solve

the large-scale SPH-SWE model. Furthermore, the model framework proposed is simple, and any modification can be made easily to this algorithm.

For this study, Algorithm 1 was developed to solve the problem of data analysis and realize the CPU-OpenMP parallel computing.

**Algorithm 1.** Calculation framework of the SPH-SWEs model. This algorithm is needed to read the particles data (include fluid particles/virtual particles/open boundary particle/riverbed particles).

*Read parameters*

*Output initial data of the model*

*Mesh riverbed particles and calculate fluid particles and the net water depth* ! *Loop* 1 "

*Search particles* ! *Loop* 2 "

*do t* = 0 →*total\_number\_of\_timesteps*

*Step 1: Calculate the water depth of fluid particles* ! *Loop* 3 " .

$$\rho\_l = \sum\_{j} m\_j \mathbf{W}\_l (\mathbf{x}\_l - \mathbf{x}\_{j\*} h\_l) \quad h\_l = h\_l \left(\frac{\rho\_0}{\rho\_l}\right)^{1/d\_{lm}}$$

*Step 2: Calculate time water depth of fluid particles and the speed gradient* ! *Loop* 4 " .

$$\nabla n\_{l} = \sum\_{j} V\_{j} (n\_{l} - n\_{j}, h\_{l}) \nabla \mathcal{W}\_{l} (\mathbf{x}\_{l} - \mathbf{x}\_{j}, h\_{l}) \qquad (n = |d/\mu/\sigma|)$$

*Step 3: Calculate time increments.* <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *CFL <sup>N</sup>* min *i* = 1 *hi ci*+ν*i* 

*Step 4: Calculate accelerations of fluid particle,corrections of riverbed gradients,speeds, and displacements*

$$\begin{array}{rcl} \left\{ \begin{matrix} \boldsymbol{L}oop & \boldsymbol{5} \end{matrix} \right\} . \stackrel{\scriptstyle \boldsymbol{\overline{a}}\_{\boldsymbol{f}}}{\boldsymbol{a}} & = & \frac{\boldsymbol{\overline{\underline{\boldsymbol{\varepsilon}}}} + \boldsymbol{\overrightarrow{\overline{\boldsymbol{\nu}}}} \cdot \boldsymbol{\overrightarrow{\underline{\boldsymbol{k}}}} \cdot \boldsymbol{\overrightarrow{\underline{\boldsymbol{\nu}}}} + \boldsymbol{\overrightarrow{\underline{\boldsymbol{\Gamma}}}} \cdot \boldsymbol{\underline{\boldsymbol{\Gamma}}} \boldsymbol{\underline{\boldsymbol{\Gamma}}}}{\boldsymbol{1} + \boldsymbol{\underline{\boldsymbol{\Gamma}}} \boldsymbol{b} \cdot \boldsymbol{\underline{\boldsymbol{\Gamma}}} \boldsymbol{b}} \boldsymbol{\underline{\boldsymbol{\Gamma}}} b\_{\boldsymbol{\underline{\boldsymbol{\Gamma}}}} \end{matrix} \boldsymbol{\underline{\boldsymbol{\Gamma}}} b\_{\boldsymbol{\underline{\boldsymbol{\Gamma}}}} \stackrel{\scriptstyle \boldsymbol{\underline{\boldsymbol{\Gamma}}}}{\boldsymbol{1}} + \boldsymbol{\overrightarrow{\underline{\boldsymbol{\Gamma}}}} \boldsymbol{\underline{\boldsymbol{\Gamma}}} \boldsymbol{\underline{\boldsymbol{\Gamma}}} $$

*Step 5: Calculate displacements of open boundary particles. Step 6: Fluid particle division.* ν*<sup>k</sup>* = *c*<sup>ν</sup> *dN dk* <sup>ν</sup>*<sup>N</sup> <sup>c</sup>*<sup>ν</sup> <sup>=</sup> *AN* -*M <sup>k</sup>* <sup>=</sup> <sup>1</sup> *Ak Step 7: Calculate fluid particle of the riverbed* ! *Loop* 1 " .

$$\begin{cases} \;d\_l = \sum\_{j} d\_j \mathbf{W}\_l (\mathbf{x}\_l - \mathbf{x}\_{j\prime} \boldsymbol{h}\_l) V\_j & S\_l = \sum\_{j} \mathbf{W}\_l (\mathbf{x}\_l - \mathbf{x}\_{j\prime} \boldsymbol{h}\_l) V\_j \\\;\;\;\;\;\;\;correction: \; d\_l = \frac{d\_l}{S\_l} \end{cases}$$

*Step 8: Search partilces* ! *Loop* 2 " . *end do*

The calculation of each time step includes five parts, which can be summarized as follows:


In this paper, multiple two-dimensional arrays were used to store the information of fluid particles, riverbed particles, virtual particles, and open boundary particles. When calculating the parameters related to the above five parts, each particle can be calculated separately, and there is no dependency between particles.

Therefore, the CPU-OpenMP parallel computing can be realized (since each cycle is for all fluid particles, the calculation amount is known, so the schedule in the cycle configuration is set to Static mode) and the SPH-SWE model with a large number of particles can also be calculated.

#### *3.1. Fluid Particle Riverbed Calculation*

When calculating the riverbed particles in each time step, considering that the bed range and the smooth length of the riverbed particles are the same and there is no relationship between them, the calculation of the riverbed fluid particles was carried as follows.

Firstly, the grid division was completed, then the calculation of the riverbed fluid particles was conducted. It was verified that the calculation efficiency has no advantage, and additional array was needed to store riverbed particle information and increase data reading time. Please see below Algorithm 2 for the calculation framework adopted to achieve this task.


Because the calculated variables were two arrays, *h\_t and sum\_h\_t*, it was not difficult to realize OpenMP parallel, and multiple threads were set to calculate the riverbed fluid particles at the same time, following these formulae:

$$\begin{array}{ll} h\\_t(i) &= \sum\_j hb(j) \cdot \mathcal{W}\_i(\mathbf{x}\_i - \mathbf{x}\_j, hb(j)) \cdot Vol\_b b(j) \\ \displaystyle \sum\_j \mathcal{W}\_i(\mathbf{x}\_i - \mathbf{x}\_j, hb(j)) \cdot Vol\_b b(j) \end{array} \tag{11}$$

After the calculation, the CSPM [79] was corrected, as shown in the next formula:

$$h\_{-}(i) = \frac{h\_{-}t(i)}{\underline{\omega m\\_h\\_t(i)}}\tag{12}$$

#### *3.2. Particle Search*

In this paper, the particle search technique [80,81] was used as a separate module to prepare for the calculation of parameters such as water depth, acceleration, and velocity.

In the particle search, the mesh was firstly divided, and then the mesh area of each particle search was calculated; in order to ensure the symmetry of particle interaction, less than *2hi* or *2hj* was used as the judgment condition of *i* effective particles.

The specific steps of the particle search technique adopted [82] (see Figure 5) are described as follows:


**Figure 5.** Framework of particle search [81].

To avoid high time consumption caused by repeated particle search in the meshless SPH-SWE model, Algorithm 3 was produced.

In the open source code, the information of particles was stored in a three-dimensional array, and the grid was divided by a maximum smooth length of *2hmax* (by adopting this, all particles can be stored into a cell, causing failure of memory allocation; however, the particle search method mentioned above solved this issue).

#### *3.3. Water Depth Calculation*

According to the Newton-iteration method, the water depth of particles was calculated, and the maximum number of iterations and the iteration-errors were taken as a criteria to terminate the iterations. In order to reduce the calculation time while ensuring the accuracy of the results, the maximum number of iterations of each particle was generally set to 50, and the iteration error was setup to 10−3. Nevertheless, different approaches can be selected according to different calculation models.

Before each calculation step, the water depth and smooth length of each particle were guessed. At the same time, the smooth length *h* and the correction coefficient α*<sup>i</sup> <sup>k</sup>* were re-calculated and updated. In the same time step, the updated smooth length was then used for the sub-sequent SPH interpolation. The calculation framework is displayed in Algorithm 4.


After each water depth calculation, each particle was judged on whether the error requirements were met, and re-calculation was then completed for the particles that did not meet them. After this iteration cycle, the water depth, speed of the water, and volume were constantly updated.

**Algorithm 4**: Water depth calculation. 1. *Stage 1: Guess for density and smoothed length* 2. *!\$OMP PARALLEL DO PRIVATE(private variable),&* 3. *!\$OMP& SHARED(shared variable),DEFAULT(none),SCHDULE(static)* 4. *do i* = 1 → *total\_number\_of\_fluid particles* 5. *if particle\_i is valid then* 6. *1a: rhop*(*i*) = *rhop*(*i*) + *dt* ·*rhop*(*i*) · *ar*(*i*) 7. *1b: h*\_var(*i*) = *h*\_var(*i*) − (*dt*/*dm*)·*h*\_var(*i*)·*ar*(*i*) 8. *endif* 9. *enddo* 10. *!\$OMP END PARALLEL DO* 11. *CALL particle search() %Search particles* 12. *Stage 2: Calculate depth* 13. *do while ((maxval(resmax) .gt. Minimum error) .and. (Iterationtimes .lt. max)* 14. *!\$OMP PARALLEL DO PRIVATE(private variable),&* 15. *!\$OMP& SHARED(shared variable),DEFAULT(none),SCHEDULE(static)* 16. *do i* = 1 → *total\_number\_of\_fluid particles* 17. *if particle\_i is valid then* 18. *CALL PURE fluid particle(i,rhop\_sum(i),alphap(i))* 19. *CALL PURE virtual particle(i,rhop\_sum(i),alphap(i))* 20. *CALL PURE open boundary particle((i,rhop\_sum(i),alphap(i))* 21. *%Calculate next step's water depth and the smooth length* 22. ϕ*<sup>k</sup> <sup>i</sup>* = *rhop*(*i*) ρ*k i* − *rhop*\_*sum*(*i*) ⎡ ⎢⎢⎢⎢⎣ - *j mjWi*(*xi* − *xj*, *hi*) ⎤ ⎥⎥⎥⎥⎦ 23. α*<sup>k</sup> <sup>i</sup>* = *alphap*(*i*) ⎡ ⎢⎢⎢⎢⎣<sup>−</sup> <sup>1</sup> ρ*idm* - *j mjrij dWi drij* ⎤ ⎥⎥⎥⎥⎦ 24. ρ*k*+<sup>1</sup> *<sup>i</sup>* = *rhop*(*i*) <sup>1</sup> <sup>−</sup> <sup>ϕ</sup>*<sup>k</sup> i* ϕ*k <sup>i</sup>* +ρ*<sup>k</sup> i* α*k i* 25. *hi* = *ho* <sup>ρ</sup><sup>0</sup> ρ*k*+<sup>1</sup> *i* 1/*dm* 26. *endif* 27. *enddo* 28. *!\$OMP END PARALLEL DO* 29. **enddo**

#### *3.4. Velocity Calculations*

In order to calculate the acceleration ( *a* ) caused by internal force, the gradient of velocity and water depth had to be calculated, and the kernel function was adopted to complete this task as shown below:

$$\nabla p\_i = \sum\_j V\_j (p\_i - p\_j) \times \overleftarrow{\nabla} W\_i (\mathbf{x}\_i - \mathbf{x}\_j, h\_i) \text{ (\$p\_i = d\$ / \$u\$ / v\$)}\tag{13}$$

where *pi* is the depth/velocity of particle *i*. The calculation conducted for this step is displayed in Algorithm 5.

If a variable time-step was implemented, the next step was to calculate the time-step according to Equation (9). The calculation framework of the time step included a loop and no sub-routine. It was found that the speed-up of time step parallel computing was less than 2 to perform serially variable time step calculations.


*3.5. Calculation of Fluid Particle Acceleration, Riverbed Scouring, Speed, and Displacement*

The acceleration ( *t* ) caused by the internal force was firstly calculated, followed by the riverbed gradient and the total acceleration ( *a* ). After these calculations, the velocity and the displacement of fluid particles needed to be regularly updated as shown in the calculation framework Algorithm 6.

**Algorithm 6**: Calculation of acceleration, velocity, and position. The Lagrangian equation of motion for a particle i is d/dt ∂L/(∂v\_i )-∂L/(∂x\_i )=0, where the Lagrangian functional L is defined in term of kinetic energy K and potential energy π as L = K-π, where π is a function of particles position but not velocity.

```
1. Stage 1: Calculate →
                       t i(ax(i)/ay(i))
```

```
2. !$OMP PARALLEL DO PRIVATE(private variable),&
```

```
3. !$OMP& SHARED(shared variable),DEFAULT(none),SCHDULE(static)
```

```
4. do i = 1 → total_number_of_fluid particles
```

```
5. if particle_i is valid then
```

```
6. 1a. use Riemann solution to calculate →
                                           t i
```

$$\begin{array}{ll} 12. & \\ 13. & \text{Stage 2: Calculate } \nabla b\_{l} = \sum\_{i} b\_{j} \widetilde{\nabla} \mathcal{W}\_{l} (\boldsymbol{x}\_{i} - \boldsymbol{x}\_{j}, h\_{i}) V\_{j} \end{array}$$

$$
\begin{array}{ccccccccc}
\dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\
\dots & & & & & & & \dots & \dots & \dots \\
\dots & & & & & & & & \dots \\
\end{array}
$$

$$14. \text{
\$\text{Stage 3: Calculate } \stackrel{\rightarrow}{a}\_{I} = \frac{\stackrel{\rightarrow}{s} + \stackrel{\rightarrow}{\vec{\nu}}\_{i} \cdot \stackrel{\rightarrow}{\vec{k}}\_{i} \cdot \stackrel{\rightarrow}{\vec{\nu}}\_{i} + \stackrel{\rightarrow}{\vec{\pi}}\_{i} \cdot \stackrel{\rightarrow}{\nabla}b\_{i}}{\text{I} + \stackrel{\rightarrow}{\nabla}b\_{i} \cdot \stackrel{\rightarrow}{\nabla}b\_{i}} \text{V}b\_{i} - \stackrel{\rightarrow}{\mathbf{f}}\_{I} + \stackrel{\rightarrow}{\mathbf{\tilde{S}}}\_{f,I}$$


If the open boundary was adopted, the displacement was updated, and it was checked whether it became a fluid particle or a buffer zone; if the case under investigation had a particle splitting zone, the fluid particles that meet the conditions identified were split [37]. Achieved this landmark, the calculation process was then completely repeated according to Algorithm 1.

#### **4. Applications**

#### *4.1. Validation 1: 2-D Dry Bed Dam Break with Particle Splitting*

In order to test the performance of the new computing SPH-SWE framework according to the CPU-OpenMP parallel computing, the open source case "*2-D dry bed dam break with particle splitting*", referred to as DBDB case [67,68], was considered. Initial conditions of this open source DBDB case were setup as follows: area of 2.6 m × 1.2 m, initial fluid particle layout of 0.8 m × 1.2 m, and initial water depth equal to 0.15 m, as shown in Figure 6.

**Figure 6.** 2-D dry bed dam break with particle splitting (DBDB) case setup.

The spacing of the fluid particles was 0.015 m × 0.015 m, 0.01 m × 0.01 m, and 0.005 m × 0.005 m, respectively. Table 3 shows the number of fluid particles, virtual particles, and riverbed particles selected for the 3 cases tested (considering that no inflow and outflow conditions were set, the number of open boundary particles was maintained equal to 0).


**Table 3.** Particle numbers for each case tested (unit: pcs).

The particles arrangements in the three cases were calculated using the open source code [67,68] and the CPU-OpenMP parallel code (in terms of 1000 particles per thread and 2000 particles per thread) to ensure that the calculation results were consistent and comparison of each model's performance could be completed.

For t = 0, the instantaneous burst occurs and the fluid flows downstream. Figure 7 shows flow velocities for T = 1.2 s under the three fluid particle arrangements displayed in Table 3. It can be seen from Figure 7 that the more fluid particles there were, the more water flow characteristics and velocity distribution characteristics after dam break were seen. Figure 8 shows the comparison between numerical and experimental results for the three cases tested in Table 3.

It can be concluded that with the increase of fluid particles (Case 3), the error between the numerical and the experimental results was smaller.

Hence, bigger is the number of fluid particles computed, and more valuable are the water depth and velocity calculations for each time step and location across the dam system, therefore providing more support for the formulation of dam break mitigation plans. This also reflects the superiority of SPH-SWE model in dealing with large deformation and free surface problems.

**Figure 7.** Results of 2-D Dry-Bed Dam Break with particle splitting at 1.2 s obtained from the simulations with: (**a**) 4374 particles; (**b**) 9801 particles; (**c**) 38,801 particles.

**Figure 8.** Numerical vs. experimental results comparison (velocities on the **left** and water depths on the **right**).

The time and speedup of total time *t*(*s*) (Equation (14)) required for calculating water depth, acceleration, speed, and displacement were quantified according to open source code and CPU-OpenMP parallel code (under different thread configurations) as shown in Table 4.

$$t(s) = \frac{t(k)}{t(c)}\tag{14}$$

where *t*(*k*) represents the running time of open source code and *t*(*c*) represents the running time of the parallel code. In Table 4, *R*(*t*) represents the particle search time; *C*(*t*) represents the time to calculate the water depth; *A*(*t*) represents the time to calculate the acceleration, speed, and displacement; *T*(*t*) represents the total running time of the case (including *t*(*k*) and *t*(*c*)); the time unit is always seconds (s). It can be seen from the results that the number of particles computed was higher, and the parallel calculation was larger, based on CPU-OpenMP (Figure 9).

**Figure 9.** Speedup of total runtime for validation 1.

CPU-OpenMP allocated a thread according to 1000 fluid particles in parallel and calculated case 2 and case 3, in 26.7 s and 91.36 s, respectively.

In parallel computing, *S*(*p*) (the speedup ratio) and *E*(*p*) (parallel efficiency) were important indexes to evaluate the parallel effect. *S*(*p*) was the ratio between the serial time and the multi-core parallel time when threads calculate (*p*) and solved the iteration at the same time, as follows:

$$S(p) = \frac{T\_s}{T\_p} \tag{15}$$

where, *Ts* is the time spent by a single processor in the serial mode; *Tp* is the time spent by threads (*p*) in the parallel mode. *E*(*p*) (parallel efficiency) is the ratio of the acceleration ratio to the number of CPU cores used in the calculation (and *E*(*p*)≤ 1), indicating the average execution efficiency of each processor. When the acceleration ratio was close to the number of cores, the parallel efficiency was higher, and the utilization rate of each thread was higher, as calculated below.

$$E(p) = \frac{S(p)}{p} \tag{16}$$


**Table 4.** Run time in different configurations for 2-D Dry-Bed Dam Break with particle splitting.

In order to check the parallel effect of CPU-OpenMP, case 3 was calculated with different threads. The acceleration ratio and parallel efficiency of different threads are shown in Table 5, and the performances of parallel algorithms are displayed in Figure 10.


**Table 5.** Case 3 speed-up and parallel efficiency under different threads.

According to Table 5 and Figure 10, as the number of enabled threads increased, the speedup ratio, parallel efficiency and calculation time were affected by the following trends: (1) the calculation time decreases with the increase of threads, but when the number of online processes exceeded 10, the time-consuming reduction speed changed from fast to slow reaching towards a balance; (2) the acceleration ratio increased all the time, but the improvements varied, and after 10 threads, the increase rate was from fast to slow; (3) the parallel efficiency decreased with the increase of threads, but the decrease rate fluctuated. When calculating the number of 10 threads, the parallel efficiency started to be less than 90%; therefore, case 3 could allocate one thread according to 5000 particles, with the acceleration ratio of 7.47 and the parallel efficiency of 93.35%.

**Figure 10.** CPU-OpenMP model parallel computing performances.

#### *4.2. Validation 2: 2-D Dam Break with A Rectangular Obstacle Located in the Downstream Area*

In order to solidify the accuracy and advantages of the new SPH-SWE model proposed calculated by CPU-OpenMP, a second case has been considered for validation. This case involved the dam break flow with a rectangular obstacle located in the downstream area as shown in Figure 11 [2–83].

**Figure 11.** Scheme of the second model used for validation of the new SPH-SWE model proposed [2–83].

For this second validation scheme, the fluid particles were arranged according to the particle spacing of 0.01, 0.005, and 0.002 m. Figures 12 and 13 display the results for t = 0.74 s and t = 1.76 s for the same particle spacing used by Gu et al., [2] (0.01 m) and increasing number of fluid particles involved.

In Table 6, it is possible to check the number of particles involved in each simulation, and it can be noticed that the speed up of total time (*t*) slightly increased with the rise of particles (using 8 threads in all three processes) (Figure 14).

**Table 6.** Particle numbers for each case tested (unit: pcs). Remarks: In all three cases, 8 threads were used for parallel calculation.


**Figure 12.** Velocity distribution at 0.74 s for Case 4 (**a1**), 5 (**b1**) and 6 (**c1**) displayed in Table 6.

**Figure 13.** Velocity distribution at 1.76 s for Case 4 (**a2**), 5 (**b2**) and 6 (**c2**) displayed in Table 6.

**Figure 14.** Speedup of total runtime for validation 2.

Table 7 unveils values of *R*(*t*), which represents the particle search time; *C*(*t*), which represents the time to calculate the water depth; and *A*(*t*), which represents the time to calculate the acceleration, speed, and displacement. It can be seen from the results (Figure 14) that when the number of particles computed was higher, the parallel calculation was slightly larger, based on CPU-OpenMP.


**Table 7.** Run time in different configurations for 2-D Dam Break with a rectangular obstacle downstream.

The simulation results are consistent with the previous case to demonstrate the improvement made by the proposed SPH-SWE calculation framework. When analyzing the same parameters as in the paper by Gu et al., [2], 12,423 fluid particles were involved. Results displayed in Figure 15 confirmed that the agreement between numerical results and experimental datasets improved even when simulating an increase of the number of fluid particles. However, the improvements do not involve the entire domain because in some positions (for example (a) H4 gauge, x = 1.3 m–1.8 m; (b) H2 gauge, x = 3.50 m–4.5 m) there are still minor inaccuracies (caused by the truncation error of the kernel function and the processing of boundary particles in the SPH method), even with the model with 323,145 fluid particles; hence, future work will focus on improving the algorithm to progress the calculation accuracy.

**Figure 15.** Comparison between experimental and numerical results water depth datasets: (**a**) gauge H4 in [2] and [80]; (**b**) gauge H2 in [2] and [80].

#### **5. Conclusions**

Vacondio et al. [34–37] made an open source version called SWE-SPHysics, which has been optimized and adapted based on the hydrodynamics investigated by other researchers during the last decade [38–47]. However, despite continuous progress, there was still a limitation related to the computational efficiency when the number of particles to simulate is very large, and this aspect still needs to be improved.

To fill this gap, in this study, a new solution to the SPH-SWE model introduced by Vacondio et al. [34–37] was proposed, and it was validated against two open source case studies of a *2-D dry-bed dam break with particle splitting* [67,68] and of a *2-D dam break with a rectangular obstacle downstream* [2,83]. To test the computing performance against the first case study, when involving large numbers of particles, three cases, involving different particles numbers, were tested (case 1—4374 particles; case 2—9801 particles; case 3—38,801 particles). Furthermore, this paper adopted the parallel computing method of CPU-OpenMP that is applied to a single machine and a multi-core to calculate the new SPH-SWEs framework.

By applying this CPU-OpenMP method, results have confirmed that the computing speeds of case 1/case 2/case 3 were increased by 4.01 times/6.14 times/7.17 times, respectively, to compute the new solution framework of the SPH-SWE model proposed to the open source case study previously mentioned [67,68]. According to the new solution framework of SPH-SWE model, case 3, characterized by the highest number of particles, was also calculated by using different threads. It was found that the speedup ratio can reach 7.47 when the parallel efficiency was more than 90%, which fully proves the good calculation performance of the CPU- OpenMP parallel new SPH-SWE model. Additionally, in the new solution framework of the SPH-SWE model proposed, particle search was used as a separate module for parallel computing, which greatly improved the computing efficiency and could replace the meshless SPH-SWE model in the open source code [67,68].

Therefore, using CPU-OpenMP parallel computing demonstrated that the SPH-SWE model new framework can accurately and timely simulate the flood evolution after a dam break.

In future works, the SPH-SWE model can be put into existing clusters to achieve more threads and further improve the calculation efficiency. Furthermore, this would enable the possibility of introducing more effective new algorithms into the SPH-SWE model (i.e., debris flow or water pollution modules) in order to expand its application. Continuous development of technology aids the improvement of new tools to design and inspect more accurate solutions, and this is an area in continuous development that needs to be addressed to support local and national authorities in making decisions to mitigate drastic effects generated by dam break.

**Author Contributions:** All authors contributed to the work. Conceptualization, S.G. and L.T.; methodology Y.W., M.R. and T.Y.; validation, Q.Z.; formal analysis, S.G., Y.W., M.R. and Z.X.; investigation, L.T., S.G., and M.R.; resources, S.G.; data curation, L.T. and M.R.; writing—review and editing, P.C., X.W., and M.R.; visualization, Q.Z.; supervision, S.G., M.R.; project administration, S.G.; funding acquisition, S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the following projects: National Key R&D Program of China (No. 2017YFC0404303), National Natural Science Foundation of China (No.51869025, 51769028, 51868066), Qinghai Science and Technology Projects (No. 2018-ZJ-710), Youth Fund of Qinghai University (Grant No. 2017-QGY-7), National Key Laboratory Project for Water Sand Science and Water and Hydropower Engineering, Tsinghua University (Grant No. sklhse-2018-B-03), Beijing Institute of Structure and Environment Engineering Fund(Grant No. BQ2019001).

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).
