Operating Effects of the Three Gorges Reservoir on the Riverbed Stability in the Wuhan Reach of the Yangtze River

Abstract

Reservoir construction may modify the downstream flow and sediment transport, and correspondingly result in adjustments in morphodynamics of a river, especially riverbed instability. Based on hydrological datasets from 2003 to 2019 during the post-dam period using the topographic data of 57 fixed cross-sections in the Wuhan reach of the Yangtze River, we calculated the indexes representing the channel stability. Moreover, considering the effects of flow, sediment concentration, grain size of sediment, and water depth, we propose a method for calculating the equilibrium values of cumulative erosion and the lateral migrate intensity of thalweg in this paper, and the method combines with the delayed response model (DRM) to comprehensively analyze the variations in the longitudinal and lateral stabilities of the riverbed. The results revealed that the channel has been obviously eroded in the downstream reach, resulting in a 76% decrease in sediment discharge after the impoundment of the Three Gorges Reservoir (TGR). Specifically, in the past 17 years, the cumulative erosion in the Wuhan reach of the Yangtze River reached 1.72 m, while the bankfull depth increased by 1.87 m. The lateral migrate intensity of thalweg increased in response to Coriolis force, with an increase of 22.3%. Taken together, the results show that the proposed formula can effectively simulate the variation process of channel stability, and it also quantifies the extent of the influence weight of interannual flow and sediment regimes. The morphodynamics adjustments in the channel stability of our studied reach were closely related to the previous five-year flow and sediment regimes, implying that channel evolution may lag behind the changes in flow and sediment discharges.

1. Introduction

Construction of the reservoirs may alter the downstream natural hydrological process from both the incoming flow and sediment regimes, such as a reduction in the peak flow and increased discharge in dry seasons, as well as decreased amplitude in annual and interannual variations of runoff, decreased sediment discharge, drastic reduction in the downstream sediment concentration, and decreased sediment particle size [1]. The variations in the characteristics of flow and sediment transport can lead to adjustments in channel morphology, thus resulting in riverbed instability [2,3]. Riverbed deformation is mainly divided into longitudinal and lateral variations [4,5]. The study of Lyu et al. [6] showed that the downstream riverbed of the Yangtze River has been scoured and cut-down to varying degrees, and that the thalweg has generally been cut-down after the impoundment of the TGR. Gierszewski et al. [7] found that the downstream channel of Wloclawek Dam has been scoured, and such scour has been inhibited, the river bank erosion has intensified, and the channel has widened since the roughening of the riverbed. It can be seen that the riverbed stability of the downstream channel is closely related to the reservoirs, and it is necessary to adopt systematic evaluation methods to study the influence mechanism of reservoir operation on downstream riverbed stability.

In addition, in the process of riverbed stability adjustments reaching equilibrium in the lower reaches of the reservoir, when upstream hydrological conditions and other external conditions change, rivers cannot quickly respond to the changed conditions, and they could reach a new equilibrium state by means of scouring and silting changes or other morphological adjustments [8], which is referred to as a delayed response. The weighted averaging and moving averaging methods, as well as the delayed response model (DRM) recently developed and continuously updated, are adopted to describe this delayed response [9]. Liang et al. [10] discovered that the previous sectional morphology had an effect on channel geometry in the Lower Yellow River using the weighted averaging method, and then drew a conclusion that the fluvial system had a “memory effect”. Moreover, Xia et al. [11] used the moving averaging method to build an empirical relationship between bankfull discharge and the average incoming sediment coefficient (ratio of sediment to flow discharge) during the flood season in the Lower Yellow River. In contrast, inspired by the rate law of Graf [12], Wu et al. introduced the concept of river geomorphology to illustrate a point that the self-adjustment rate of a channel is proportional to the differences between its current and equilibrium values. Then, the DRM was developed, and has been widely applied to the Yellow River Basin and some other rivers [13,14,15,16].

The Three Gorges Project (TGP) is the world’s largest integrated hydro-project. Since its operation in 2003, many scholars have carried out studies on the morphological self-adjustments and the delayed response for the Jingjiang reach downstream of the dam. For instance, Zhang et al. [17] took the sectional area corresponding to the annual dry season water level of the Shashi reach as the calculation index and then introduced the concept of “annual combined flow”, and their model used the piecewise function iteration method, with a correlation coefficient (R²) of 0.7–0.8. Moreover, Lyu et al. [18] introduced the sediment transport capacity factor and considered the impacts of sediment discharge in the flood season, to establish a delayed response model for calculating the channel sediment discharge, with a correlation coefficient (R²) greater than 0.93. In addition, the DRM has also achieved good results in simulating the sediment deposition in the TGR [19], the change of shoal area in the Jingjiang reach, and so on [20].

The previous analysis of riverbed evolution was mostly based on the research on the channel geometry of alluvial rivers reaching equilibrium [21,22,23,24,25,26,27,28,29]. For instance, as for the influence factors of thalweg migration, it is generally believed that the thalweg migration intensity increases with the increase in sediment concentration and channel slope [30]. In contrast, it is believed that for meandering rivers, when the annual flow discharge suddenly changes dramatically, the migration intensity will also increase [31]. Moreover, Wyżga et al. [32] studied the water level changes under a certain flow according to the water level–flow relationship, and then explored the influence of riverbed degradation and incision on the lateral channel migration. However, due to the fact that morphological evolution cannot occur freely in some rivers based only on natural factors, because there are some strong external constraints including human activities, equilibrium-based flow and sediment regimes may not be fully applicable to the study of the riverbed evolution characteristics. Furthermore, there are not many studies on the systematic evaluation methods and parameters of riverbed stability. Therefore, it is meaningful and important to exactly explain the variation characteristics of non-equilibrium riverbed stability and the delayed response phenomenon with the existing equilibrium theories of riverbed evolution.

The current study was conducted mainly for the following purposes: (1) to calculate the longitudinal and lateral stability indexes of the riverbed in the studied reach from 2003 to 2019, and to analyze the external factors for their variations, (2) to develop a method for estimating the equilibrium values of the representative indexes of the longitudinal and lateral stability variations of the riverbed after an analysis of several external factors, and (3) by establishing a quantitative relationship between the current values and the equilibrium ones, to analyze the variation characteristics of the longitudinal and lateral stability of the riverbed in the studied reach after the impoundment of the TGR.

2. Study Area

The study area in this work is located in the Wuhan reach of the Middle Yangtze River, which is approximately 264 km upstream of the Poyang Lake estuary in the Middle Yangtze River. It extends from Shamao Mountain in the Hannan District of Wuhan city to Yangluo Town (electric tower) in Xinzhou District, spanning a total length of approximately 70.3 km (Figure 1). The channel mainstream flows from Shamao Mountain to Zhuankou, then to the left branch of Baishazhou, before flowing down after passing through Sheshan. Then, it flows through the right branch of Tianxingzhou and continues to flow down along Yangluo on the left bank after converging with the water from the left branch near the water gap at the tail of Tianxingzhou sandbar.

As shown in Figure 1 and

Table 1. Essential information of the two sub-reaches in the Wuhan reach.

Channel Reach	Start and End Points	Channel Type	Channel Length (km)	Maximal/Minimal Spacing between Two Sections (km)
Upper sub-reach	Shamao Mountain–Sheshan	Straight braided channel	35.0	2.32/0.62
Lower sub-reach	Sheshan–Yangluo	Slightly curved braided channel	35.3	2.82/0.83

, the Wuhan reach, extending from the cross-section CZ51 to CZ57, covers 57 cross-sections surveyed by the Changjiang Water Resources Commission (CWRC) after the impoundment of the TGR. The Wuhan reach is divided into upper and lower sub-reaches. The upper sub-reach is the straight braided channel of Baishazhou sandbar extending from CZ51 to HLZ09+1, with its main branch being the left branch of Baishazhou, while the lower sub-reach extends from HLZ09 + 1A to CZ57, and is a slightly curved braided channel of Tianxingzhou sandbar with its main branch being the right branch of Tianxingzhou. Wuhanguan is located 1.5 km below the Hanjiang River estuary in the Wuhan reach (Figure 1), where Hankou Hydrological Station is located. The number of years of hydrological data collection of this station is higher, and these data can be used as the incoming flow and sediment regimes for studying the Wuhan reach in this work.

From 2003 to 2007, in order to combat possible catastrophic floods, most areas of the Wuhan reach completed various bank revetment projects, including river beach flood control and comprehensive improvement projects (Figure 1), which can effectively ensure the stability of channel width and have basically controlled the large-scale bank erosion.

The main sources of runoff and sediment in the Wuhan reach are the main stream and tributaries in the Upper Yangtze River. Since the dam closure began in 2003, significant changes have occurred in the downstream flow and sediment regimes, and the channel geometry has been evolving constantly in the Wuhan reach, as well as the runoff and sediment of Hankou station have been affected. To investigate the variation characteristics of the studied reach, we collected the post-flood cross-sectional profiles of 57 specified locations and the hydrological data of Hankou Hydrological Station.

Table 2. Sources of measurements.

Data Type	Station (Location)	Period of Record	Sources
Monthly discharge	Hankou Hydrological Station	1954–2019	CWRC
Monthly water level	Hankou Hydrological Station	1954–2019	CWRC
Monthly sediment concentration	Hankou Hydrological Station	1954–2019	CWRC
Monthly suspended sediment gradation	Hankou Hydrological Station	1954–2019	CWRC
Surveyed profiles	57 fixed cross-sections within the Wuhan reach	2003–2019	CWRC

shows the data sources, measurement time, and other specific information.

3. Methods

3.1. Methods for Evaluating Channel Stability

To analyze the channel evolution, various indexes of bankfull channel morphology (Figure 2) were calculated to reflect the geometric characteristics of the channel cross-sections. These indexes can directly reflect the relevant changes in the channel caused by the variations of incoming flow and sediment regimes. The variations in bankfull width and depth can show the variation characteristics of the lateral and longitudinal stability of the riverbed.

The total length of the Wuhan reach in the Yangtze River in this study is 70.3 km, and the channel morphology of each typical cross-section is different (Figure 2). The variation characteristics of the channel morphology of a certain cross-section cannot represent the variation characteristics of the whole reach. Therefore, a study of the variations of the whole reach is necessary. The method combining the geometric average with the weighted average of the cross-section spacing based on logarithmic transformation was used in this paper [33]. Assuming that the measured number of cross-sections of a reach is N, the index of the overall channel morphology in the whole reach can be expressed as:

$$\overline{{\mathrm{G}}_{\mathrm{bf}}}=\exp \{\frac{1}{2\mathrm{L}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}-1}\left(\ln {\mathrm{G}}_{\mathrm{bf}}{}^{\mathrm{i}+1}+\ln {\mathrm{G}}_{\mathrm{bf}}{}^{\mathrm{i}}\right)\times \Delta {\mathrm{x}}_{\mathrm{i}}\}$$

(1)

where, $\overline{{\mathrm{G}}_{\mathrm{bf}}}$ is the average channel morphology index (bankfull width, $\overline{{\mathrm{B}}_{\mathrm{bf}}}$, depth, $\overline{{\mathrm{H}}_{\mathrm{bf}}}$, etc.),${\mathrm{G}}_{\mathrm{bf}}$ is the channel morphology index of the ith section, L is the total length of the studied reach, and ∆x is the distance between two adjacent sections (i, i+1).

The water level of each cross-section was estimated using the water surface line. By viewing the multi-year water level–discharge curve at Hankou Hydrological Station, the channel corresponding to the water level at the flow rate of 45,000 m³/s was taken as the bankfull channel, and the average bankfull channel morphology index was calculated according to the water level of the cross-sections. The corresponding bankfull level of the most upstream and downstream cross-sections was 21.48 and 20.42 m, respectively.

The external condition of the bankfull level should be considered in the calculation of bankfull channel indexes. Two indexes, including cumulative erosion and lateral migrate intensity of thalweg (in this paper we called it the distance from the thalweg to the initial point), were also considered in this work. Their variations not only fully reflect the variation characteristics of lateral and longitudinal stability, but also avoid the restriction of external conditions.

The riverbed cumulative erosion was calculated according to the measured data of the fixed cross-sections in the studied reach. The mean bed level of a cross-section of previous years was calculated first. The mean bed level of the cross-sections in the initial year was used as the initial mean bed level, and the difference between the mean bed level of a cross-section of previous years and the initial value was used as the cumulative erosion of the section in the same year. As a result, each section had a cumulative erosion value every year. Since the cumulative erosion of some sections was negative in a certain year, the arithmetic mean of cumulative erosion values of fixed cross-sections of previous years was taken as the cumulative erosion of the whole studied reach.

The distance from the thalweg to the initial point was also calculated according to the measured data of fixed cross-sections in the studied reach. The distance between the thalweg point and the initial point on the left bank of each cross-section of previous years was calculated, and then the overall distance from the thalweg to the initial point of the studied reach in each year was obtained using a method similar to that used to calculate the average channel morphology index.

3.2. Methods for Calculation of the Changes in Morphological Parameters

When the upstream incoming flow and sediment regimes change, the process of scouring and silting adjustments of the whole channel to reach a new equilibrium state often lags behind the variations in these upstream incoming flow and sediment regimes, which requires a time-dependent approach. Wu [9] created a delayed response model (DRM) centered on the concept of riverbed evolution based on the variable rate principle in physics. This model has been well-applied to simulate the self-adjustment process of a river after the upstream incoming flow and sediment regimes change. The basic equation of the model is as follows:

$$\mathrm{dy}/\mathrm{dt}=\mathsf{\beta }(\mathrm{y}-{\mathrm{y}}_{\mathrm{e}})$$

(2)

where y and y_e are the actual value of the characteristic variable at t and the equilibrium value to which it tends, and β is the characteristic parameter of the riverbed evolution adjustment rate. The larger the β, the faster the adjustment rate of the characteristic variable, and vice versa.

Equation (2) is integrated, and the multi-step model is as follows:

$${\mathrm{y}}_{\mathrm{n}}=(1-{\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}}){\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{(\mathrm{e}}^{-(\mathrm{n}-\mathrm{i})\mathsf{\beta }\Delta \mathrm{t}}{{\mathrm{y}}_{\mathrm{ei}})+\mathrm{e}}^{-\mathrm{n}\mathsf{\beta }\Delta \mathrm{t}}{\mathrm{y}}_{\mathrm{e}0}$$

(3)

where y_n is the value of channel index y at the end of the nth ∆t time step, y_ei is the equilibrium value of y in the ith time step, i and n are the serial number and total number of time steps respectively, and y_e0 is the initial value of y and i = 0.

In Equation (3), the index of channel y_n represents the weighted mean of the calculated equilibrium values, y_ei, of previous years, the weight of the initial year (i = 0) is ${\mathrm{e}}^{-\mathrm{n}\mathsf{\beta }\Delta \mathrm{t}}$, and that of the ith year (i = 1, 2,…n) is ${(1-\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}}{)\mathrm{e}}^{-(\mathrm{n}-\mathrm{i})\mathsf{\beta }\Delta \mathrm{t}}$. The closer it is to the current year, the greater the corresponding weight will be. The corresponding influence weight ${(1-\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}})$ of the current year (i = n) is the largest. The sum of the influence weights of each year is equal to 1 in Equation (3), and ${(1-\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}}){\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{e}}^{-(\mathrm{n}-\mathrm{i})\mathsf{\beta }\Delta \mathrm{t}}{+\mathrm{e}}^{-\mathrm{n}\mathsf{\beta }\Delta \mathrm{t}}=1$.

Li and Wu [13] found that the influence weight, ${\mathrm{e}}^{-\mathrm{n}\mathsf{\beta }\Delta \mathrm{t}},$ of the initial year (i = 0) was greater than that of the previous year, ${(1-\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}}{)\mathrm{e}}^{-(\mathrm{n}-\mathrm{i})\mathsf{\beta }\Delta \mathrm{t}}$(i = 1), when $\mathsf{\beta }>0.693$, which is contrary to the principled view that the impact on the channel morphology decreases with time after a change in the upstream incoming flow and sediment regimes. This then improved the DRM by changing the influence weight of the initial year into ${(1-\mathrm{e}}^{-\mathsf{\beta }}{)\mathrm{e}}^{-\mathrm{n}\mathsf{\beta }}$ to eliminate the dependence of the model on the initial value, so that the influence weight of the previous flow and sediment regimes on the current reach is distributed from large to small when the time step is from near to far. The improved DRM can be expressed as follows:

$${{\mathrm{y}}_{\mathrm{n}}=(1-\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}}){\displaystyle \sum }_{\mathrm{i}=0}^{\mathrm{n}}{(\mathrm{e}}^{-(\mathrm{n}-\mathrm{i})\mathsf{\beta }\Delta \mathrm{t}}{\mathrm{y}}_{\mathrm{ei}})$$

(4)

However, the drawback of this model is that the sum of weights of each year is less than 1, and ${(1-\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}}){\displaystyle \sum }_{\mathrm{i}=0}^{\mathrm{n}}{\mathrm{e}}^{-(\mathrm{n}-\mathrm{i})\mathsf{\beta }}$ = ${1-\mathrm{e}}^{-(\mathrm{n}+\mathrm{i})\mathsf{\beta }\Delta \mathrm{t}}<1$ may make the calculated equilibrium value of this index larger. Thus, the improved weighted normalized DRM [34] can be calculated as follows:

$${\mathrm{y}}_{\mathrm{n}}=\frac{{1-\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}}}{{1-\mathrm{e}}^{-(\mathrm{n}+1)\mathsf{\beta }\Delta \mathrm{t}}}{\displaystyle \sum }_{\mathrm{i}=0}^{\mathrm{n}}{\mathrm{e}}^{-(\mathrm{n}-\mathrm{i})\mathsf{\beta }\Delta \mathrm{t}}{\mathrm{y}}_{\mathrm{ei}}$$

(5)

In this formula, the influence weights of each year, $\frac{{1-\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}}}{{1-\mathrm{e}}^{-(\mathrm{n}+1)\mathsf{\beta }\Delta \mathrm{t}}}{\mathrm{e}}^{-(\mathrm{n}-\mathrm{i})\mathsf{\beta }\Delta \mathrm{t}},$ can meet the general rule that the farther the upstream incoming flow and sediment regimes are from the current year, the smaller the influence weight on a certain characteristic evolution of the channel will be, and it can meet the requirement that the sum of the weights of each year is 1.

Assuming that each calculation index is only affected by the current year’s flow and sediment regimes (i.e., n = 0, △t = 1 year), the parameters in the formula are fitted and the calculated values are compared to the measured values. Then, the correlation coefficient (R²) between the calculated and measured values and their relative error MNE can be calculated as follows:

$${\mathrm{R}}^2=\frac{{\left[{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}({\mathrm{f}}_{\mathrm{mi}}-\overline{{\mathrm{f}}_{\mathrm{m}}}\left)\right({\mathrm{f}}_{\mathrm{ci}}-\overline{{\mathrm{f}}_{\mathrm{c}}})\right]}^2}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{({\mathrm{f}}_{\mathrm{mi}}-\overline{{\mathrm{f}}_{\mathrm{m}}})}^2{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{({\mathrm{f}}_{\mathrm{ci}}-\overline{{\mathrm{f}}_{\mathrm{c}}})}^2}$$

(6)

$$\mathrm{MNE}=\frac{1}{\mathrm{N}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}}\left|\frac{{\mathrm{f}}_{\mathrm{ci}}-{\mathrm{f}}_{\mathrm{mi}}}{{\mathrm{f}}_{\mathrm{mi}}}\right|\times 100\%$$

(7)

where f_mi and f_ci are the measured and calculated values of the characteristic variables at the ith time step, and $\overline{{\mathrm{f}}_{\mathrm{m}}}$ and $\overline{{\mathrm{f}}_{\mathrm{c}}}$ are the mean of measured and calculated values of the characteristic variables in the whole calculation period (N time steps).

In this calculation, the measured results are first substituted into the characteristic variable y_n on the left of Equation (5). Given that time step Δt = 1 year, and the number of preceding affected years, n, is gradually increased from 0, the iterative calculation of Equation (5) is carried out, and the parameters in the model formula are fitted by the simulated annealing algorithm, and then the calculated values are compared to the measured values of the characteristic variables. The correlation coefficient (R²) at different values of the preceding affected years n is calculated. n, in which the R² is stable around the maximum value, is regarded as the number of years that the calculation index of this reach is best affected by changes in the antecedent incoming flow and sediment regimes (the actual number of affected years is n + 1). In this study, the measured data from 2003 to 2016 were used to fit the parameters in the model formula, and the accuracy of the model formula was verified by the measured data from 2017 to 2019.

The cumulative riverbed erosion and distance from the thalweg to the initial point were selected as the characteristic parameters of the longitudinal and lateral stability changes of the studied reach, and the corresponding calculation methods were proposed according to the DRM. When the incoming flow and sediment regimes are known, there will be a corresponding equilibrium state in any period of time. However, due to the continuous change in the incoming flow and sediment and the delayed response of riverbed evolution, in a fixed period of time it is difficult to observe the equilibrium state in the riverbed evolution. Therefore, the key to calculating the two characteristic parameters in this study was to determine the calculation expression of their equilibrium values.

When the river is in a state of sediment transport equilibrium (the sediment discharge of the channel is equal to the sediment carrying capacity, and the channel does not change in scouring and silting), the formula of the sediment carrying capacity can be established:

$$\mathrm{S}*=\mathrm{k}{\left(\frac{{\mathrm{U}}^3}{\mathrm{gh}\mathsf{\omega }}\right)}^{\mathrm{m}}$$

(8)

where k is the coefficient, U is the velocity at a cross-section (m/s), g is the standard gravity (N/kg), h is the average water depth (m), $\mathsf{\omega }$ is the settling velocity of sediment (m/s), and m is a positive exponent.

Considering the flow continuity formula, the flow resistance formula, and the Goncharov formula (d < 0.15 mm):

$$\mathrm{Q}=\mathrm{BhU}$$

(9)

$$\mathrm{U}=\frac{1}{\mathrm{n}}{\mathrm{h}}^{2/3}{\mathrm{J}}^{1/2}$$

(10)

$$\mathsf{\omega }=\frac{1}{24}\frac{{\mathsf{\gamma }}_{\mathrm{s}}-\mathsf{\gamma }}{\mathsf{\gamma }}\mathrm{g}\frac{{\mathrm{d}}^2}{\mathsf{\upsilon }}$$

(11)

Through the simultaneous solution method, Equation (8) can be transformed into:

$$\mathrm{S}*=\mathrm{k}{(\frac{24\mathsf{\gamma }\mathsf{\upsilon }\mathrm{J}}{{\mathrm{g}}^2({\mathsf{\gamma }}_{\mathrm{s}}-\mathsf{\gamma }){\mathrm{n}}^2}\cdot \frac{\mathrm{Q}}{{\mathrm{Bd}}^2{\mathrm{h}}^{2/3}})}^{\mathrm{m}}$$

(12)

where Q is the average flow (m³/s), d is the sediment particle size (mm), h is the average water depth (m), B is the average width (m), γ_s and γ refer to the volume weight of sediment and water respectively (N/m³), υ, J, and n refer to the viscosity coefficient of water movement, slope, and roughness of the river channel respectively, k is the coefficient, and m is a positive exponent. This equation can fully represent the hydrodynamic factor in the expression of cumulative erosion.

In fixed reaches, γ_s, γ, g, and υ can be taken as constants. Huang et al. [26] found that the channel slope in the studied reach was smaller and the variation range was small, and that the channel roughness was also stable after observations in recent decades. A large number of revetment projects in the studied reach can effectively ensure the stability of the channel width. Based on the above research results, this study can define the sediment carrying capacity factor as ${\mathrm{F}}_{\mathrm{i}}={\mathrm{Q}}_{\mathrm{i}}/{{\mathrm{d}}_{\mathrm{i}}}^2{{\mathrm{h}}_{\mathrm{i}}}^{2/3}$, which can approximately reflect the intensity of the sediment carrying capacity and has an important impact on the equilibrium value of the channel scouring and silting state. Moreover, the channel scouring and silting are also closely related to the sediment discharge, so the sediment discharge factor should be considered when establishing a calculation formula of cumulative erosion. However, the dimension of the sediment carrying capacity factor, F, is different from the sediment concentration, S, and this paper used a ratio of the two: the factor F_i^a/S_fi^b approximately reflects the size relationship between them. Therefore, the calculation formula of the equilibrium value representing the cumulative erosion of the Wuhan reach is as follows:

$${\mathrm{G}}_{\mathrm{ei}}=\mathrm{K}{{\mathrm{F}}_{\mathrm{i}}}^{\mathrm{a}}/{{\mathrm{S}}_{\mathrm{fi}}}^{\mathrm{b}}=\mathrm{K}{\left(\frac{{\mathrm{Q}}_{\mathrm{i}}}{{{\mathrm{d}}_{\mathrm{i}}}^2{{\mathrm{h}}_{\mathrm{i}}}^{2/3}}\right)}^{\mathrm{a}}/{{\mathrm{S}}_{\mathrm{fi}}}^{\mathrm{b}}$$

(13)

where G_ei is the equilibrium value of the annual cumulative erosion in the ith year. At present, the incoming flow and sediment regimes in the flood season are the main factors causing riverbed erosion and siltation in the Wuhan reach, so Q_i, S_fi, and d_i are the average flow in the flood season (m³/s), the average sediment concentration in the flood season (kg/m³), and the median diameter of suspended sediment in the flood season (mm) in the ith year respectively, which can be obtained from the measured hydrological data of the previous years of Hankou Station. Considering that the bankfull depth can fully reflect the longitudinal geometric characteristics of the channel, h_i is taken as the bankfull depth in the ith year, m, and K, a, and b are the empirical parameters obtained by regression.

Based on Equation (5) and given that ∆t = 1 year, the cumulative erosion of the Wuhan reach can be calculated with the weighted normalized DRM as follows:

$${\mathrm{G}}_{\mathrm{n}}=\frac{{1-\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}}}{{1-\mathrm{e}}^{-(\mathrm{n}+1)\mathsf{\beta }\Delta \mathrm{t}}}{\displaystyle \sum }_{\mathrm{i}=0}^{\mathrm{n}}{[\mathrm{e}}^{-(\mathrm{n}-\mathrm{i})\mathsf{\beta }\Delta \mathrm{t}}\mathrm{K}{\left(\frac{{\mathrm{Q}}_{\mathrm{i}}}{{{\mathrm{d}}_{\mathrm{i}}}^2{{\mathrm{h}}_{\mathrm{i}}}^{2/3}}\right)}^{\mathrm{a}}/{{\mathrm{S}}_{\mathrm{fi}}}^{\mathrm{b}}]$$

(14)

Li [35] found that the Coriolis force is an important parameter in the study of water bodies moving on the Earth’s surface. The Earth’s rotation produces Coriolis force, which affects rivers around the world, except the equator. When the water bodies in the northern (southern) hemisphere move along the meridian, they are affected by Coriolis force to the right (left), which results in the direction of the motion of water bodies moving toward the right (left). Therefore, this paper used Coriolis force as a hydrodynamic factor to study the variations of distance from the thalweg to the initial point in the studied reach.

The mathematical expression of Coriolis force is as follows:

$$\mathrm{F}=2\mathrm{mv}\mathsf{\omega }\sin \mathsf{\phi }$$

(15)

where F is the Coriolis force, m is the mass of the object, v is the speed of the object, ω is the rotational angular velocity for the Earth, and φ is the latitude at which the object is located.

The mass of water flow in motion is:

$$\mathrm{m}=\frac{\mathrm{V}\mathsf{\gamma }}{\mathrm{g}}$$

(16)

where V is the volume of water (m³) and γ is the volume weight of water (N/m³).

Given that the area of the river bank (one side) contacted by the flow in the reach is A (m²) and the average width of the reach is B (m), then:

$$\mathrm{m}=\frac{\mathrm{AB}\mathsf{\gamma }}{\mathrm{g}}$$

(17)

The Coriolis force per unit area is:

$$\mathrm{F}=2\mathrm{v}\mathsf{\omega }\mathsf{\gamma }\frac{\mathrm{B}}{\mathrm{g}}\sin \mathsf{\phi }$$

(18)

where v = Q/Bh, Q is the average flow (m³/s), and h is the average water depth (m). Equation (18) can be changed to:

$$\mathrm{F}= \frac{2\mathsf{\omega }\mathsf{\gamma }\sin \mathsf{\phi }}{\mathrm{g}}\cdot \frac{\mathrm{Q}}{\mathrm{h}}$$

(19)

where ω, γ, and g are fixed values. φ of the location of the Wuhan reach almost remains unchanged. The Coriolis force factor can be defined as F_i = Q_i/h_i, which can approximately reflect the intensity of the Coriolis force. In addition, the thalweg migration in the cross-section is accompanied by a certain degree of channel degradation. Therefore, the sediment discharge factor should also be considered when studying the thalweg migration. Similarly, because the dimension of the Coriolis force factor, F, is different from that of the sediment concentration, S, this paper used a ratio of the two, that is, the factor F_i^a/S_fi^b approximately reflects the size of the relationship between them. The calculation formula of the distance from the thalweg to the initial point is as follows:

$${\mathrm{V}}_{\mathrm{ei}}=\mathrm{K}{{\mathrm{F}}_{\mathrm{i}}}^{\mathrm{a}}/{{\mathrm{S}}_{\mathrm{fi}}}^{\mathrm{b}}=\mathrm{K}{\left(\frac{{\mathrm{Q}}_{\mathrm{i}}}{{\mathrm{h}}_{\mathrm{i}}}\right)}^{\mathrm{a}}/{{\mathrm{S}}_{\mathrm{fi}}}^{\mathrm{b}}$$

(20)

where V_ei is the equilibrium value of the distance from the thalweg to the initial point in the ith year, similar to the calculation of that of the cumulative erosion. Q_i and S_fi are the average flow (m³/s) and the average sediment concentration (kg/m³) in the flood season in the ith year respectively, h_i is the average bankfull depth in the ith year (m), and K, a, and b are the empirical parameters obtained by regression.

Based on Equation (5) and given that ∆t = 1 year, the distance from the thalweg to the initial point of the Wuhan reach can be calculated with the weighted normalized DRM as follows:

$${\mathrm{V}}_{\mathrm{n}}=\frac{{1-\mathrm{e}}^{-\mathsf{\beta }\Delta \mathrm{t}}}{{1-\mathrm{e}}^{-(\mathrm{n}+1)\mathsf{\beta }\Delta \mathrm{t}}}{\displaystyle \sum }_{\mathrm{i}=0}^{\mathrm{n}}{[\mathrm{e}}^{-(\mathrm{n}-\mathrm{i})\mathsf{\beta }\Delta \mathrm{t}}\mathrm{K}{\left(\frac{{\mathrm{Q}}_{\mathrm{i}}}{{\mathrm{h}}_{\mathrm{i}}}\right)}^{\mathrm{a}}/{{\mathrm{S}}_{\mathrm{fi}}}^{\mathrm{b}}]$$

(21)

4. Results

4.1. Changes in Flow and Sediment Regimes

In this paper, the flow and sediment data of Hankou Hydrological Station were collected according to four stages: before the impoundment of the TGR (1954–2002), during the cofferdam-supported power generation of the TGR (2003–2006, when the water storage level was 135 m), in the initial operation period of the TGR (2007–2008, when the water storage level was 156 m), and in the normal operation period of the TGR (2009–2019, when the water storage level was 175 m). The statistics are detailed in Figure 3 and

Table 3. Statistics of the flow and sediment discharges of previous years at Hankou Hydrological Station.

Period	Runoff (108 m3)			Sediment Discharge (108 kg)
	Total Annual Average	Annual Average in Flood Season	Proportion of Flood Season (%)	Total Annual Average	Annual Average in Flood Season	Proportion of Flood Season (%)
1954~2002	7147	5247	73.4	39,670	34,820	87.8
2003~2006	6734	4763	70.7	13,320	11,580	86.9
2007~2008	6589	4618	70.1	10,750	9077	84.4
2009~2019	6892	4692	68.1	9578	7664	80.0

. In general, although the annual runoff was almost unchanged, the TGP operation has greatly affected the annual sediment discharge and the annual law of runoff distribution. It can be seen in Table 3 that the average annual runoff at Hankou Hydrological Station in the four stages was 7147 × 10⁸, 6734 × 10⁸, 6589 × 10⁸ and 6892 × 10⁸ m³/year respectively, and the operation of the TGR had little impact on it. Nevertheless, the annual sediment discharges at Hankou Hydrological Station declined drastically due to TGR operation, from 39,670 × 10⁸ kg/year in the first stage to 9578 × 10⁸ kg/year in the fourth stage, with a decrease of 76%. Moreover, the distributions of runoff and sediment throughout the year were different. From the four different stages, the proportion of runoff and sediment discharge in the flood season in total declined in turn, indicating that the runoff and sediment discharge decreased during the flood season but increased in the non-flood season, under the influence of the TGR operation.

4.2. Morphological Adjustment

4.2.1. Variation in Vertical Channel Stability

According to the measured data of 57 fixed cross-sections in the Wuhan reach from 2003 to 2019, the bankfull depth and cumulative erosion after impoundment of the TGR could be calculated. The results are shown in Figure 4.

Figure 4 shows the increasing trend of the cumulative erosion and bankfull depth in the Wuhan reach after the impoundment of the TGR. The former reached 1.72 m, while the latter increased by 1.87 m in the past 17 years. Through the changes in incoming flow and sediment regimes of the previous years in Section 4.1, the main reason for significant degradation in the longitudinal profile is that there was little change in the upstream flow discharge, and a gradual reduction in sediment discharge. The variations of the thalweg longitudinal profiles along the riverbed of the Wuhan reach in 2003 (during the cofferdam-supported power generation of the TGR), 2008 (in the initial operation period of the TGR), 2014, and 2019 (in the normal operation period of the TGR) are shown in Figure 5. It can be seen that the riverbed thalweg elevation along the Wuhan reach significantly decreased with the increase in years after the impoundment of the TGR. Compared to the cofferdam-supported power generation and the initial operation period of the TGR, the riverbed thalweg elevation along the Wuhan reach decreased more greatly in the normal operation period of the TGR, which implied that the operation of the upstream TGR changed the characteristics of flow and sediment transport, and caused morphological adjustments of the downstream channel, thus resulting in longitudinal instability of the Wuhan reach. At present, the number of years of measurement data for fixed cross-sections after the impoundment of the TGR is small, and whether the longitudinal profiles of the Wuhan reach still maintain the erosion trend will require further study by increasing the years of measurement data in the future.

4.2.2. Variation in Lateral Channel Stability

Figure 6 and Figure 7 show the variations in the bankfull width and distance from the thalweg to the initial point in the Wuhan reach over the years after the impoundment of the TGR: (1) In 2003–2007, the average bankfull width of the Wuhan reach changed greatly from 1695 m in 2003 to 1676 m in 2004, and then to 1695 m in 2005, 1678 m in 2006, and up to 1711 m in 2007. After the completion of various revetment works in the Wuhan reach by 2007 (Figure 1), the bankfull width decreased gradually year by year and remained relatively stable, with a maximum annual variation of only 8 m. (2) In contrast, the distance from the thalweg to the initial point in the same period exhibited a significant increasing trend year by year, from 763.80 m in 2003 to 934.17 m in 2019, with an increase of 22.3%. That is to say, the overall thalweg of the studied reach has gradually shifted to the right over the years.

This study also selected two typical cross-sections in the Wuhan reach (HLZ06 and HLZ12-1) to observe the lateral morphological adjustments in 2003, 2008, 2014, and 2019, as shown in Figure 8.

Figure 8 shows that the left and right banks of HLZ12-1 remained stable and the channel width was basically unchanged due to the complete revetment works built on both sides (Figure 1), while a retreat of approximately 110 m occurred on the right bank of HLZ06 during 2003–2019 due to the lack of revetment projects. This indicates that the overall channel width in the Wuhan reach basically stabilized after the impoundment of the TGR due to the control effect of built revetment projects in most areas.

Furthermore, although the change trends of the left and middle parts of the two typical cross-sections are different, the right side has experienced an obvious erosion trend over the years. The thalweg point of the cross-section shifted substantially to the right from 2003 to 2019, and the change range was approximately 1 km. Li [35] found that the upper, middle, and lower reaches of the Yellow River experienced a natural phenomenon of shifting to the right bank, and explained that this was because of the Coriolis force of the water flow. Based on this mechanical law, the right bank of rivers in the northern hemisphere is seriously scoured. Coriolis force is also the cause of this phenomenon in the Wuhan reach.

To date, some studies have suggested that the direction of thalweg migration in meandering reaches is wavy, and there is no continuous trend in this migration toward the left or right [36], which is different from the characteristics of thalweg migration in the Wuhan reach, a straight braided channel. Similar to the study of cumulative erosion, the number of years of measurement data for the fixed cross-sections after the impoundment of the TGR is small, and whether the thalweg of the Wuhan reach continues to shift to the right will also require further study by increasing the years of the measurement data in the future.

4.3. Simulations of Channel Stability Indexes Based on the Improved DRM

According to the cumulative erosion and distance from the thalweg to the initial point of the Wuhan reach in 2003–2016 obtained in Section 4.2, the results calculated by the improved weighted normalized DRM (Equations (14) and (21)) are shown in Figure 9 and Figure 10, and the fitted parameters and corresponding correlation coefficients (R²) are shown in

Table 4. Fitting parameters, correlation coefficient, MNE, and p-values of the longitudinal and lateral stability indexes of the Wuhan reach using the DRM.

Indexes	K	a	b	β	R2	MNE	p
Cumulative erosion	0.003	0.08	2.17	0.14	0.93	19%	$<0.01$
Distance from thalweg to initial point	0.09	1.10	0.43	0.01	0.87	2.1%	$<0.01$

.

It can be seen from Figure 9 that good results were obtained when calculating the cumulative erosion and distance from the thalweg to the initial point of the Wuhan reach using two DRMs. By studying the incoming flow and sediment regimes of the previous five years (n = 4), it showed a high correlation between the calculated values and the measured ones. With the increase of n, the correlation coefficient R² remained essentially unchanged and fluctuated around the maximum value, with R² values of 0.93 and 0.87, respectively. The morphodynamics adjustments in the channel stability of our studied reach were closely related to the previous five-year flow and sediment regimes. Moreover, the accuracies of Equations (14) and (21) were verified by the measurements from 2017 to 2019 (Figure 10), implying that channel evolution may lag behind the changes in flow and sediment discharges, and this computational expression can be used to predict the variation trend of riverbed longitudinal and lateral stability with changes in hydrological conditions.

5. Discussion

According to the fitting results of the parameters in Table 4, it can be seen that: (1) The calibrated exponent a of hydrodynamic factors (sediment carrying capacity and Coriolis force) and the calibrated exponent b of the sediment discharge factor in the two models are all positive. The forms of Equations (14) and (21) in the two models reveal that the variations in channel morphology are positively correlated with the flow dynamics of the upstream water, and the larger the sediment concentration in the inflow, the smaller the riverbed degradation capacity, which is consistent with the basic laws in riverbed evolution. (2) By comparing the exponents of the different factors in the two models, the exponent of the sediment discharge factor in Equation (14) was much greater than that of the sediment carrying capacity factor, verifying that an important factor for significant degradation occurring in the longitudinal profile is the gradual decrease in the upstream sediment discharge. However, the exponent of the Coriolis force factor in Equation (21) was larger than that of the sediment discharge factor, which is similar to the research conclusion of Li [35]. The natural phenomenon that the river regime shifted to the right bank in the Wuhan reach is mainly because of the action of Coriolis force.

A large number of studies have shown that the previous flow and sediment regimes have an important impact on the riverbed evolution [9,15,16,37], but how to scientifically and reasonably explain such an impact in the physical sense is still a big challenge. The calculation methods proposed in this paper can quantify the influence weights of the previous flow and sediment regimes to a certain extent. Since the calculation accuracy of the cumulative erosion model is better than that of the distance from the thalweg to the initial point, the cumulative erosion was taken as an example for discussion below.

Table 5. Comparison of the influence weights of the incoming flow and sediment regimes of previous years on the cumulative erosion in the DRM (n = 4).

Note: The figures in the red circle represent the influence weights of the incoming flow and sediment regimes in 2006 for each year from 2006 to 2010.

lists comparisons of the influence weights of the incoming flow and sediment regimes of previous years on the cumulative erosion in the DRM. When n = 4, the calculation expression of the cumulative erosion using the DRM is as follows:

$${\mathrm{G}}_4=\frac{1-{\mathrm{e}}^{-\mathsf{\beta }}}{1-{\mathrm{e}}^{-5\mathsf{\beta }}}{\mathrm{G}}_{\mathrm{e}4}+\frac{(1-{\mathrm{e}}^{-\mathsf{\beta }}){\mathrm{e}}^{-\mathsf{\beta }}}{1-{\mathrm{e}}^{-5\mathsf{\beta }}}{\mathrm{G}}_{\mathrm{e}3}+\frac{(1-{\mathrm{e}}^{-\mathsf{\beta }}){\mathrm{e}}^{-2\mathsf{\beta }}}{1-{\mathrm{e}}^{-5\mathsf{\beta }}}{\mathrm{G}}_{\mathrm{e}2}+\frac{(1-{\mathrm{e}}^{-\mathsf{\beta }}){\mathrm{e}}^{-3\mathsf{\beta }}}{1-{\mathrm{e}}^{-5\mathsf{\beta }}}{\mathrm{G}}_{\mathrm{e}1}+\frac{(1-{\mathrm{e}}^{-\mathsf{\beta }}){\mathrm{e}}^{-5\mathsf{\beta }}}{1-{\mathrm{e}}^{-5\mathsf{\beta }}}{\mathrm{G}}_{\mathrm{e}0}$$

(22)

Then, the specific percentages of the G₄ occupied by the right five items can represent the influence weight of hydrological conditions in each year. For example, the influence weight of the current year’s hydrological conditions can be expressed as $(1-{\mathrm{e}}^{-\mathsf{\beta }}){\mathrm{G}}_{\mathrm{e}4}/(1-{\mathrm{e}}^{-5\mathsf{\beta }}){\mathrm{G}}_4$ and that of previous conditions can be expressed in the same way.

To sum up, the average influence weights on the riverbed stability in the current year and the previous four years were 32%, 24%, 19%, 14%, and 11%, respectively (Table 5). Table 5 not only indicates that the changes in the riverbed stability were mostly caused by the previous hydrological conditions (the total weights accounted for 68%), but also proves that the more recent hydrological conditions exert a greater impact, in accordance with the basic understanding of channel evolution.

In addition, Table 5 shows that the peak value of the influence weight of the flow and sediment regimes in the current year appeared in 2006 (70%), significantly higher than that in other years. In the following years, 2007–2010, the influence weight in 2006 was still the largest. This is because the year 2006 was the most typical dry year after the impoundment of the TGR; thus, the sharp decline in the incoming flow and sediment naturally has a considerable influence on the riverbed stability. Therefore, the model calculation method established in this paper can sensitively reflect the response process of riverbed stability variations to sudden changes in inflow and sediment regimes. However, it should be noted that for the different characteristic parameters of riverbed stability variations, the previous incoming flow and sediment regimes may have different impacts on them. The reason is that when different characteristic parameters of riverbed stability variations are adjusted under the action of flow and sediment regimes, due to the different river boundary conditions, the sensitivity to flow and sediment variations and human activities, as well as their self-adjustment rates, are different. The faster the adjustment rate, the smaller influence the flow and sediment regimes will have on it, and vice versa.

6. Conclusions

The variation characteristics of the longitudinal and lateral stability of the Wuhan reach in the Yangtze River were comprehensively analyzed in this study, based on the flow and sediment data of Hankou Hydrological Station and the measured topographic data of 57 fixed cross-sections in the Wuhan reach in 2003~2019. The main conclusions are as follows:

(1) The channel has been obviously eroded in the downstream reach, resulting from a 76% decrease in sediment discharge after the impoundment of the Three Gorges Reservoir (TGR). The cumulative erosion and average bankfull depth in the Wuhan reach exhibited a significant growth trend in the span of 17 years, with the former reaching 1.72 m, and the latter increasing by 1.87 m. After the year 2007, with the completion of various revetment works, the average bankfull width in the Wuhan reach remained stable. The distance from the thalweg to the initial point increased in response to the Coriolis force, with an increase of 22.3%.

(2) A method for calculating the equilibrium values of cumulative erosion and distance from the thalweg to the initial point was established after an analysis of the effects of flow, sediment concentration, sediment gradation, and water depth, combined with the improved weighted normalized DRM calculation formula. The parameters in this method were calibrated with the measured values from 2003 to 2016, and good results were obtained for them. The morphodynamics adjustments in the channel stability of our studied reach were closely related to the previous five-year flow and sediment regimes. Moreover, the accuracy of the equations were verified by measurements from 2017 to 2019, implying that channel evolution may lag behind the changes in flow and sediment discharges, and this computational expression can be used to predict the variation trend of riverbed longitudinal and lateral stability with changes in hydrological conditions.

(3) By analyzing the comparisons about the influence weights of the incoming flow and sediment regimes on riverbed stability in the DRM over previous years, it was found that the previous hydrological conditions have significant effects on the changes, and the recent hydrological conditions exert the greatest impacts. This phenomenon conforms to the basic understanding of channel evolution. The calculation methods established in this paper can reflect the extent of the influence of the interannual variations of flow and sediment regimes on the change process of riverbed stability to a certain extent, and provide technical support for channel improvement.