Vertical Distribution of Suspended Sediment Concentration in the Unsaturated Jingjiang Reach, Yangtze River, China

: The Rouse formula and its variants have been widely used to describe the vertical distribution of the sediment concentration in sediment-laden ﬂows in equilibrium. Han’s formula extends the Rouse formula to non-equilibrium regimes, where the diffusive ﬂux is still assumed to be Fickian. The turbulent ﬂow and suspension regimes downstream of a mega-reservoir, e.g., the Three Gorges Reservoir, usually exhibit fractal and unsaturated properties, respectively. To characterize the non-Fickian dynamics of suspended sediment and the non-equilibrium regime in natural dammed rivers, this study proposes a new formula for the concentration proﬁle of unsaturated sediment based on the Hausdorff fractal derivative advection–dispersion equation. In addition, we ﬁnd that the order of the Hausdorff fractal derivative is related to the sizes of the sediment and the degrees of non-equilibrium. Compared to Rouse and Han’s formulae, the new formula performs better in describing the sediment concentration proﬁles in the Jingjiang Reach, approximately 100 km below the Three Gorges Dam.


Introduction
Stream channels transport water and sediment from the upper watersheds. With the rapid socio-economic development in the middle-and lower-Yangtze River Basin, China, the authorities have constructed a series of mega-reservoirs in the upper reach. Among them, the Three Gorges Reservoir (TGR) may be the most famous one. Upstream damming regulates runoff and traps sediment, therefore significantly altering the natural flow and sediment regimes in the lower reaches [1,2], resulting in channel degradation, bed armoring, bar reshaping, streambank erosion [3][4][5], and eco-environmental issues [6,7].
Under the long-term influence of low-sediment-concentration flows, the riverbed will suffer from scouring, undercutting, or bank widening, and river patterns may change, which in turn affects human beings [8,9]. Before the impoundment of the TGD in 2003, the middle Yangtze River was considered a quasi-equilibrium channel with limited bed deformation and fluctuations in the suspended sediment concentration (SSC) [10]. With the completion and operation of a series of reservoirs on the main stem and tributaries of the upper Yangtze River, the lower reaches, e.g., the Jingjiang (JJ) Reach, have experienced a significant change in sediment transport regimes, mainly from quasi-equilibrium to unsaturation [9,[11][12][13]; this change has resulted in bed scouring [14], bank collapse [4], and a decrease in the SSC [15].
Since the Schmidt diffusion theory was introduced in 1925, the vertical distribution of suspended sediment has been studied for a long time by many scholars based on various theories [16]. Notably, the Rouse equation has been widely used [17], which motivated various improvements by a number of scholars. Based on laboratory data, Van Rijn and Wang and Qian [18,19] discussed the value range of the correction parameter β, a coefficient related to the diffusion of sediment particles. Umeyaina [20] studied the vertical distribution of suspended sediment in uniform open-channel flows based on the mixing length theory. Cao et al. [21] developed a model for the vertical distribution of the SSC in open-channel flows by using the turbulent burst theory. Fu et al. [22,23] quantified the suspended sediment in solid-liquid two-phase flows with a kinetic theory-based model. Zhong et al. [24] derived a suspended sediment transport equation based on a solid-liquid two-phase flow model. All of the above studies focused on the vertical distribution of the SSC in an equilibrium state, but dammed river channels are typically not in equilibrium. Dou [25] proposed a preliminary theoretical system for non-equilibrium sediment transport. Han et al. [26] modified the Rouse formula by adding a non-equilibrium coefficient c, which represents the deviation of the vertical profiles of the SSC from the Rouse result. Ma et al. [27] improved Han's equation by reducing the number of unknown parameters from three to two.
One of the implicit hypotheses of most previous studies was Fick's first law, in which the diffusive flux is proportional to the gradient of concentration. To address the anomalous diffusion in suspension, e.g., non-local vertical jumps induced by turbulent bursts, Chen et al. [28] proposed a fractional dispersion equation for sediment suspension. Recently, researchers employed fractal derivatives to study complex fluid flows and solute transport [29][30][31][32][33]. Nie et al. applied both fractional and fractal derivatives to describe the anomalous diffusion of suspended sediment in many classic laboratory experiments [34][35][36]. Sun et al. [16] proposed a Hausdorff fractal derivative-based advection-dispersion equation (HADE) model, where the Fickian diffusive flux in the Rouse formula was replaced with a fractal derivative, rescaled using a constant diffusivity. This paper aims to extend Sun et al.'s [16] work to non-equilibrium sediment regimes and develop a novel equation that can be applied to an unsaturated reach below reservoirs, such as the Jingjiang Reach, which is approximately 100 km below the Three Gorges Dam.

Study Area
The Yangtze River stands out as the third longest river in the world, and the Jingjiang (JJ) Reach (length of 347.2 km from the Zhicheng to Chenglingji gauges) is located in the middle Yangtze River. This study focuses on the sediment suspension regime in the JJ Reach from the Zhicheng (ZC) to Jianli (JL) gauges, which is approximately 243 km in length, as shown in Figure 1. ZC is approximately 100 km below the Three Gorges Dam (TGD). In addition to ZC and JL, the Shashi (SS) gauge is present in this reach and is approximately 92 km below ZC. Three tributary channels, namely, Songzikou, Taipingkou, and Ouchikou, connect the Yangtze River with Dongting Lake. Ouchikou divides the JJ Reach into upper and lower segments, and the riverbeds in these segments are composed of gravel-sand and sand, respectively [37,38].
The TGR began operation after June 2003. The initial storage level was 135 m during dry seasons, then increased to 156 m in October 2006 and 175 m in November 2008 [39]. The authorities kept building cascade dams upstream, including the Wudongde, Baihetan, Xiluodu, and Xiangjiaba reservoirs, which began operation after 2012 [1].

Data
Field survey data from the JJ Reach were provided by the Bureau of Hydrology, Changjiang Water Resources Commission, and they include flow discharge, flow velocity, water depth, sediment concentration, and size gradation of suspended sediment on vertical profiles obtained at the ZC, SS, and JL gauges. The duration of the measurement period varied at the three gauges, as shown in Table 1

Data
Field survey data from the JJ Reach were provided by the Bureau of Hydrology, Changjiang Water Resources Commission, and they include flow discharge, flow velocity, water depth, sediment concentration, and size gradation of suspended sediment on vertical profiles obtained at the ZC, SS, and JL gauges. The duration of the measurement period varied at the three gauges, as shown in Table 1

Hausdorff Fractal Space
Many previous investigations of sediment suspension in steady and equilibrium sediment-laden flows are based on the Schmidt equation [16,28]: where ω denotes the sediment settling velocity (L/T); S is the volumetric concentration of sediment; y represents the vertical distance from the streambed [L]; and ε sy is the vertical component of the sediment diffusion coefficient [L 2 /T]. Equation (1) indicates that the vertical distribution of the SSC in the equilibrium regime is determined by the balance between two opposite movements: downward settling and upward diffusion. By assuming ε sy = ε m = κu * (1 − y/h)y where ε m is the fluid eddy viscosity; κ is the Von Karman constant; u * is the shear velocity [L/T]; and h is the water depth [L]. Rouse [17] provided an analytical solution to Equation (1): where a is a given height [L]; S * is a reference concentration at a given height a [L] above the streambed; h is the water depth [L]; κ is the Von Karman constant; and u * is the shear velocity [L/T]. ω /κu * is called the suspension index.
The Rouse equation adopts the passive scalar hypothesis and is therefore suitable for dilute suspension, such as in Jiangjiang cases. However, it apparently suffers from fatal deficiencies, such as zero and infinite values of concentration at the water surface and streambed, respectively. More crucially, the Rouse equation is not applicable to non-equilibrium states in which the downward settling of suspended sediment cannot be balanced by the upward turbulent diffusion in Equation (1). To address this issue, Han et al. [26] extended the suspension index to cω /κu * . Here c is called the non-equilibrium parameter, which is assumed to be a function of the ratio of the depth-averaged concentration S to the carrying capacity S * , namely c = 1 − f S S * [27]. c > 1 and c < 1 represent the unsaturation and over-saturation regimes, respectively. However, the parameter c actually represents the deviation of the non-equilibrium SSC profiles from that in Equation (2). Due to the limitation of Equation (2) (e.g., lack of non-Fickian diffusion consideration, as explained in Chen et al. [28] and Nie et al. [36]), Han's equation does not always accurately express the vertical profile of the SSC, sometimes even contradicting real-world measurements [41]. Figure 2a sketches a diagram of sediment suspension with the flow depth divided into eight layers. Statistically, any particles from the bottom layer have a chance to jump to any upper layers. In equations based on Fickian law, e.g., Rouse or Han's equations, it is assumed that the turbulence-induced upward movement of particles is controlled by the Central Limit Theorem. In other words, the probability of a suspended particle jumping upward decays exponentially with the jump length [28]. However, it can be argued that non-local particle jumps are common in an open-channel flow due to the existence of coherent turbulent structures. As a more general expression, the probability of a particle moving a given distance upward decreases as a power law with an order of α (0 < α ≤ 1, where α = 1 indicates that Fick's first law is valid).
where γ is called the Hausdorff space coefficient, which is related to the fractal dimension of the space.

Hausdorff Fractal Deviation for Non-Equilibrium SSC Vertical Profile
To address the anomalous diffusion and the non-equilibrium regime in the JJ Reach after the TGR impoundment, we propose a Hausdorff fractal derivative, as follows.
First, we modify the Schmidt equation, as follows: and rearrange Equation (4) as: An alternative way to ensure that α ≡ 1, i.e., Fick's first law, is always valid is to transform the Euclidean geometric space to the Hausdorff fractal space, as shown in Figure 2. The metric transformation is based on two assumptions: (a) the physical laws are invariant; and (b) the influence of anomalous fluctuations on physical behavior is equivalent to that of the fractal space transform [42]. With these assumptions, the Hausdorff fractal derivative can be expressed as where γ is called the Hausdorff space coefficient, which is related to the fractal dimension of the space.

Hausdorff Fractal Deviation for Non-Equilibrium SSC Vertical Profile
To address the anomalous diffusion and the non-equilibrium regime in the JJ Reach after the TGR impoundment, we propose a Hausdorff fractal derivative, as follows.
First, we modify the Schmidt equation, as follows: where q sy is the vertical net sediment flux at the position y [L/T]. Following Chen et al., 2013 [28], we substitute the depth-variable ε sy for a depth-averaged or depth-constant ε s : and rearrange Equation (4) as: ∂s ∂y where q ss = q sy /ε s . It is assumed that q ss linearly changes with the water depth, reaching a value of zero at the bed and a maximum at the water surface: q ss = q s (y/h) (7) where q s is constant and q s = q ss when y = h. Equation (6) is a non-homogeneous linear first-order differential equation, and its general solution can be expressed as: y a q s y h e ω εs (y−a) + C where C is an integral constant. We set: where D = q s e We express e ω εs (y) = κ · e ω εs (a) using first-order Tayler's formula and κ = 1 + ω ε s (y − a). Since the suspended sediment in the JJ Reach is mainly composed of wash-load with particle size d < 0.062 mm [43], therefore ω << ε s and κ ≈ 1.
By substituting Equation (11) into Equation (8), we obtain: According to the boundary condition, s = s a when y = a; therefore, C = 1. Equation (12) becomes: Equation (13) describes the vertical distribution of the non-equilibrium SSC, and its Hausdorff fractal derivative is: where γ is the Harsdorf fractal coefficient and 0 <γ ≤1.

Comparison of Formulae in Describing Vertical Sediment Profiles
Numerous formulae can be used when simulating the vertical distribution of the SSC in open channel flows. In this study we validate the applicability of Equation (14) by comparing the calculations to the measurements with various flow discharges and particle sizes at the ZC, SS, and JL gauges in the JJ reach.
In Figure 3, we compare the performance of three formulae, namely, Rouse equation, Han's equation, and Equation (14), in simulating the real-world vertical distribution of the SSC in the JJ Reach. The measured data in Figure 3a,b are from the ZC gauge; the data in  (14) performs the best in simulating the vertical distribution of suspended sediment under a non-equilibrium sediment regime. Due to the relatively high sediment concentration in the near-bed area and the low sediment concentration in the upper water layers, there is a remarkable "tailing phenomenon" in the near-bed zone (<0.2 h) of the vertical profiles in non-equilibrium sediment regimes, for which the Rouse formula is not at all applicable. Han's formula can better describe the tails; however, it sometimes underestimates the sediment concentration in the upper water layers, especially for finer particles in the near-surface zone (>0.8 h). This limitation is very likely because Han's formula neglects the anomalous diffusion of suspended sediment. As indicated by Chen et al. [28], turbulent bursts may sweep the bed, grab sediment particles, and non-locally jump to the upper water layers. The orders of the Hausdorff fractal derivative in Figure 3 are as low as 0.1-0.3, indicating that the upward anomalous diffusion ability in the non-equilibrium regime might be stronger than that quasi-equilibrium conditions.

Parameters That Influence γ
We analyzed the parameters that can affect the order of the Hausdorff fractal derivative γ. Figure 4 shows the relationship between γ and the particle size d.

Parameters That Influence γ
We analyzed the parameters that can affect the order of the Hausdorff fractal de ative γ. Figure 4 shows the relationship between γ and the particle size d.  The fitting curves in Figure 4 indicate that the order of the Hausdorff fractal der tive γ increases with particle size d in all four periods. The small γ values for fine parti indicate that these particles are more likely to make longer vertical jumps in turbu bursts, resulting in super-diffusion in Euclidean geometric space, which needs to be r cally transformed into Hausdorff fractal space to validate Fick's first law. In contrast, coarser particles are more controlled by gravitational settling and less affected by tur lent bursts.
In addition to particle size d, the degree of non-equilibrium should also be an portant parameter that influences the vertical profiles, as discussed above. In this stu The fitting curves in Figure 4 indicate that the order of the Hausdorff fractal derivative γ increases with particle size d in all four periods. The small γ values for fine particles indicate that these particles are more likely to make longer vertical jumps in turbulent bursts, resulting in super-diffusion in Euclidean geometric space, which needs to be radically transformed into Hausdorff fractal space to validate Fick's first law. In contrast, the coarser particles are more controlled by gravitational settling and less affected by turbulent bursts.
In addition to particle size d, the degree of non-equilibrium should also be an important parameter that influences the vertical profiles, as discussed above. In this study, we calculated the degree of non-equilibrium ξ using the following formula: where S is the depth-averaged concentration and S * is the depth-averaged concentration before the impoundment of the TGD, which is equivalent to the carrying capacity under equilibrium conditions. According to the long-term data in the pre-TGD period, the JJ Reach was in a quasi-equilibrium state, i.e., the ξ value was close to 0. Negative and positive ξ values indicate unsaturation and supersaturation regimes, respectively. The calculation of the carrying capacity is based on Zhang's formula [44]. Figure 5 indicates that the order of the Hausdorff fractal derivative γ decreases as the absolute value of the degree of non-equilibrium increases, approaching −1. Figure 5 also indicates that the higher non-equilibrium degrees are usually associated with smaller particle size d. The smallest value of γ occurred when ξ = −0.973 and d = 0.028 mm. Again, this trend may be related to the enhanced non-local particle jumps in the non-equilibrium regime. where S is the depth-averaged concentration and * S is the depth-averaged concentration before the impoundment of the TGD, which is equivalent to the carrying capacity under equilibrium conditions. According to the long-term data in the pre-TGD period, the JJ Reach was in a quasi-equilibrium state, i.e., the ξ value was close to 0. Negative and positive ξ values indicate unsaturation and supersaturation regimes, respectively. The calculation of the carrying capacity is based on Zhang's formula [44]. Figure 5 indicates that the order of the Hausdorff fractal derivative γ decreases as the absolute value of the degree of non-equilibrium increases, approaching −1. Figure 5 also indicates that the higher non-equilibrium degrees are usually associated with smaller particle size d. The smallest value of γ occurred when ξ = −0.973 and d = 0.028 mm. Again, this trend may be related to the enhanced non-local particle jumps in the non-equilibrium regime. Equation (16) provides the best linear fit:

Sensitivity Analysis
In Figure 6, we evaluated the sensitivity of the vertical distribution of the SSC to the selection of γ values. The case chosen here is the field measurement from the ZC gauge on 12 July 2006, for which d 50 = 0.0038 mm. The best fit in Figure 6 is obtained with γ = 0.3. Figure 6 indicates that a small γ, e.g., 0.1, in Equation (14) will lead to underestimation of the concentration in the near-bed zone and overestimation in the upper water layers. In contrast, if γ is too large, e.g., 0.6, in Equation (14), the concentration in the near-bed zone will be overestimated, with underestimation in the upper water layers. The estimated concentration values at 0.8 h will decrease to 36% and 88% of the measured values (0.255 kg/m 3 ) for γ = 0.6 and γ = 0.4, respectively. This indicates that the calculated profiles for the fine particles are sensitive to γ values.

Sensitivity Analysis
In Figure 6, we evaluated the sensitivity of the vertical distribution of the SSC to the selection of γ values. The case chosen here is the field measurement from the ZC gauge on 12 July 2006, for which d50 = 0.0038 mm. The best fit in Figure 6 is obtained with γ = 0.3. Figure 6 indicates that a small γ, e.g., 0.1, in Equation (14) will lead to underestimation of the concentration in the near-bed zone and overestimation in the upper water layers. In contrast, if γ is too large, e.g., 0.6, in Equation (14), the concentration in the near-bed zone will be overestimated, with underestimation in the upper water layers. The estimated concentration values at 0.8 h will decrease to 36% and 88% of the measured values (0.255 kg/m 3 ) for γ = 0.6 and γ = 0.4, respectively. This indicates that the calculated profiles for the fine particles are sensitive to γ values. When the downward settling of the suspended sediment cannot be balanced by upward turbulent diffusion, a net sediment flux qs will be generated in the vertical direction. Figure 7 shows that the vertical distribution of the SSC is also sensitive to the vertical net sediment flux qs. When the absolute values of qs approach zero, the sediment regime is in quasi-equilibrium, the concentration in the near-bed zone is underestimated, and the concentration in the upper water layers is overestimated. When qs = −0.0005, Equation (14) becomes nearly identical to the Rouse formula. In Figure 6, qs = −0.0005 provides the best fit to the measurements at both 0.8 h and 1.0 h. However, when qs is taken as −0.0006 and −0.0008, the estimated concentration values at 0.8 h are reduced to 94% and 86% of the measured value (0.255 kg/m 3 ), respectively. When the downward settling of the suspended sediment cannot be balanced by upward turbulent diffusion, a net sediment flux q s will be generated in the vertical direction. Figure 7 shows that the vertical distribution of the SSC is also sensitive to the vertical net sediment flux q s . When the absolute values of q s approach zero, the sediment regime is in quasi-equilibrium, the concentration in the near-bed zone is underestimated, and the concentration in the upper water layers is overestimated. When q s = −0.0005, Equation (14) becomes nearly identical to the Rouse formula. In Figure 6, q s = −0.0005 provides the best fit to the measurements at both 0.8 h and 1.0 h. However, when q s is taken as −0.0006 and −0.0008, the estimated concentration values at 0.8 h are reduced to 94% and 86% of the measured value (0.255 kg/m 3 ), respectively. Fractal Fract. 2023, 7, x FOR PEER REVIEW 12 of 14 Figure 7. The SSC profiles with various vertical net sediment flux qs (calculated using Equation (13)). The other parameters are κ = 0.4, εs = 0.13, ω = 1.45 × 10 −5 , u* = 0.105, γ = 0.3.

Conclusions
Due to the operations of a series of mega-reservoirs in the upper reachof the Yangtze River, the sediment regime in the JJ Reach has changed dramatically, displaying fractal and unsaturated characteristics. Based on the Hausdorff fractal derivative, which transforms the Euclidean geometric space to the fractal space, we propose a new equation to describe the vertical distribution of the SSC in the JJ Reach. The new equation performs better than the Rouse and Han's formulae, avoiding their boundary-value deficiencies. In addition, we find that the order of the Hausdorff fractal derivative γ is related to the sizes of sediment particles and the degrees of non-equilibrium. We still need further data from field measurements and experiments to obtain a more robust prediction of γ in the future.  (13)). The other parameters are κ = 0.4, ε s = 0.13, ω = 1.45 × 10 −5 , u * = 0.105, γ = 0.3.

Conclusions
Due to the operations of a series of mega-reservoirs in the upper reachof the Yangtze River, the sediment regime in the JJ Reach has changed dramatically, displaying fractal and unsaturated characteristics. Based on the Hausdorff fractal derivative, which transforms the Euclidean geometric space to the fractal space, we propose a new equation to describe the vertical distribution of the SSC in the JJ Reach. The new equation performs better than the Rouse and Han's formulae, avoiding their boundary-value deficiencies. In addition, we find that the order of the Hausdorff fractal derivative γ is related to the sizes of sediment particles and the degrees of non-equilibrium. We still need further data from field measurements and experiments to obtain a more robust prediction of γ in the future.