3D Numerical Simulation of Gravity-Driven Motion of Fine-Grained Sediment Deposits in Large Reservoirs

Abstract

Deposits in dam areas of large reservoirs, which are commonly composed of fine-grained sediment, are important for reservoir operation. Since the impoundment of the Three Gorges Reservoir (TGR), the sedimentation pattern in the dam area has been unexpected. An integrated dynamic model for fine-grained sediment, which consists of both sediment transport with water flow and gravity-driven fluid mud at the bottom, was proposed. The incipient motion driven by gravity in the form of fluid mud was determined by the critical slope. Shallow flow equations were simplified to simulate the gravity-driven mass transport. The gravity-driven flow model was combined with a 3D Reynolds-averaged water flow and sediment transport model. Solution routines were developed for both models, which were then used to simulate the integral movement of the fine-grained sediment. The simulated sedimentation pattern agreed well with observations in the dam area of the TGR. Most of the deposits were found at the bottom of the main channel, whereas only a few deposits remained on the bank slopes. Due to the gravity-driven flow of fluid mud, the deposits that gathered in the deep channel formed a nearly horizontal surface. By considering the gravity-driven flow, the averaged error of deposition thickness along the thalweg decreased from −13.9 to 2.2 m. This study improved our understanding of the mechanisms of fine-grained sediment transport in large reservoirs and can be used to optimize dam operations.

1. Introduction

The impoundment of large reservoirs results in numerous changes, including decreases in water flow velocity, sediment deposition, and sorting [1]. In the reservoir, coarser sediments tend to be deposited in the upstream regions of the reservoir, advancing downstream steadily but slowly in the form of a delta, whereas fine sediments reach the reservoir in the form of suspended load, by means of homogeneous or stratified dense flow [2].

Sedimentation is an important issue that receives a significant amount of attention in reservoir construction and operations [3]. The excess of sedimentation in a reservoir leads to sediment entrainment in waterway systems and hydraulic schemes [4]. During the demonstration phase and early operational periods of the Three Gorges Reservoir (TGR), several hydraulic research institutes studied the problem of reservoir sedimentation [5,6]. After the impoundment, the sediment influx conditions of the TGR changed significantly, and the actual reservoir sedimentation conditions were different from the predictions made during the demonstration phase [3,7]. With the increase in the water level in front of the dam during flood seasons, the trap efficiency increases [8,9]. In the TGR, deposited sediment is mostly fine sediment with median diameters of less than 0.01 mm, which was not forecasted by previous studies [7]. In addition, fine-grained sediments with median sizes of approximately 0.004–0.010 mm have been observed in the dam area as early as the initial stage of the impoundment [10]. This sediment was deposited due to the significant increase in water depth of more than 100 m that occurred after the water storage operation. Although reservoir deposition is common, the sedimentation pattern in the dam area of the TGR is unusual. Generally, in the dam area, the distribution of sediment concentration is uniform, which would cause cross-sectional uniform deposition in the dam area, and then lead to an uneven surface across the section [11,12]. However, the deposits of fine-grained sediment in the dam area of the TGR have mainly accumulated in the deep channel and have a flat surface [13,14]. This has led to a rapid rise in the bed surface in the deep channel despite the small total amount of siltation. The elevations of the bottom discharge hole and deep discharge hole of the Three Gorges Dam (TGD) are 55 and 90 m, respectively. The bottom elevations of the water inlets of the power plant are 108 m. However, the surface elevation of the deposits in the deep channel before the TGD reached 60 m in September 2006, which was 20 m higher than that before the first impoundment in June 2003 [14]. Furthermore, the deposition has been persistently increasing. From 2003 to 2017, the total amount of accumulated sediment deposition in the near-dam area was 1.594 × 10⁸ m³; among them, the channel below 90 m elevation accounts for 73% of the total siltation. The average deposition thickness in the deep channel before the dam is 33.5 m, with a maximum deposition thickness of 63.7 m [15]. Of all the deposits, 0.47 × 10⁸ m³ was deposited during the period of 2008–2013, with a maximum thickness of 16.8 m [3]. Until recently, the elevation of the deposits was lower than that of the power plant inlets [16].

The rapid rise in the silt surface has greatly increased the probability of sediment entering the turbines. Hydraulic turbines may be subject to abrasion and corrosion by the heavier specific gravities of quartz and albite particles in sediment passing though, which will shorten its service life. On the other hand, the accumulation of sediment in the deep channel is conducive to the discharge of sediment from the reservoir. As the deposits have a certain fluidity (demonstrated below), the sediment can be discharged from the reservoir more easily by optimizing the dispatching measures [16].

The unique sedimentation pattern near the TGD requires a reasonable explanation based on scientific research. Models depicting the physical processes based on dynamic mechanisms are efficient research tools that would help to optimize dam operations and improve the theory of sediment transport in large reservoirs.

Many factors affect the cross-sectional deposition pattern in the reservoir, such as the sediment concentration, gradation, and lateral distributions of the flow velocity and water depth. The local river regime and cross-sectional shapes are also important. Han [17] ascribed the sedimentation pattern in which most deposits accumulate in the deep channel to the nonuniform crosswise distribution of sediment concentration and gradation. Jia et al. [18] proposed that the sedimentation pattern in the dam area of the TGR should be related to both the three-dimensional flow structure and the particle characteristics. During the early stage of reservoir impoundment, the medium size of suspended sediment before the dam was about 0.003–0.005 mm and the dry bulk density within the top 1 m of the deposits was about 0.4–1.0 t/m³ [19]. Particles deposited in the dam area are so fine that their initial dry bulk density is small, and they are able to move as a fluid. Previous simulations of the sedimentation pattern in the dam area of the TGR did not consider the fluid characteristics of the deposits, and their predictions were significantly different from the observations [11,12,20]. The differences indicate the significance of considering the fluid characteristics of the deposits. By coupling a 1D gravity-driven flow model, which depicts the motion of sediment deposits, with a traditional 1D flow and sediment transport model, the accumulation processes of fine-grained sediments in the dam area of the TGR had been preliminarily identified [21,22]. However, more work is necessary to determine the mechanisms of the fine-grained sedimentation pattern in large reservoirs.

Considering the sedimentation features in the dam area of the TGR, we proposed a physical-based numerical model for fine-grained sediment deposition in large reservoirs. A method of modeling fine-grained sediment movement that includes gravity-driven motion in reservoirs was developed and applied, along with a 3D Reynolds-averaged numerical model of water flow and sediment transport, to simulate the uncommon sedimentation pattern in the dam area of the TGR. Our method can be applied to other reservoirs with similar sedimentation problems.

2. Materials and Methods

2.1. Features of Fine-Grained Sediment Deposition in Deep Reservoirs

The Three Gorges Dam (TGD), which is located in the Three Gorges reach along the trunk of the Yangtze River in China, is a key backbone project for flood control, navigation, and water resources development and has a total installed capacity of 18.2 GW and an average annual output of 84.7 billion kWh. The dam is located in Sandouping, 40 km upstream from Yichang City, as shown in Figure 1. The normal pool level of the Three Gorges Reservoir (TGR) is 175 m (relative to the Chinese water level zero point), which corresponds to an aggregated storage capacity of 39.3 billion m³. The TGR has been in operation since June 2003. The construction of large hydraulic projects such as the TGD has generated problems in navigation, environment, and ecology, all of which are related to sediment transport processes, including advection, diffusion, deposition, and erosion.

Figure 1. Location of the TGD on the Yangtze River. (The left part of the figure is based on [23]).

Figure 1. Location of the TGD on the Yangtze River. (The left part of the figure is based on [23]).

Observations of the sedimentation pattern in the dam area of a typical section (S34#, see Figure 2 for the location) of the TGR have shown that the sediment mainly accumulated in the deep channel with a nearly horizontal surface (Figure 3) in the early years of reservoir operation. Conversely, deposits of fine-grained sediment in reservoirs typically have a relatively uniform distribution across the section. The physical and dynamic characteristics of the sediment deposits in the dam area can be used to understand this difference.

Before impoundment, the median size of the riverbed materials near the dam was approximately 0.10 mm. With the initial impoundment of the TGR, the water level rose to 135 m, and large amounts of siltation occurred in the river reach before the dam. In October 2003, the observed median particle size of the bed materials had decreased significantly to 0.02 mm. With further deposition, the median size had further decreased to 0.008 mm in May 2006. After the 156 m water storage operation of the TGR, the observed median size of the bed materials in the dam area was 0.005 mm in December 2006 and October 2007 [10]. The silt particles deposited in the dam area are extremely fine and require a long time to compact. Before compaction, the low-density deposits (the bulk density is approximately 1.00–1.20 t/m³) in the form of sludge are called fluid mud. The fluid mud has significant fluidity and can flow or deform under its own weight. If the surface slope exceeds a certain value, the fluid mud will move down the slope. This may be the main reason that most of the deposits were found in the deep channel of the TGR and that the surface of the sludge is horizontal.

Therefore, the sedimentation process of fine sediment in deep areas of reservoirs can be summarized as several physical processes, including the dynamic process of sediment transport and settlement in flowing water, gravity-driven fluid flow of deposits, and compaction and consolidation.

The gravity-driven fluid flow of fine deposits is similar to that of fluid mud. Suppose that the surface slope of the fluid mud layer is $J_c$. If $J_c>{\tau }_B/\left({\gamma }_s\mathsf{\Delta }z\right)$ (where ${\tau }_B$ is the Bingham ultimate shear force, and ${\gamma }_s$ and $\mathsf{\Delta }z$ are the weight and thickness of the fluid mud layer, respectively), the fluid mud would move down the slope under its own weight, as shown in Figure 4. Field observations have shown that the deposits in the deep channel are thick and that their surface is nearly horizontal in the dam area of the TGR (as shown in Figure 3). The main reason for this phenomenon is that the suspended sediment that was deposited during the early period of the TGR impoundment was so fine (median size of 0.004–0.01 mm) that the resulting deposits had a low initial dry bulk density and could move like fluid mud. The mud-like deposits moved downward of the surface slope due to gravity and accumulated in the deep channel of the riverbed, forming an uncommon sedimentation pattern in the early stage of the impoundment of the large reservoir.

2.2. Numerical Simulation Method for Transport and Settlement of Fine-Grained Sediment

Based on the analysis presented above, the whole processes of fine-grained sediment in the deep reservoir can be summarized as transport and settlement, gravity-driven fluid flow, and consolidation. The 3D Reynolds-averaged numerical model that depicts flow, sediment transport, and sediment settlement in a reservoir can be found in the related literature [20]. Here, we briefly describe the method.

The governing equations for water flow are given as:

$$\frac{\partial u_i}{\partial x_i}=0$$

(1)

$$\frac{\partial u_i}{\partial t}+\frac{\partial \left(u_i{u}_j\right)}{\partial x_j}=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+\frac{1}{\rho }\frac{\partial {\tau }_{ij}}{\partial x_j}+f_i$$

(2)

where $u_i$ (i = 1,2,3) are velocity components; $f_i$ is the gravitational force per unit volume; $\rho$ is the density of water; and $p$ is the pressure. In this model, turbulence stresses ${\tau }_{ij}$ are calculated with the standard $k-\epsilon$ turbulence model.

The 3D nonuniform suspended sediment transport equation can be expressed in Cartesian coordinates as:

$$\frac{\partial s_k}{\partial t}+\frac{\partial }{\partial x}\left(u{s}_k\right)+\frac{\partial }{\partial y}\left(v{s}_k\right)+\frac{\partial }{\partial z}\left(w{s}_k\right) =\frac{\partial }{\partial x}\left({\epsilon }_s\frac{\partial s_k}{\partial x}\right)+\frac{\partial }{\partial y}\left({\epsilon }_s\frac{\partial s_k}{\partial y}\right)+\frac{\partial }{\partial z}\left({\epsilon }_s\frac{\partial s_k}{\partial z}\right)+\frac{\partial }{\partial z}\left({\omega }_{sk}{s}_k\right)$$

(3)

where $u$, $v$, and $w$ are the flow velocities in the $x$, $y$, and $z$ directions, respectively; $s_k$ is the sediment concentration of group $k$; ${\omega }_{sk}$ is the settling velocity of group $k$; and ${\epsilon }_s$ is the sediment diffusion coefficient.

The equation of bed evolution caused by suspended sediment transport (deposition) is

$${\gamma }_s^{\prime }\frac{\partial Z_b}{\partial t}={\displaystyle \sum }_1^{N_s}{\omega }_{sk}\left(s_{kb}-s_{kb*}\right)$$

(4)

where ${\gamma }_s^{\prime }$ is the dry bulk density of the deposits; $s_{kb}$ and $s_{k{b}^{*}}$ are the real and equilibrium sediment concentrations of group $k$ near the bed, respectively; $Z_b$ is the bed level; and $N_s$ is the total number of sediment groups.

2.3. Numerical Simulation Method for Gravity-Driven Motion of Fine-Grained Sediment Deposits

Solid consolidation can be determined using existing methods. Here, we mainly discuss the method of simulating the gravity-driven flow of the fine-grained sediment deposits. Newly deposited fine-grained sediment deposits are a non-Newtonian fluid and are generally considered to be a Bingham fluid. Before movement, it is necessary to overcome the Bingham ultimate shear force. In this paper, a critical slope is used as a criterion for the destabilizing flow of fine-grained sediment deposits. The shallow water flow model is used to describe the gravity-driven flow process of the destabilized deposits.

When the surface slope of fine-grained sediment deposits $J_c$ satisfies the condition $J_c>{\tau }_B/\left({\gamma }_s\mathsf{\Delta }z\right)$ (where ${\tau }_B$ is the Bingham ultimate shear force, and ${\gamma }_s$ and $\mathsf{\Delta }z$ are the weight and thickness of the fluid mud layer, respectively), the deposits will move down the slope under their own weight [24]. Fine-grained deposits with different bulk densities have different Bingham ultimate shear forces. Here, we apply the formula for the Bingham ultimate shear force on fluid mud on a muddy shore [25],

$${\tau }_B=\mathsf{\alpha } \times \exp \left(\beta \times {\gamma }_s\right)$$

(5)

where $\mathsf{\alpha }$ and $\beta$ are empirical parameters; $\mathsf{\alpha }$ = 9 × 10⁻⁹; and $\beta$ = 16.719 based on previous research [25].

After losing stability, the fine-grained deposit will flow. As the thickness of the deposits is small compared to the horizontal scale of the reservoir, we assume that the flow of the fine-grained deposit can be described by the 2D shallow flow equations.

Continuity equation:

$$\frac{\partial Z}{\partial t}+\frac{\partial M}{\partial x}+\frac{\partial N}{\partial y}=q$$

(6)

Momentum equation:

$$\frac{\partial M}{\partial t}+\frac{\partial uM}{\partial x}+\frac{\partial vM}{\partial y}=-gh\frac{\partial Z}{\partial x}+D\left(\frac{{\partial }^{2}M}{\partial x^2}+\frac{{\partial }^{2}M}{\partial y^2}\right)-\frac{g{n}^{2}M\sqrt{u^2+v^2}}{h^{\frac{4}{3}}}+q{u}_0$$

(7)

$$\frac{\partial N}{\partial t}+\frac{\partial uN}{\partial x}+\frac{\partial vN}{\partial y}=-gh\frac{\partial Z}{\partial y}+D\left(\frac{{\partial }^{2}N}{\partial x^2}+\frac{{\partial }^{2}N}{\partial y^2}\right)-\frac{g{n}^{2}N\sqrt{u^2+v^2}}{h^{\frac{4}{3}}}+q{v}_0$$

(8)

where $Z$ is the surface elevation of the fine-grained deposits; M = uh; N = vh; h is the thickness of the deposits; u and v are the velocities in the x and y directions, respectively; q is the source/sink of the fine deposit in a unit area, which can be calculated according to the instability as $q=\left(\mathsf{\Delta }z-\mathsf{\Delta }{z}_c\right)/dt$, where $\mathsf{\Delta }{z}_c$ is the maximum thickness that can remain stable, and $\mathsf{\Delta }{z}_c=\frac{{\tau }_B}{\left(J_c{\gamma }_s\right)}$; $g$ is the acceleration due to gravity; D is the kinematic viscosity; and n is the integrated Manning’s roughness. As gravity is the main driving force of the flow, the diffusion terms in Equations (7) and (8) can be neglected.

2.4. Numerical Simulation Procedure

Consistent with the 3D flow and sediment transport model, the gravity-driven flow equations of fine deposits are discretized based on the finite volume method, which has the advantage of conservation. The SIMPLEC algorithm is used to solve the problem. The solution process is shown in Figure 5 and illustrated as follows.

Step (1): Use the 3D numerical model to simulate the flow, sediment transport, and sediment deposition and to obtain the distribution of the deposition thickness $\mathsf{\Delta }z$.

Step (2): Based on the silt thickness and terrain slope, determine whether the deposits in a cell are stable or not. If $J_c>{\tau }_B/\left({\gamma }_s\mathsf{\Delta }z\right)$, the deposit is unstable and the source term q will be calculated; otherwise, the deposit is stable and q is zero.

Step (3): Solve the deposit flow equations and apply the SIMPLEC algorithm to obtain the surface elevation of the deposits after being redistributed.

Step (4): Repeat steps (1)–(3) until the end of the calculation.

3. Results

Based on the topographic, hydrological, and sediment data measured in the initial stage of the TGR impoundment, a 3D Reynolds-averaged water flow and sediment transport numerical model in the dam area of the TGR (from Miaohe to the dam, see Figure 2) was established. The simulated area was approximately 15 km long. The initial terrain for the simulation was measured in June 2003 (Figure 6). Flow and sediment data at Miaohe from June 2003 to November 2008 were used for the inflow conditions. The water level in front of the dam was used for the outflow conditions. The numbers of grid points in the longitudinal, lateral, and vertical directions were 198, 97, and 18, respectively.

3.1. Flow Field and Suspended Sediment Concentration

Figure 7 shows the calculated flow field for a flow discharge of 16,700 m³/s and a pool level of 135 m above sea level. Strong 3-D flow features can be observed in the bends near the dam area, and secondary currents are significant. The velocity vectors near the bottom point to the convex bank, whereas the surface vectors point to the concave bank. Based on the calculated flow velocities, the suspended load concentration can be calculated by solving the convection–diffusion equation. Figure 8 shows the calculated suspended load concentrations on typical cross-sections. The sediment concentration decreases much more rapidly in wider valley areas. Due to the impact of the secondary currents shown in Figure 7, the maximum concentration occurs at the convex side, whereas the minimum occurs at the concave side. Figure 9 shows a comparison of the observed and simulated velocities and suspended load concentrations on typical cross-sections (see Figure 2 for the locations of the sections) for a flow discharge of 24,600 m³/s. In general, the 3D numerical model accurately predicted the flow patterns and suspended load concentrations in the dam area of the TGR.

3.2. Features of Sedimentation Pattern

The river reach in front of the TGD is a wide valley of crystalline rocky hills with significant slopes on both banks. Fine-grained sediment settled in this area after the impoundment of the TGR. If the thickness of the deposits in a certain location exceeds a critical value, the deposits will flow down the slope under their own weight. Based on the verifications of the flow velocity and sediment concentration (Figure 9), Figure 10 shows the sediment thickness near the dam calculated by traditional 3D flow and sediment transport modeling (without considering the gravity-driven fluid flow of the deposits), and Figure 11 shows the results calculated using the proposed method. Although the total amounts of sediments are consistent, the sediment distributions are quite different. The sediment thickness calculated by the traditional method, which does not consider the fluid flow of the fine-grained deposit, is relatively uniform in map view (Figure 10). In contrast, the proposed model (Figure 11) predicts that most of the deposits will accumulate in the deep channel and that a few deposits will remain on the flat terraces on both sides, which is consistent with the observations.

Figure 12 shows comparisons between the observed and calculated sedimentation on several typical cross-sections (section locations are shown in Figure 2). The cross-sectional sedimentation pattern that considers gravity-driven motion is more similar to the observations than the pattern that does not consider the motion. This again confirms the significance of the gravity-driven flow of fine deposits.

Figure 13 compares the surface elevations along the thalweg between the observations and the results calculated with and without considering the gravity-driven deposit flow. It can be seen that the calculated results that consider the deposit flow match the observations quite well. Except the point at the entrance, the maximum error at the other points is 5.1 m, which is located at section S31# (see Figure 2 for the location). Based on measured data, the averaged and maximum deposition thickness along the thalweg are 22.22 and 57.3 m, respectively. Based on simulating results without considering the gravity-driven flow, the two values are 8.32 and 31.5 m correspondingly. The corresponding values based on simulating results that consider the gravity-driven flow are 24.43 and 56.5 m. By considering the gravity-driven flow, the averaged error of deposition thickness along the thalweg decreases from −13.9 to 2.2 m.

4. Discussion

4.1. Difference between Calculation and Measurement

Although the comparison strongly confirmed the necessity and advantage of considering the gravity-driven movement of the deposits, there was some discrepancy between the observations and simulating results.

Near the dam, as shown in Figure 12a and Figure 13, the depth of the deposit at the thalweg is greater than that of the observation. This is probably because this section is near the discharge hole, so the fluid deposits may be scoured and flushed by the flow. The effect of discharge flow was not simulated precisely in this study, which may have caused the deviation. Figure 14 shows the computed velocity distribution in front of the dam. In this work, the outlets of the power plant were generalized as a single outlet area, which would flatten the flow speed nearby. Another possible reason is the terrain flattening caused by mesh interpolation in the numerical simulation. Without terrain flattening, some of the deposit may be blocked by local terrain undulation, reducing the amount of deposits that gathered in the deep channel. It is easy to imagine that those cross-sections with wider sides, more undulating, relatively gentle side slopes would suffer greater impact from terrain flattening. As the reach of the near-dam area gradually widens from Miaohe to the dam, the terrain flattening error is generally more obvious at the sections before the dam, which can explain the relatively larger discrepancy within 5 km before the dam than that of the upstream results in Figure 13.

4.2. Shortcomings of Present Method

In this work, shallow water equations are applicable to a Newtonian fluid. For a Bingham fluid, when integrating the continuity and momentum equations with depth to obtain the 2D governing equations, the nonlinear relationship between the shear stress and strain would generate significant differences with the shallow water equations [26]. If the shallow water equations are used, the resistance term would be more complicated than that of Manning’s formula [27]. Therefore, the modeling of gravity-driven flow of fine-grained deposits should be improved further. In addition, due to the discretization, the topographic data used in the numerical simulation are not consistent with the real topography. This also contributed to the discrepancy between the simulated results and observations. However, considering the gravity-driven fluid flow of the deposits, the simulation results of the overall sedimentation trend and distribution showed good agreements with the observations.

5. Conclusions

The deposition processes of fine-grained sediment in large reservoirs can be summarized as the settling of suspended sediment, gravity-driven fluid flow of the deposits, and sludge consolidation and compaction. Based on the physical mechanisms, we analyzed the dynamic causes of the patterns of fine-grained sedimentation in the dam area of the TGR during the initial stage of impoundment. The sediment deposited in the dam area of the TGR during this stage is mostly fine-grained sediment with a grain size generally between 0.004 and 0.010 mm. These fine-grained sediment deposits have a low dry bulk density and can move downslope toward the bottom of the riverbed due to the effect of gravity. These fluid deposits finally gather in the deep channel and have a horizontal surface, which is an uncommon feature of sedimentation in large reservoirs. The rapid increase in the surface elevation of the silt in front of the dam threatens the operations of the dam.

By coupling a traditional 3D water flow and sediment transport model with a newly proposed numerical method for the gravity-driven flow of fine-grained sediment deposits, a new method was developed and used to simulate the sedimentation pattern in the dam area of the TGR. The results showed that sediment would be deposited relatively uniformly in the dam area without considering the gravity-driven flow. When the gravity-driven fluid flow of the fine-grained deposits is considered, the distribution of the sediment deposits was more consistent with the observations than the results that do not consider gravity-driven flow. The simulation results showed that the deposits mainly accumulated in the deep channel with a flat surface, which agreed well with the observations. This study revealed the key role of the gravity-driven flow of deposits in simulating the deposition of fine-grained sediment in large reservoirs.