Rivers are open systems with clear spatial boundaries which exchange material and energy with the outside world. As a medium for water and sediment transport, they transport water and sediment from the earth's surface to the ocean, playing an important role in shaping the earth's surface morphology (Knighton, 1998). With the downstream movement of water and sediment, the riverbed undergoes continuous erosion and sedimentation, manifested at the micro level as the gradual approach of sediment concentration to sediment transport capacity and at the macro level as the development of river channel morphology towards a certain hydraulic geometry (Qian et al., 1987; Wu and Zheng, 2015). As a result, rivers can approach a relatively equilibrium state over a long period. Generally, in the process of adjustment from the non-equilibrium state to the equilibrium state, the initial adjustment rate is relatively large. As the system gradually approaches the equilibrium state, the adjustment rate will decrease gradually. Furthermore, when the system approaches equilibrium infinitely, the adjustment rate approaches zero. This characteristic in the fluvial system can be represented by the rate law or exponential decay equation. For example, Graf (1977) introduced the rate law into the study of fluvial processes and used an exponential decay formula to describe the temporal variation process of gully length, achieving good simulation results. The study by Surian et al. (2003) also indicates that the process of riverbed incision under the influence of human activities can be fitted by using an exponential decay function. However, the simple exponential decay formula cannot reflect the dynamic adjustment process of characteristic variables in response to external environmental controls.

Changes in riverbed erosion and deposition are the results of the cumulative movement of sediment particles. Macroscale riverbed deformation cannot be accomplished at one stroke, in fact, this is a slow and continuous process. Therefore, the riverbed evolution always delays behind the changes in water and sediment conditions (Wu et al., 2006; Wu et al., 2007; Zhang et al., 2010). The research of Xia et al. (2016) indicates that in the Jingjiang reach of the Middle Yangtze River, there is a good correlation between the variations in bankfull area (as well as water depth) and the water and sediment conditions of the previous five years in the flood season. Wu et al. (2007, 2012) derived a Delayed Response Model (DRM) based on the spontaneous adjustment principle and rate law of fluvial systems, including the universal integral mode, single step analytical mode, and multistep iterative mode, providing an effective simulation method for the non-equilibrium riverbed evolution process (Wu et al., 2008; Zheng et al., 2015; Lyu et al., 2019; Wang et al., 2020b; An and Fu, 2022). However, the DRM can only simulate the temporal adjustment process of characteristic variables and does not involve spatial adjustment.

Non-equilibrium sediment transport is the fundamental reason for riverbed evolution. Simon et al. (2006) analysed the evolution characteristics of rivers with excess flow energy and indicated that during the process of river adjustment, flow energy will gradually dissipate. Driving factors such as stream power and energy slope will decay nonlinearly over time, while hindering factors such as roughness, critical shear stress, and riverbed sediment particle size will increase nonlinearly over time. Jing et al. (2020) provided a mathematical description and theoretical model of the spatial evolution of internal feedback induced by external disturbances in alluvial rivers based on statistical mechanics. They simulated the adjustment process of the river channel after the construction of the dam, and the results indicated that the erosion volume decreased nonlinearly in downstream direction until it disappeared. The non-equilibrium suspended sediment transport equation, as a commonly used formula in fluvial processes research, can easily and effectively describe the spatial variation of sediment concentration (Xie, 1990) and analyse the saturation level of sediment-laden flow (Liu et al., 2022). According to the non-equilibrium sediment transport equation, if the sediment concentration is less than the sediment transport capacity, the water flow will scour the riverbed to supplement sediment. In this process, the sediment concentration increases along the stream-wise direction, and the corresponding sediment concentration's spatial variation rate is positive. Whereas, if the sediment concentration is greater than the sediment transport capacity, sediment particles will deposit on the bed surface. In this situation, the sediment concentration decreases in the stream-wise direction, and therefore, the sediment concentration's spatial variation rate is negative. Non-equilibrium sediment transport equation is often used to describe the variation of sediment concentration or sediment transport rate along the river. It also provides a theoretical basis for analysing the spatial variation of riverbed erosion and deposition. Compared with the empirical formula of “more incoming more desilting” (Wang et al., 2016), it has a more rigorous theoretical basis.

The erosion and deposition evolution of the Lower Yellow River (LYR) caused by variations in water and sediment conditions (Long, 2002; An et al., 2020) and the extension of the estuary (Qian and Zhou, 1965) is a typical case of the long-term and non-equilibrium riverbed evolution process under the joint influence of upstream water and sediment inflow and downstream erosion base level. Based on the simulation results of the sediment mathematical model, Li et al. (1989) indicated that the adjustment of the longitudinal profile of the LYR is the result of the joint and interactive effects of downward deposition and headward deposition. The study by Lu et al. (2003) indicated that the longitudinal profile of the LYR is related to factors such as water flow energy, aggradation environment, and estuary extension. Meanwhile, the water flow energy and aggradation environment reflect the influence of water and sediment conditions, while the estuary extension reflects the effect of erosion base level. Since the operation of the Xiaolangdi Reservoir in 1999, the sediment concentration of the flow entering into the LYR has decreased significantly, leading to a long-term, non-equilibrium and continuous erosion process in the LYR, and the spatiotemporal adjustment process of riverbed erosion and deposition is very complex. The existing research has analysed the evolution law of the LYR since the operation of the Xiaolangdi Reservoir, such as longitudinal profile adjustment (Xia et al., 2008), cross-sectional profile adjustment (Wang et al., 2020a), sediment transport characteristics (Zhang et al., 2018), and the relationship between reservoir operation and channel erosion and deposition (Zheng, 2013). However, research on the spatiotemporal variation law and internal mechanism of riverbed evolution in the LYR is still very limited.

This research focuses on the spatiotemporal erosion and deposition process in the LYR since the operation of Xiaolangdi Reservoir. Based on the DRM for riverbed evolution, a calculation method for the simulation of the temporal process of the accumulated erosion and deposition volume in the LYR is established. Meanwhile, based on the non-equilibrium suspended sediment transport equation, a calculation method for simulating the spatial variation of the accumulated erosion and deposition volume per unit river length of the main channel in the LYR is established. Based on the established method, the spatiotemporal variation process of erosion and deposition in the LYR since the operation of the Xiaolangdi Reservoir was simulated, and the spatiotemporal characteristics of the adjustment of erosion and deposition in the river channel were revealed. The research results can provide a scientific reference for the quantitative study of the spatiotemporal process of strong non-equilibrium riverbed erosion and deposition at the downstream of dams.

The LYR (Figure 1) originates from Taohuayu on the southeast slope of the Loess Plateau and joins the Bohai Sea near Lijin. The total length of the LYR is approximately 786 km, and the river's cross-sectional profile is wide upstream and narrow downstream, with the slope of riverbed steep upstream and gentle downstream. There are a total of 7 hydrological stations in the LYR: Huayuankou (HYK), Jiahetan (JHT), Gaocun (GC), Sunkou (SK), Aishan (AS), Luokou (LK), and Lijin (LJ). According to the differences in river patterns, the LYR can be divided into three reaches. The reach upstream of Gaocun is the braided reach, which has a length of approximately 284 km, the mainstream of this reach changes significantly, and the riverbanks are prone to be eroded. The reach between Gaocun and Taochengpu is the transitional reach, with a length of approximately 184 km, and the mainstream of this reach is relatively stable. The reach between Taochengpu and Lijin is the meandering reach, which has a length of approximately 272 km, with the most stable mainstream among the three river reaches (Xia et al., 2014).

In this research, the input station is the Huayuankou station, and the measured water discharge, sediment concentration data, and coefficient of incoming sediment defined as the ratio of sediment concentration to corresponding water discharge from 2000 to 2020 were collected and organized, respectively. The output station is the Lijin station, and the water level at the discharge of 3000 m³/s at the Lijin station is taken as the erosion base level of the LYR. The erosion and deposition data are the annual measured accumulated deposition from the Huayuankou station to the Lijin station during the period of 2000-2020 (based on the deposition value in 1999), obtained by analysing the measured cross-sectional data of the LYR. The sediment particle size data in this research are based on the measured bed sediment particle size data from 7 hydrological stations in the LYR, represented by the median particle size D₅₀ of the bed sediment.

The temporal variations of the coefficient of incoming sediment at the input station and the water level at the output station are shown in Figure 2. Figure 2a shows that there is a significant decrease in the annual average coefficient of incoming sediment at Huayuankou before and after the year 2000. After 2000, the coefficient of incoming sediment remained at a low level for a long time, with an annual average value fluctuating between 0.0005 and 0.0117 kg∙s/m⁶. The multiyear averaged value is 0.0045 kg·s/m⁶, which accounts for only 15.4% of the value in previous period (1986-1999).

With the operation of the Xiaolangdi Reservoir in 1999, the LYR gradually transformed from continuous deposition to continuous erosion, and the significant change in water and sediment conditions was one of the main reasons for the transition. From Figure 2b, it can be revealed that during the erosion period (2000-2020), the water level at the discharge of 3000 m³/s at Lijin station shows a decreasing trend over time, which is prone to form the headward erosion from downstream to upstream.

Based on the measured erosion and deposition data of the LYR, the calculation method for accumulated deposition represented by Eqs. (1) and (7) is adopted to calibrate the parameters in the formula. The parameters taken to evaluate the accuracy of the calibration results are the coefficient of determination R² and the relative error MNE. The corresponding calculation formulas are listed as follows:

where y_c,i and y_m,i represent the calculated and measured values of the characteristic variable y in the ith time period, respectively, while $\overline{y_{c}}$ and $\overline{y_{m}}$ represent the average calculated and measured values of the characteristic variable y over N time periods, respectively. R² quantifies the degree of correlation between the calculated values and the measured values of the characteristic variable. The closer the value of R² is to 1, the higher the correlation between the calculated values and the measured values; MNE quantifies the average relative error between the calculated values and the measured values of the characteristic variable. The closer the value of MNE is to 0, the higher the calculation accuracy.

After the operation of the Xiaolangdi Reservoir in 1999, the state of the riverbed of the LYR began to transform from deposition to erosion. The value of accumulated deposition and erosion base level in this research are based on the value in 1999 (set as 0 for both the accumulated deposition and the water level at the discharge of 3000 m³/s at Lijin station in 1999), and the calculation period is from 2000 to 2020. Based on the measured water and sediment data and accumulated deposition data from 2000 to 2020, taking ΔT=1 a as the time interval, the following formula for calculating the annual accumulated deposition in the LYR can be obtained through a nonlinear regression method:

where the decay rate β=0.23/a.

the formula of the equilibrium value of accumulated deposition is:

where $\bar{Q}_{n}$ and $\bar{S}_{n}$ are the annual average water discharge and sediment concentration at Huayuankou station in nth time period, and ΔZ_n is the accumulated variation value of the water level at Lijin station at the discharge of 3000 m³/s over n time periods.

The comparison between the calculated and measured values of accumulated deposition in the LYR from 2000 to 2020 is shown in Figure 5. The coefficient of determination between the calculated and measured values is equal to 0.99, and the average relative error is 16.5%, indicating high calculation accuracy of the DRM.

Due to changes in incoming water and sediment conditions, as well as the downstream erosion base level, the equilibrium value of the accumulated deposition in the LYR has been constantly changing. Meanwhile, taking V_e as the target, the deposition volume V is adjusting constantly. In general, the relationship between the equilibrium value and the calculated value of accumulated deposition can be expressed as follows:

Figure 6a shows the relationship between the equilibrium value and the calculated value of accumulated deposition in the LYR, reflecting the delayed response characteristics between the accumulated deposition and the corresponding equilibrium or target value during each time period. Figure 6a reveals that the current equilibrium value of accumulated deposition from 2000 to 2020 is smaller than the calculated value at the end of the previous period. Therefore, during this period, the accumulated deposition adjusted in a decreasing direction, and the LYR experienced continuous erosion. After 2000, the erosion volume per unit runoff gradually decreased (Figure 7), indicating that with the continuous erosion, the riverbed gradually coarsened (Figure 12), and the erosion rate of the riverbed slowed down. Moreover, the difference between the equilibrium value of the current period and the calculated value at the end of previous period is decreasing in recent years, indicating that riverbed erosion is gradually approaching equilibrium (Figure 6a). Based on the Eq. (9), the temporal variations of the equilibrium slope of the generalized deposition body can be obtained (Figure 6b). From 2000 to 2020, with the continuous erosion in the LYR, the temporal variations of the equilibrium slope of the riverbed showed a decreasing trend, and the rate of decrease gradually slowed down. The decrease in the equilibrium slope is beneficial for slowing down the erosion rate. Meanwhile, with the coarsening of the riverbed (Miao et al., 2016), the adjustment of riverbed erosion is gradually approaching equilibrium, which reflects the equilibrium tendency of riverbed evolution. In addition, according to the derivation in section 3, K₂=BL. Based on the value of K₂ obtained from the calibration, the generalized river width is equal to 791 m, which is within the variation range of the bankfull width in the LYR (377-1226 m) (Wang et al., 2023).

Based on the single step mode of the DRM and referring to the formula of adjustment time of the characteristic variable in previous research (Julien, 2018), the adjustment time of the accumulated deposition volume can be expressed as:

where V is the adjustment value of the accumulated deposition, V₀ is the initial value of the accumulated deposition, and V_e is the equilibrium value of the accumulated deposition.

The adjustment amount required for the accumulated deposition V_n to reach the equilibrium value V_e,n from the initial value V_n-1 is defined as the target adjustment amount ΔV_n, which means ΔV_n = V_e,n - V_n-1. The adjustment amount of erosion and deposition during the period is defined as the adjustment completion amount $\Delta V_{n}^{\prime}, \quad \Delta V_{n}^{\prime}=V_{n}-V_{n-1}$. Based on the DRM, the equilibrium value V_e,n can only be infinitely approached, which means V_e,n cannot be reached. To study the adjustment completion degree of erosion and deposition, the ratio of the adjustment completion amount $\Delta V_{n}^{\prime}$ to the target adjustment amount ΔV_n is defined as completion degree m (the range of m is 0-1), and the time required to accomplish the completion degree is defined as the adjustment time T_m, where the erosion and deposition amount V_n = V_n-1 + m(V_e,n - V_n-1). Substituting the expression of V_n into Eq. (26), the calculation formula for adjustment time T_m can be obtained:

Based on Eq. (27) and substituting in the value of decay rate β, the adjustment time T₅₀, T₆₀, T₇₀, T₈₀, T₉₀, and T₉₅ required for achieving adjustment ratio of 0.50 (half of the target adjustment amount), 0.60, 0.70, 0.80, 0.90, and 0.95 can be calculated, respectively. With the operation of the Xiaolangdi Reservoir in 1999, the sediment load entering the LYR has sharply decreased, leading to long-term and continuous erosion. Based on the calculation results in section 4, the decay rate β=0.23/a. By substituting β into Eq. (27), the time required for the adjustment of deposition to reach different completion degree can be calculated (Table 1). For example, in the background of decreasing sediment inflow and river channel erosion, for the initial riverbed of any year, under the conditions of constant upstream water and sediment conditions and downstream erosion base level, the time required to accomplish half of the accumulated deposition adjustment is approximately 3.0 years, and the time required to reach a quasi-equilibrium state (95% of the accumulated deposition adjustment to be completed) is approximately 13.0 years.

The spatial distribution characteristics of riverbed erosion and deposition in the LYR are the result of the coupling effect of downward and headward erosion and deposition. The spatial distribution of accumulated main channel erosion and deposition per unit river length in the LYR from 2003 to 2015 can be simulated based on Eq. (19). The calibration results of the key parameters such as the erosion and deposition at the input station V_a, the sediment transport capacity variation term V_b, and the spatial adjustment index φ are shown in Table 2, and the comparison between calculated and the measured values of accumulated main channel erosion and deposition per unit river length is shown in Figure 8 (the figure only lists the results of odd years).

From Table 2 and Figure 8, it can be concluded that Eq. (19) can effectively simulate the spatial distribution of accumulated main channel erosion and deposition per unit river length in the LYR. The coefficient of determination R² between the fitted values and the measured values is higher than 0.98, and the relative error is lower than 18%. Table 2 reveals that the absolute value of V_a during the erosion period increased over time. Considering that the coefficient of incoming sediment shows a decreasing trend since 2000 and has remained at an extremely low level for a long time, the main reason for the temporal increase in the absolute value of V_a may be the water flow at the input station has been in an unsaturated state and needs to scour the riverbed to supplement sediment. It can be foreseen that if the current water and sediment conditions at the input station are maintained in the future, the absolute value of V_a will continue to increase. However, due to the coarsening of riverbed sediment, which means the increase in resistance to water flow, the increase rate of V_a may slow down.

Note that the erosion amount in the LYR tends to approach the spatial equilibrium value V_b near the lower part of the river, rather than approaching 0, and the absolute value of the spatial equilibrium value V_b shows an increasing trend over time. It can be concluded that the sediment transport capacity of the water flow in the LYR has a spatial enhancement trend to some extent, which may be due to factors such as narrowing of the cross-section, fining of bed sediment, and decrease in erosion base level. The combined effect of these factors induces that the water flow still has erosion ability near the output station of the LYR, enabling the continuous development of riverbed erosion along the river.

Figure 9 shows the spatial variations of sediment concentration and coefficient of incoming sediment along the LYR. From Figure 9, it can be concluded that the sediment concentration has an approximately linear increasing trend in space, whereas the coefficient of incoming sediment shows an exponential increasing trend along the river. Based on Eq. (13), combined with the calibration results in Table 2 and Figure 8, it can be concluded that the sediment transport capacity of the water flow in the LYR is generally higher than the sediment concentration after the operation of the Xiaolangdi Reservoir. Therefore, a widespread and continuous erosion exists from upstream to downstream, and the water flow supplements sediment through erosion along the way, causing the sediment concentration to show an increasing pattern from upstream to downstream. This is a typical case of the spatial adjustment of sediment-laden flow from non-equilibrium to equilibrium.

The Shields number reflects the relative magnitude of the shear stress exerted on sediment particles by water flow and the resistance force exerted on sediment particles. The larger the Shields number, the more easily the bed sediment is washed away by water flow. The calculation formula is:

where γ is the unit weight of clear water; γ_s is the unit weight of sediment particles; J is the slope of the river, which can be approximately replaced by the riverbed slope; and D is the particle size of bed sediment; h is the average water depth of main channel before and after flood season.

The changes in the unit weight of clear water and sediment particles are not significant, therefore these two parameters can be approximately treated as constants. Taking the median particle size of the bed sediment D₅₀ as the representative particle size, then Eq. (28) can be simplified as the following flow intensity parameter (Wu and Zheng, 2015):

where h is the average water depth in the main channel of the LYR before and after the flood season, obtained based on measured cross-sectional data of each reach in the LYR; J is the riverbed slope of the LYR, obtained by referring to relevant research (Wang et al., 2023); $\bar{D}_{50}$ is the average median bed sediment particle size in each river reach, for example, the average median particle size of riverbed sediment of the Hua-Jia reach is the average value of median particle size of riverbed sediment at Huayuankou station and Jiahetan station.

The flow intensity parameter σ can reflect the erosion capacity of water flow, and the calculation results are shown in Figure 10. According to Figure 10, it can be concluded that, in general, during the continuous erosion period after 2000, the flow's erosion capacity towards the riverbed increased from upstream to downstream in the LYR, which may be one of the reasons why the water flow of the LYR still has a certain erosion capacity when approaching the output station.

The temporal variations of spatial adjustment index of erosion and deposition φ is shown in Figure 11a. In general, the value of φ shows a decreasing trend over time. This decreasing trend indicates that the adjustment of the accumulated erosion and deposition per unit river length from the value at the input station V_a to the spatial variation term of sediment transport capacity V_b gradually slows down. To analyse the uniformity of the spatial distribution of erosion in the LYR, this research proposed the parameter L_0.5. The midpoint of erosion L_0.5 is defined as the position where the accumulated erosion originating from Huayuankou (input station) reaches half of the total erosion in the LYR. The smaller the L_0.5, the more concentrated the riverbed erosion distributed in the upper part of the river, and the more nonuniform the spatial distribution of riverbed erosion, whereas the larger the L_0.5, the more uniform the spatial distribution of riverbed erosion. Based on the values of V_a, V_b, and φ listed in Table 2, it is not difficult to obtain the annual calculated value of L_0.5 through integration, and the calculation result of L_0.5 is shown in Figure 11b.

The year 2000 when erosion began in the LYR is defined as the initial year. As shown in Figure 11b, L_0.5 shows an overall increasing trend over time, indicating that with the continuous development of riverbed erosion, the spatial distribution of accumulated erosion in the LYR will gradually become uniform. When there is less sediment and more water entering the LYR, the riverbed undergoes erosion from upstream to downstream, and the sediment concentration in the water gradually recovers from upstream to downstream. Compared with the lower part of the LYR, the upper part of the LYR has a greater degree of sediment unsaturation, therefore the erosion capacity of water flow on the riverbed is more significant, forming a decreasing spatial pattern of accumulated erosion in recent years.

Corresponding to the spatial pattern of accumulated erosion, the riverbed composition exhibits a severe coarsening phenomenon (Figure 12a), and the median particle sizes of bed sediment also shows a decreasing trend from upstream to downstream (Figure 12b). In addition, the upper part of the LYR has a relatively larger amount of erosion, therefore, its riverbed slope has a decreasing trend with the continuous development of erosion. Under the joint effect of the above factors, when erosion develops to a certain degree, the boundary conditions of the lower part of the LYR are more prone to be eroded than those of the upper part of the river. Therefore, in recent years, when the accumulated erosion per unit river length maintained a decreasing pattern from upstream to downstream in the LYR, the difference in erosion amounts between the upper and lower parts of the LYR decreased gradually. This process reflects the equilibrium tendency of riverbed evolution.

(1) Based on the DRM and the generalized equilibrium geometry of the deposition body, the method for calculating the accumulated deposition in the LYR has been established. The coefficient of determination R² between the calculated and measured values from 2000 to 2020 reached 0.99, indicating good calculation results. With the continuous erosion of the riverbed by water flow, the erosion rate of the riverbed in the LYR shows a decreasing trend.

(2) Based on the non-equilibrium sediment transport equation, the formula for calculating the spatial distribution of accumulated main channel erosion and deposition per unit river length has been derived. The formula can effectively simulate the spatial distribution characteristics of erosion and deposition amount in the LYR. The coefficient of determination R² between the fitted values and the measured values from 2003 to 2015 is between 0.98 and 0.99, with an average relative error of 6.2% for the research period.

(3) Since the operation of Xiaolangdi Reservoir, the equilibrium riverbed slope of the LYR has been decreasing over time, and the rate of decrease has gradually slowed down. For the initial riverbed of any year, if the upstream water and sediment conditions and the downstream erosion base level remain constant, the time required to complete half of the accumulated deposition adjustment is approximately 3.0 years, and the time required to complete 95% of the adjustment is approximately 13.0 years.