Numerical Simulation of Imported Sediment in a Stilling Basin

Sediment buildup at the bottom of a stilling basin can result in premature drainage of spillway structures and can even lead to dam failure in severe cases. Such failures pose ecological and human safety hazards to downstream areas. To evaluate the sudden discharge and potential dam failure associated with sediment buildup, we developed a two-dimensional two-phase flow simulation model built on a particle-based force balance equation. We compared the flow patterns and energy dissipation effects in the stilling basin at different inlet flows (2, 3, 4.5, and 6.75 m2/s), and the subsequent bottom deposition was compared across different sand discharge mass flow rates (0.1, 0.2, and 0.3 kg/s). The results show that the turbulent energy increased with the increasing inlet unit width flow rate. When more vortices were generated and the flow velocity was reduced significantly, the energy dissipation was more effective. The sediment deposition at the bottom of the stilling basin gradually increased with the decrease of inlet unit width flow and the decrease of the sediment mass flow rate. Meanwhile, at a fixed inlet shape, the change in inlet unit width flow had little effect on the maximum sedimentation height at the bottom of the basin. In addition, the average deposition rate at the bottom of the stilling basin was positively correlated with the inlet sedimentation concentration, and the correlation coefficient could be as high as 0.97. In this two-phase flow method, the error of the simulated value over the theoretical value was less than 10%. This simulation of sediment deposition at the bottom of the stilling basin provides a practical reference for dam managers.


INTRODUCTION
More than 60% of China's Loess Plateau was once subject to severe soil erosion, which caused riverbed uplift and erosion in the lower reaches of the Yellow River (Shi and Shao, 2000;Xin et al., 2012). The widespread nature of soil erosion means that different measures are taken to mitigate different dimensions of damage (Cerdà et al., 2009). Check dams are one such mitigation measures, and they have been constructed on streams to trap soil and to retain flood waters (Ran et al., 2008). Usually, these key dams are equipped with a spillway and other discharge structures that are connected to a dissipation pond (Tsujino et al., 2010). Stilling basins are an example of a common energy dissipation facility in water conservancy projects that reduce both the energy of discharged water and the loss of downstream equipment (Xie et al., 2016).
Many scholars (Cheng and Liu, 2011;Wobus et al., 2011;Qiu et al., 2012;Wu and Mu, 2012;Javan and Eghbalzadeh, 2013) have developed simulations to provide a reference for basin shape optimization in actual projects. For example, Speziale and Ngo (1988), Zheng et al. (2010), and Luo et al. (2012) used the RNG K-εturbulent flow model to test the design of the dissipation basin. Other scholars have used hydraulic model tests to study the hydraulic characteristics of energy stilling basins. Li et al. (2015) analyzed the water leap pattern in shallow water cushion stilling basins and concluded that the inlet shape of the stilling basin influences the depth of the water cushion. Liu (2012) realized that an increase in the length of the stilling basin reduces the fluctuation of water flow out of the basin, and therefore mitigates downstream scouring. Zhang and Zhao (2015) deduced the relationship between the coefficient of head loss along the hydraulic jump and the local head loss coefficient, and ascertained that the percentage of local head loss in the hydraulic jump area increases with the increased Froude number.
The average annual erosion of the Loess Plateau is as high as 10,000 km 2 year −1 (Shi and Shao, 2000;Li et al., 2019). When the flood water level is higher than the spillway elevation, upstream sediment is carried out of the release building by the flood water. However, sediment buildup at the bottom of the dissipation basin can lead to premature drainage of the spillway structures ( Figure 1), which in severe cases can lead to dam failure. Such failures pose ecological and human safety hazards to downstream areas. The gas-liquid-solid flow can be assessed with computational fluid dynamics (CFD) using the discrete particle method (DPM). The volume of fluid tracking is expressed as the volume of the fluid (VOF). Li et al. (1999) used a combination of CFD, DPM, and VOF to simulate gas-liquid-solid flow in fluidized beds. Chen et al. (2012) developed a CFD-DPM model to study the behavior of gas/ solid flow in the airways of patients with chronic obstructive pulmonary disease (COPD). Li et al. (2013) combined the VOF and DPM multiphase flow models to develop a model describing the gas-liquid two-phase flow in a top-bottom blowing steelmaking oxygen converter.
In sum, many scholars have studied the hydraulic characteristics of the stilling basin, but few scholars have considered the sediment deposition in the bottom of the basin. This study analyzes how the design of flood check dams affects sediment deposition in stilling basins. More specifically, we used the CFD-VOF-DPM model to simulate the deposition condition of the stilling basin. We analyzed sedimentation at the bottom of the stilling basin and determined the relationship between sedimentation volume and boundary condition. The amount of sand discharged from the stilling basin will inform the future flood control design of Loess Plateau check dams.

Control Equation
The RNG k-ε turbulence model is based on the improvement and modification of the standard model. The RNG k-ε adds an additional correction term to the ε equation, which improves the accuracy of calculating the rapid flow conditions.
In this study, we used the RNG turbulent flow model and the DPM discrete phase of FLUENT 16.0 software for numerical calculations. VOF was used for free interface tracking (Thinglas and Kaushal, 2008). The continuity equation, momentum equation and K, and the εequation of the RNG turbulent flow model can be expressed as follows: where ρis the average density of the volume fraction, μ is the molecular viscosity coefficient, t is the time, x i is the spatial coordinate in the i direction, u i is the velocity component in the i direction, B i is the force per volume, G k is the turbulent kinetic energy induced by the mean velocity gradient, k is the turbulent kinetic energy, ε is the turbulence dissipation rate, and σ ε and C 1ε are the constant values used in the turbulence model (σ ε 0.718, C 1ε 1.68).

DPM Model
The orbits of discrete-phase particles (Akhtar et al., 2007) were solved in FLUENT by integrating the differential equation for the forces acting on the particles in the Rasch coordinate system. The Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 659855 equilibrium equation for the forces acting on the particles (particle inertia various forces acting on the particle) in the Cartesian coordinate system had the following form (x-direction): Re where F D represents the mass traction of particles, u represents the fluid phase velocity, u p represents the particle velocity, g x represents the gravity acceleration in the x direction, m is the hydrodynamic viscosity, r represents the fluid density, r p is the particle density, d p represents the particle diameter, and Re represents the relative Reynolds number.
The traction coefficient C D can be expressed as follows: For spherical particles, a 1 , a 2 , and a 3 in the above equation are constants for a range of Reynolds numbers. The other forces included in the force balance equations for particles may be important in some cases. The most important of these "other" forces is the "apparent mass force" (additional force), which is the additional force caused by the acceleration of the fluid around the particle, and it is expressed as follows: When r >r p , the apparent mass force cannot be ignored. The additional forces due to the fluid pressure gradient present in the flow field are as follows: According to the soil particle classification of test soils derived from the U.S. Department of Agriculture (USDA), the Loess Plateau has a relatively large proportion of silty soil (0.02-0.002 mm). The particle diameter of the simulated soil was set at 0.02 mm using homogeneous sediment particles. The outlet section is assumed to be a fully developed turbulent flow. The inlet section is a given water-level value, and the turbulence energy K is used in the empirical formula. The turbulence intensity I is calculated using the following equation: where urepresents the average flow velocity, I is the turbulence length scale (I 0.07 L), L represents the hydraulic diameter, and C μ represents an empirical constant (C μ 0.09).
Sediment particles lose some kinetic energy when colliding with the wall. To account for this, we defined the geometric bottom of the model as a reflective wall. The recovery coefficient after collision indicates the kinetic energy loss of the sediment particles. We set the normal phase and tangential recovery coefficients of the particles colliding with the wall, assuming the same amount of energy is lost in each collision (Ye, 2019). The specific formulas are as follows: where e n represents the normal phase recovery coefficient of sediment particles colliding with the wall, e t represents the tangential phase recovery coefficient of the collision between sediment particles and the wall, and θrepresents the direction of sediment particles colliding with the wall and the angle of the collision surface. In the DPM model, the deposition rate is equal to the ratio of wall deposition mass flow rate to the area, according to the particle mass balance equation. The equation is as follows: where R A represents the sediment deposition rate, m p represents the mass flow rate of the particle stream, and A face represents the area of the wall at the particle impact boundary.

Solutions
In this test, the finite volume method is used as the control equation, and the second-order implicit scheme is used in the discretization of the time item. The PISO algorithm is used to solve and control the coupling of velocity and pressure in the equation. The VOF method is used to track and simulate the free surface and two-phase flow of air and water. The free water surface is established by the geometrical reconstruction scheme.

Computational Domain and Boundary Conditions
Using drone aerial photography technology, Figure 2 is taken from the top of the stilling basin. The stilling basin is located at The main ditch of Guandigou is 18 km long, with an average slope gradient of 1.15%, a "V" shaped ditch profile, a gully density of 5.34 km/km 2 , and an elevation of 820-1,180 m. The topography of the watershed is fragmented, the gullies are crisscrossed, and soil erosion is a serious issue. The average soil erosion modulus is 14,000 t/(km 2 a) along the middle and lower reaches of the Wuding River. The model was built at a 1:1 scale and consisted of a steep slope, a dissipation pond, a gradual section, and a tail channel ( Figure 3). As shown in Figure 3, the overflow weir connected with the steep slope was 24 m high, so the starting point of the steep slope was 24 m higher than the ground elevation, with a slope drop ratio of 1:5. The steep slope was 120 m long, the stilling basin was 11 m long and 2.2 m deep, the gradient section was 8 m long, and the tail channel was 20 m long, with a slope drop ratio of 0.025.
The model's geometry consisted of a quadrilateral grid defined by ICEM CFD 16.0. Grid accuracy varied based on proximity to the extremes of the calculation domain. The grid near the bottom of the calculation domain was more precise (0.02 mm), while the grid near the top of the calculation domain was less precise (0.10 mm). We accounted for this variation because the steep slope had only a thin layer of water flow. The entire grid consisted of 285,100 cells ( Figure 3). Its grid average mass was greater than 0.92, and no negative volume appeared. This means that the model had great mesh quality. Simulations were carried out using the Fluent 16.0 commercial package. The uncoupling arithmetic method was used to separate and solve model equations: pressure-velocity coupling was defined by the pressureimplicit with splitting of operators (PISO) algorithm. The left side of the calculation domain represented the point of traffic entry, while the right side of the calculation domain represented the pressure outlet. The top of the calculation field represented the pressure, and the bottom of the calculation domain represented the reflective wall.

Working Conditions
The test is based on a spillway designed for a 20-year flood. The bottom width of the check dam's spillway trapezoidal section is 10.5 m, and the side slope ratio is 1:1. The unit width flow rate of the spillway inlet is the ratio of the maximum spillway flow rate to the bottom width of the inlet section. The change in the spillway discharge flow rate is in 1.5 times increments (Table 1). To obtain the maximum flow rate of the spillway, we consulted the existing check dam design information. The calculation formula is as follows: where Q DF represents the maximum flow rate of the spillway, Q F represents the 20-year design flood flow (105 m 3 s −1 ), V SF represents the stagnant flood volume (550,000 m 3 ), and W F represents the maximum flood releasing capacity of the spillway (100,000 m 3 ). The annual sand discharge of the reservoir does not exceed 40%. The average annual sand discharge is 20% (Xin et al., 2012;Wang et al., 2020). The erosion modulus of the Loess Plateau is 1000t/(km 2 · year −1 ). The formula for annual sediment discharge is as follows: where M 0 is the annual sand discharge (kg), F is the check dam control watershed area (F 5 km 2 ), and K is the annual average erosion modulus [k 1,000 t/(km 2 year −1 )]. Assume that the annual sand discharge mass is obtained during one rainstorm and all the sediment is discharged by the spillway. The unit width sand mass flow rate is then obtained by dividing the sand discharge mass flow rate by the  bottom width of the inlet section. Different unit width sand discharge mass flow rates were defined for each working condition, including CK (clear water), CL 1 (M 0.1 kg s −1 ), CL 2 (M 0.2 kg s −1 ), and CL 3 (M 0.3 kg s −1 ). The median particle size of the sediment in the Yellow River field sub-high sand-bearing flood is generally 0.01-0.03 mm (Wang et al., 2020). We used homogenous sand with a grain size of 0.02 mm for the simulation.

Verification of the Mathematical Model
Jiuyuangou watershed is located in a temperate semiarid region. The average precipitation is 469 mm, with rainfall mostly concentrated over 6-9 months, falling in heavy rainfall events. The main function of the check dam stilling basin is to dissipate the energy of the rising water in front of the dam during heavy rainstorms, thus protecting the downstream farmland from being washed away. Data on the check dam's stilling basin were obtained from the stilling basin flood control center, which provided data on the sequent water depths when the flow rates were 2.5, 5, 7.5, and 10 m 3 /s. Because these flow rates are low, the sediment content added to the stilling basin was ignored in the calculations. Table 2 shows the measured and empirical formulas for calculating values at different flow rates, and the empirical formula is shown in Eq. 21. The error of the conjugate water depth ratio of the measured value to the calculated value is within 6%. Because the error is small, this indicates that the true value can be replaced by the calculated value of the empirical formula. Second, Figure 4A shows the simulated and empirical formula-calculated values under CK 3 processing. The horizontal coordinate is the Froude number (Fr), and the vertical coordinate is the conjugate bathymetry ratio. The errors of the simulated and empirical formula-calculated values are within 10%. This shows that the operating conditions are reliable in clear water and that these can be used for other operating conditions as well. Finally, the sediment deposition at the bottom of the pool could not be obtained accurately due to the unstable water flow. Therefore, to verify the reliability of the sediment deposition simulation, the settling velocity of the modeled (ω m ) single-particle sediment was compared with the existing empirical (ω c ) equation (Camenen, 2007). Camene proposed a simple, reliable, and general equation for the settling velocity of the particles that accounts for their shape and roundness. The vertical coordinate of Figure 4B shows the ratio of calculated (ω c ) to simulated values (ω m ), which is lower than 1.3. This indicates that the simulated results are similar to those theoretical values calculated by the empirical formula. where q is the pre-jump section unit width flow (m 2 · s −1 ), h'is the pre-jump water depth (m), h''is the post-jump water depth (m), and η is the conjugate water depth ratio.

Changes in the Energy Dissipation Rate and Hydraulic Jump Pattern
To calculate the energy dissipation rate (K j ) of the stilling basin at different flow rates, we compared the energy of the initial section (E 0 ) at the bottom of the stilling basin to the energy of the exit section (E 1 ) at the tail channel. The total energy of each section can be expressed asE i Z i + pi c + v 2 i 2g . Here, Z i is the potential energy, p i is the average pressure, and v i is the average flow velocity of the cross section. The energy dissipation rate of each section can be expressed as η E1−E0 E0 . The energy dissipation rate increased with increasing Fr (Table 3), and the energy dissipation rate was greater in clear than in sandy water. The reason is that the sand-bearing water flow accelerates the flow velocity of the slope surface water flow so that the flow velocity of the water flow into the stilling pool is faster. At a constant unit width flow rate, the energy dissipation rate increased as the sand discharge mass flow rate increased. This shows that under the test conditions, the sand concentration of the flow affected energy dissipation: when Fr > 4.5, K j was greater   Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 659855 than 50% and when Fr < 4.5, K j was less than 50%. Figure 5 illustrates the water phase cloud at the bottom of the stilling basin for different inlet unit width flows. As the inlet unit width flow increased, the flow pattern at the bottom of the basin destabilized, thereby increasing energy dissipation. At a unit width flow rate less than 3 m 2 /s, the flow pattern was more stable. Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 659855 7 Change in the Sediment Deposition Height at the Bottom of the Stilling Basin Figure 6 shows the sediment deposition at the bottom of the stilling basin. The sediment was first deposited at the end of the pool bottom. Then the deposition height and length gradually increased. After reaching the maximum deposition height of 0.2 m, the deposition length continued to increase until the pool bottom was covered ( Figures 6D2,D3). In addition, as the mass flow rate increased and the inlet unit width flow rate decreased, the maximum deposition rate at the bottom of the pool increased. As the inlet unit width flow rate decreased, the flow pattern of the basin bottom gradually stabilized ( Figure 5). This indicates that the more stable the flow pattern of the pool bottom, the lower the Fr and the greater the sediment deposition at the bottom of the pool.

Analysis of Factors Influencing the Amount of Sedimentation on the Basin Floor
The effects of dynamic and sediment parameters on the wall deposition rate were analyzed using the Pearson correlation analysis (Figure 7). The inlet flow rate (q 0 ) was found to be strongly and positively correlated with the inlet velocity (v 0 ), turbulent kinetic energy (E), and shear force (τ), with a correlation coefficient as high as 0.92. The deposition rate (R A ) at the bottom of the pool was negatively correlated with the inlet flow, inlet velocity, turbulent energy, and shear, with correlation coefficients greater than 0.73. The deposition rate was positively correlated with the mass flow rate of sand discharge (M), inlet sediment concentration (C 0 ), and average sediment concentration (C). The correlation coefficient between sedimentation rate and inlet sediment concentration was as high as 0.97. Thus, the relationship between sedimentation rate and inlet sediment concentration was linearly fitted, with R 2 0.94 (Figure 8). The inlet sediment concentration is the ratio of the mass flow rate of discharged sand to the inlet flow rate. The final equation for the relationship between the average deposition rate and the mass flow rate and inlet unit width flow rate is as follows: where y A is the deposition rate (kg · m −2 ), M is the mass flow rate of sediment discharge (kg · s −1 ), q 0 is the inlet unit width flow rate (m 2 · s −1 ), C 0 is the inlet sediment concentration per unit area (kg · m −2 ), and H 0 is the water inlet height (m).

Energy Dissipation Effect of the Stilling Basin
The energy dissipation rate is an important index for measuring the effect of energy dissipation. Based on the analysis of the energy dissipation rate in a previous article, many scholars believe that Fr has a significant effect on the initial energy dissipation rate. Li et al. (2018) believes that the larger the Fr, the higher the energy dissipation rate. Zhang et al. (2017) proposes that when Fr is 4.5-9.0, the energy dissipation effect improves because the water jump is stable. Sun et al. (2019) further refined the range of the effect of Fr on the initial energy dissipation rate and concluded that when Fr is 2.55-4.5, the energy dissipation rate (Kj) ranges from 20 to 45%, and when Fr is 4.5-9, the energy dissipation rate (Kj) ranges from 45 to 85%. Here, we found that  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 659855 the variation in the energy dissipation rate (Kj) between Fr in the pre-leap (E0') and the post-leap section of the dissipation pool basically conformed to this law, even under different working conditions. When Fr was > 4.5, the energy dissipation rate Kj was greater than 45%. When Fr was < 4.5, the energy dissipation rate Kj was less than 47%. Also, in Figure 9, it is shown that the higher the inlet flow, the higher the turbulent kinetic energy, the more vortices are generated, and the greater the energy dissipation rate. We further believe that the more vortices, the more significant the flow rate reduction and the better the energy dissipation effect. This is consistent with the conclusion of other experts as well. Tan et al. (2020) believed that the surface vortex zone is significantly larger than the bottom vortex zone in the second-stage force elimination pool, and that the water flow is more stable and the energy dissipation rate is larger than that in the single-stage pool. Dong et al. (2016) found that a noncomplete wide tail pier force pool can produce a large transverse velocity gradient and that this transverse velocity gradient produces additional turbulent shear and lateral flow. In this way, the energy dissipation effect is improved, as desired.

Sediment Deposition at the Bottom of the Stilling Basin
The maximum lift force generated by turbulent pressure fluctuations acting on the bottom of the stilling basin results in poor sediment stability (Bowers and Tsai, 1969). If the lifting force is greater than the gravity of the sediment particles, no sediment will be deposited. If the lifting force is less than the gravity of the sediment particles, the sediment will be deposited on the bottom of the basin. The bottom surface of the water flow is gradually transformed from a relatively smooth basin to a submerged sediment topography with a certain undulating height, and the resulting additional friction at the bottom leads to a relative increase in water flow resistance (Zhu, 1982;Wei 2013).
As the inlet unit width flow rate increased, the siltation height at the bottom of the basin did not exceed 0.2 m. The effect of inlet unit width flow on the maximum silt thickness was insignificant. Chen et al. (2019) predicted the deposition of raw Yangtze River water into the basin's transfer pond, concluding that when the shape of the inlet is fixed, the sediment deposition distribution does not change much. Guo (2017) analyzed the sediment deposition pattern of the river and determined that the maximum deposition of the river would not exceed 0.03 m, and the deposition in most areas was below 0.017 m. Our results are consistent with Chen's conclusion, suggesting that a fixed inlet shape results in decreased significance of flow pattern on maximum sedimentation height at the bottom of the basin. Flow is the main factor affecting the sand discharge rate at the outlet of the stilling basin. The flow rate has a significant effect on the amount of sediment siltation at the bottom of the stilling basin. By simulating sediment deposition in inverted siphons, Bagchi (2012) concluded that the deposited sediment gradually decreases as the flow rate increases.

CONCLUSION
In this study, we analyzed the effects of different flow rates, sand discharge on the flow pattern, and sediment deposition at the bottom of the stilling basin. The VOF-DPM model was used for the simulation. After validating the simulations with empirical equations, we analyzed the energy dissipation effect and deposition law of the stilling basin under different inlet flow rates (2, 3, 4.5, and 6.75 m 2 /s) and different sand discharge mass flow rates (0.1, 0.2, and 0.3 kg/s). When the mass flow rate was 0.3 kg/s, the deposition height at the bottom of the basin remained at 0.2 m even as the inlet unit width flow rate increased. With an inlet unit width flow rate of 2 m 2 /s, deposition occurred as far as 11 m from the stilling basin discharge point. If the maximum instantaneous lift force generated by the turbulent flow pressure fluctuations acting on the bottom plate is less than the sediment's own gravity, the sediment will be deposited. Deposition along the bottom of the basin increases the water flow resistance, and the correlation coefficient between the sediment deposition rate at the bottom of the basin and the inlet sediment concentration was as high as 0.97. Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 659855