Analysis of the effect of unsteady reservoir flow and sand excavation on the channel conditions of the Yangtze River (Yibin-Chongqing)

This paper systematically describes the current research methods for analyzing unsteady flows in an open channel. In particular, it focuses on elaborating a discrete method to solve the one-dimensional Saint–Venant equations for unsteady flow in an open channel, compares the advantages and disadvantages of different discrete methods, and proposes a new implicit iterative mathematical model applicable to the calculation of the multiwavelength unsteady flow in the Yangtze River from Yibin to Chongqing. This paper also simulates three typical unsteady flow processes and analyzes the effect of sand excavation on the water level in the upper reaches of the Yangtze River from Yibin to Chongqing considering flood, normal-water, and low-water conditions. The spreading characteristics of the daily adjusted unsteady flow were obtained in the Shuifu-Chongqing river reach. Under the conditions of the daily minimum unsteady discharge flow rate, the designed water level of the reach from Yibin to Chongqing increases by 0.19 m compared to the natural water level. With the increase in concentrated sand excavation in the upper reaches of the Yangtze River, a large number of shingle beaches have disappeared, and the riverbed has been drastically lowered. Affected by the excessive and unordered sand excavation in the river channel, the water level of the upper reaches of the Yangtze River from Chongqing to Yibin has decreased by 0.5–2.0 m. This work improves the Preissmann implicit difference method to effectively simulate the propagation characteristics of unsteady flow in long reaches. The daily adjustment of unsteady flow in the upper reaches increases the river flow during the low-water period so that the designed water level under the action of this daily unsteady flow adjustment increases by approximately 0.19 m compared to that under normal conditions. A large number of point bars have been damaged in the wake of sand excavation in long reaches,resulting in a water level decline of approximately 0.5–2.0 m. This work improves the Preissmann implicit difference method to effectively simulate the propagation characteristics of unsteady flow in long reaches. The daily adjustment of unsteady flow in the upper reaches increases the river flow during the low-water period so that the designed water level under the action of this daily unsteady flow adjustment increases by approximately 0.19 m compared to that under normal conditions. A large number of point bars have been damaged in the wake of sand excavation in long reaches,resulting in a water level decline of approximately 0.5–2.0 m.


Introduction
In studies on unsteady flow, there were few notable breakthroughs until the Saint-Venant [1] equations were established. Saint-Venant conducted research to derive the Saint-Venant partial differential equation for unsteady flow. Massau [2] proposed the use of the method of characteristics to solve the equation for unsteady flow in an open channel. Richardson [3] used partial differential equations to analyze unsteady flow, which is considered the first application of numerical simulation in hydraulics. Pereira [4] studied the fourthorder Runge-Kutta scheme, which has good stability. Pappenberger [5] conducted a stability limit test on the Preissmann four-point eccentric implicit scheme based on the Hydraulic Engineering Center's River Analysis System (HEC-RAS). Their results show that at the limit value, the Preissmann four-point eccentric implicit scheme is still very stable. Freitag [6] employed the improved Thomas chasing method to solve the Preissmann fourpoint eccentric implicit scheme and achieved good results. With a limited number of research methods, research on this topic has been slow over the years, and few major breakthroughs have been made.
Sneider et al. [7] established an unsteady flow level prediction model for dam discharge to the lower reach. Theiling et al. [8] studied the effect of dam adjustment in the upper reach of the Mississippi River on the flow conditions of the waterway. Sokolov [9] studied the celerity of long waves in channels with complex cross sections. Toro [10] studied how the unsteady flow discharge of the Lagdo Hydropower Station causes a decrease in incoming flow in the lower reach during the low-water period. French [11] combined the Manning-Chezy formula (they assumed that it is applicable) with the Saint-Venant equations, deducing that the unsteady flow rate-water level relationship forms a loop-rating curve. Kundzewicz et al. [12] demonstrated wave attenuation under unsteady flow. Tu and Graf [13] believed that the increasing time of unsteady flow is approximately one-third of the cycle when the unsteady flow propagates for one wavelength. Many studies have shown that wave attenuation and deformation occur during the spreading of unsteady flow Lamberti and Pilati [14]. The wave amplitude decreases as the spreading distance increases. Through mathematical analysis, Kundzewicz and Dooge [15] provided an expression for the attenuated amplitude of unsteady flow waves.
The Yangtze River (Yibin-Chongqing) is 356 km in length. Research results have shown that human activities impact the variation in water and sediment conditions in the upper reaches of the Yangtze River, (Marwan et al. [16], Matos et al. [17], Cheng et al. [18]). Wang et al. [19] showed that the average annual sediment in the upper reaches of the Yangtze River has decreased compared with that before the completion of the reservoir group. In recent years, a large amount of sand and gravel has been extracted in the Yibin to Chognqing reach of the Yangtze River. In addition, the negative effects of sand and gravel mining on river channels have been studied (Xiao et al. [20], Xiao et al. [21]). Lu et al. [22] analyzed the impact of sand mining on the sedimentation of the Three Gorges Reservoir. Wang et al. [23] investigated the relation between sand mining and runoff-sediment discharge and concluded that the sediment transport intensity was lower in the sand-mined river reach. The previous studies all focused on local river reaches, whereas results concerning long reaches of the upper Yangtze River have rarely been reported.
The Xiangjiaba Hydropower Station is the last-stage power station in the Jinsha River hydroelectric cascade development. Ma [24] and Liu [25] showed that the unsteady flow caused by daily adjustments in the dam release changes the natural flow conditions of the river downstream. Currently, studies on the unsteady flow of an open channel in long reaches are still based on the mathematical model of water flow. For the study of unsteady flow downstream of Xiangjiaba, which is limited by the test conditions, research within 100 km downstream of Xiangjiaba was selected, but this selection cannot fully indicate the spreading features of the unsteady flow. Zhang [26] showed that sand mining has damaged the channel conditions to different degrees. Until now, few studies have been conducted on the unsteady flow and sand excavation in the Yangtze River (Yibin-Chongqing) given new boundary conditions.
The current study has eight sections. Section 2 presents an overview of the studied river reach. Section 3 presents a review of the numerical model and solution method. Sections 4 and 5 present the propagation characteristics of unsteady flow and the impact of unsteady flow on the water level of the river section. Section 6 presents the gravel excavation pit characteristics in the upstream reach of the Yangtze River based on field measurements and the impact of sand mining on the water level of the channel during different water level periods. Finally, Sect. 7 describes and discusses the results. In this paper, the one-dimensional unsteady flow mathematical model is improved. Studies on the daily changes in the unsteady flow in the reaches downstream of the dam and the characteristics of the water level drop in these reaches after sand excavation can provide technical support for the implementation of channel improvement projects in the upper reaches of the Yangtze River.

Overview of the studied river reach
The Jinsha River is 356 km in length from Shuifu (upstream channel 1076 km) to Jiangjin of Chongqing (720 km). At present, the river reach from Shuifu to Yibin is approximately 31 km long. As a Class IV channel, it has dimensions of 2 m × 40 m × 300 m (waterway depth × waterway width × bending radius). Fleets of 220 kW + 2 × 300 t can navigate this reach year round. For the Yibin-Jiangjin channel, the technical level reaches Class III, and the maintenance dimensions are 2.9 m × 50 m × 560 m. During the low-water period, this reach can accommodate 1,000-ton ships and their fleets. Figure 1 shows the area of research in the upstream reach of the Yangtze River. The mileage of the main site is shown in Fig. 2. 3 Review of the numerical model and solution method

Comparison of the Saint-Venant equation solution method
Currently, due to calculation accuracy, the instantaneous flow regime method and the microamplitude wave theory typically used to solve the Saint-Venant equations are no longer used. According to the complicated calculation requirements, the method of characteristics needs to classify and solve a variety of flow regimes and requires very complete data. Thus, this method is seldom used. The discrete finite difference method is intuitive and simple and thus has been widely used. Anderson [27] presented a systematic summary of the finite difference method and divided this method into the explicit difference method and the implicit difference method. The unsteady flows of the Yibin-Chongqing river reach exhibit great changes and a long spreading period. We  is not considered in the partial derivative of the distance calculation.
The influence of f j+1 i+1 can cause calculation errors that continue to increase with increasing distance and time step. Clearly, this method cannot be applied to the 500-kmlong river reach from Yibin to Chongqing.

Scheme stability problem
All explicit schemes exhibit stability problems. The solutions remain stable only when the time interval Δx and distance interval Δx meet certain conditions. For example, the stable scheme conditions of the diffusion scheme need to satisfy the relationship shown in Eq. (1): where is the flow velocity, A is the section area, and B is the section width.
Daily adjustments in the unsteady flow occur in the river reach from Yibin to Chongqing. For the conditions of incoming flow, since the time interval Δ t = 3600 s cannot be changed, there are limitations in selecting the section area Δ x. For a stable diffusion scheme, when Δ t = 3600 s For most of the river reaches with dangerous beaches, the length of the channel that obstructs navigation is far shorter than 3.6 km. To study these segments, Δx should be less than 300 m, which causes serious instability in the explicit scheme. (1)

Implicit difference method
We discuss the applicability of the Preissmann implicit difference method as an example.
(1) The Preissmann implicit difference method calculates the conditions where water flow changes only slowly. The unsteady flow of the Yibin-Chongqing river reach changes drastically during the low-water period. In addition, the Preissmann implicit difference method always linearizes the data, which is not suitable for the conditions of this river reach. For example, for Q n j , As shown in Eq. (2), the Preissmann implicit difference method ignores ΔQ 2 j in linearization. In practice, the daily adjusted unsteady flow of the Yibin-Chongqing river reach during the low-water period generally changes at ± 1000 m 3 /s, and the change is sometimes greater than 1500 m 3 /s, accounting for over 50% of the total flow. In addition, according to the characteristics of the flow discharge of the power station, such a large change is usually completed within an hour, i.e., within a unit time interval Δt . During the low-water period, the flow Q often fluctuates between 2500 and 4000 m 3 /s. It is evident that the linearization in the Preissmann implicit difference method is not applicable to hydraulic conditions in which there are dramatic changes.
(2) The Preissmann implicit difference method is very unfavorable for program debugging. To simplify the method, the continuity and energy equations are simplified as shown in Eqs. (3) and (4). The equations show that both ΔQ j and Δz j are known fixed conditions and that A 1j , B 1j , C 1j , and D 1j have no obvious physical meanings. When calculation errors occur, it is not convenient to observe the changes in coefficients caused by various parameters, so debugging is almost impossible.
In conclusion, traditional difference schemes cannot completely meet the calculation conditions. The explicit scheme has inherent errors and may amplify errors due to instability; although the implicit scheme can unconditionally reach stability and calculate over a long duration, the dramatic changes in the unsteady flow of the Yibin-Chongqing river reach cause the Preissmann implicit difference method using the chasing method to generate large errors.

Proposed implicit iterative method
The Saint-Venant equations are transformed as shown in Eqs. (5) and (6): When R = 0, the solution meets the conditions of the Saint-Venant equations. In particular, x is the cross-sectional distance, t is the time step, Q is the flow rate, A is the cross-sectional area, h is the water depth, g is the acceleration of gravity, z is the water level, and K is the flow modulus.
This method refers to the four-point grid of the Preissmann implicit scheme and uses the weight coefficient = 1 , making the scheme fully implicit. The Preissmann format grid layout is shown in Fig. 3.
The continuous difference method form is shown in Eq. (7): Equation (8) shows the kinematic difference method form: In these equations, all f n j+1 ,f n+1 j ,f n j ,Δt and Δx are known. The unknowns include the area A n+1 j+1 of section j + 1 at moment n + 1, the flow modulus K n+1 j+1 , the water level z n+1 j+1 , and the flow rate Q n+1 j+1 . When the elevation of the section is known, as long as the water level z of a certain section is known, the area A and the flow modulus K of the section can be determined.
With the assignment of values to the unknown water level z n+1 j+1 , the value of the unknown flow rate Q n+1 j+1 can be obtained through continuity Eq. (1), and the corresponding area A n+1 j+1 and flow modulus K n+1 j+1 can be determined according to the terrain data. Then, we substitute these parameters into the kinematic equation (Eq. (8)) and set limit values n(n → 0) for R. When R < n, it is deemed to satisfy the kinematic equation (Eq. (8)), and the next step can begin; if R < n is not satisfied, the unknown water level z n+1 j+1 should be reassigned. Since there are multivariate unknowns in the kinematic equation of the Preissmann scheme and more than two solutions have been obtained, it is necessary to solve these problems through linearization. In the implicit iterative method mentioned herein, the unk nown water level z n+1 j+1 is replaced with z n+1 j+1 =z n j+1 + Δz j+1 , in which z n j+1 is known and Δz j+1 is unknown, representing the water level amplitude of section j + 1. When the condition R < n is satisfied but Δz j+1 takes an incorrect value, a limit value Rz can be set based on experience to recorrect Δz j+1 , and the limit value Rz is selected such that 1 < | | R z | | < 3 . In fact, the simulation reveals that the error value Δz j+1 in the multivariate solution is generally greater than 5 m. This restriction condition can prevent errors caused by linearization and ensure a correct calculation result.

Model accuracy verification
In Sect. 3.1, the advantages and disadvantages of the difference method for understanding the Saint-Venant equations are discussed in detail. The implicit iterative method proposed in this paper has the following advantages over the existing difference method: (1) Compared to the explicit difference method, the implicit difference method is characterized by few errors, stability, and few restrictions. The implicit difference method can be used to calculate the multi-wavelength unsteady flow in the Yibin-Chongqing river reach. (2) Instead of using the implicit difference chasing method, this paper employs an iterative calculation method to eliminate the errors resulting from the linearization of the chasing method. Moreover, owing to the calculation features, the iterative method can better debug the program, and the effect of unsteady flow propagation on various indicators can be observed.

Model accuracy verification (flume test)
Hu [28] studied the characteristics of unsteady flow propagation in a flume. Here, the test was conducted in a water tank 28 m in length, 0.56 m in width, and 0.7 m in height (Fig. 4). Eight ultrasonic water gauge probes were arranged in the water tank to synchronously measure the water level in real time so that the mode could be verified in three aspects: range of variation, outlet flow, and deformation features.

Flow change verification
The model is used to calculate the relative amplitude change and outlet amplitude of the outlet flow. Then, the calculated value is compared with the measured value. The results are shown in Fig. 5a and b. As shown in the figure, this model can clearly present the changes in the outlet flow amplitude.

Outlet flow verification
The working conditions for the verification are as follows: The inlet flow varies between 15 and 40 L, the gradient of the water tank is J = 0.3%, and the waveform with a is a sine wave with a change period of 20 s. The water-taking level with a steady waveform is located at section No. 6, which is 19.05 m away from the edge of the water tank. The moment when the wave trough appears is regarded as 0 s, i.e., the start-   Fig. 6. The average, maximum, and minimum flow rates under the verification working conditions are as follows: Q=27.5L∕s , Q max = 40L∕s,Q min = 15L∕s. As shown in Fig. 6, the calculated wave peak and trough in section No. 6 have very small errors compared to the measured values. The maximum error of 0.00278 m (a percentage error of 4%) appears at the wave trough where the period is 100 s.
As shown in Fig. 7, the maximum error in the flow at the wave peak is 0.769 L/s, and the maximum percentage error is 2.5%. The maximum calculated error in the flow at the wave trough is 2.306 L/s, and the maximum percentage error is 13.3%.

Verification of deformation features
The working conditions for the verification are as follows: The inlet flow changes between 15 and 40 L; the changing period is 10 s; the waveform is a sine wave; the following sections are used for verification: No. 2 (6.235 m from the edge of the water tank), No. 4 (12.64 m from the edge of the water tank), and No. 6 (19.05 m from the edge of the water tank). The working conditions for the calculation are as follows: The inlet flow changes between 15 and 40 L; the changing period is 10 s; and the waveform is a sine wave. The glass roughness factor is 0.0095, the time interval Δt is 0.2 s, and the gap between sections is 0.06 m. Figure 8 shows that the calculated results increase quickly and sharply and then decrease gradually and

Yibin-Chongqing river reach data verification
The working conditions for verification (1)  According to Figs. 9 through 11, the water levels, waveforms, and phases in the unsteady flow verification    (1) and (2) (all errors are absolute values.) Table 1 shows that the iterative method has a higher accuracy than the chasing method in terms of the mean water level value errors and in terms of the maximum value of errors at all verification sites.

Application of the model for rapid flow
For this example, the following conditions were used: A rectangular sectional channel with an overall length L = 100 m, width B = 10 m, flow rate Q = 20 m 3 /s, and Manning coefficient n = 0.03. The boundary conditions were as follows: both upstream and downstream inflows were rapid flows [29]. The calculation parameters were as follows: Δx = 10 m, Δt = 10 s, the total calculation time was 1000 s, and the time weight factor was θ = 1 [30] (Fig. 12).
In summary, the implicit iterative model proposed herein is more accurate in the process of simulating unsteady flow in long reaches and is thus appropriate for use in the unsteady propagation and calculation in long reaches.

Daily regulated unsteady flow verification
The river reach selected for calculation starts at the Xiangjiaba Hydropower Station along the Jinsha River and the Gaochang Hydrologic Station along the Minjiang River and extends to the mainstream of the Yangtze River in Fuling, Chongqing. Using the measured terrain data, we divide all zones into sections. The section spaces were 200-500 m, and the total number of sections was 1180. By verifying the unsteady flow spreading process in the 3 typical periods mentioned above (2 low-water periods and  Fig. 13, we know that the changing process of the simulated water level is essentially consistent with the changes in the measured water level at each station under the three operating conditions. Except for some deviation in the calculated value from the measured value at a few test points, the errors are mainly within ± 10 cm, which meets the calculation requirements.

Daily decrease in the amplitude and flow rate
After the hydropower station along the Jiansha River and Minjiang River began operation, the daily amplitude and maximum hourly amplitude of the daily adjusted unsteady flow in the Shuifu-Chongqing river reach are calculated and are detailed below and are shown in Fig. 15. As shown in Fig. 16, under typical operating conditions with a large average flow rate, the flow rate

Spreading of unsteady flow in the Yibin-Chongqing river reach
The spreading speed of unsteady flow is associated with the overall flow rate, and the flow velocity mainly depends on the flow rate. This section calculates the spreading process of unsteady flow for the three typical flow rates (normal-water, low-water, and flood periods) to evaluate the unsteady flow spread along this river reach. Figure 3 shows the changes in the flow rate along the Yibin-Fuling river reach at a certain moment under the three operating conditions. The flow rate in Yibin reaches a peak at this time. For convenience of comparison, a corresponding fixed value (4000 m 3 /s for normal water and 8000 m 3 /s for flood) is subtracted from the flow rates along the path for the typical operating conditions during the normal-water and flood periods.
We know that under low-water conditions, there are approximately 3. The spreading speed is mainly related to the overall unsteady flow velocity. Therefore, studies on the spreading time are primarily conducted by identifying the relationship between the spreading time and the overall flow rate.
The fitting formula is shown in Eq. (9): where T is the arrival time and Q is the flow rate in Yibin.
For a = × s × � BC √ Ri g , the unit is m 1.5 ·s 0.5 . In particular, s is the distance from a dangerous beach or station to Yibin, B is the width of the river, C is the Chezy coefficient, R is the hydraulic radius, I is the gradient, g is the acceleration of gravity, and α is the constant term. The value of b is generally approximately 0.4. For the decrease in the flow rate amplitude, the ratio of the flow rate amplitude at a dangerous beach to that in Yibin is used to measure the decrease in the flow rate amplitude at the designated dangerous beach.
The fitting formula is shown in Eq. (10): where ΔQ n is the flow rate amplitude at the dangerous beach or station, ΔQ yibin is the flow rate amplitude in Yibin, a is a proportion, and b denotes an exponential term. For the water level amplitude, the relationship between the flow rate amplitude and water level amplitude is identified.
The fitting formula is shown in Eq. (11): where Δh is the water level amplitude; ΔQ is the flow rate amplitude; a = BC √ Ri , where B is the river width, C is the Chezy coefficient, R is the hydraulic radius, i is the gradient, and α is the constant term; and b is the exponential term.
By summarizing the spreading time, the flow rate amplitude percent decrease, and the flow rate amplitude-water level amplitude relationship at the corresponding dangerous beach and important downstream stations (Luzhou, Hejiang, Zhutuo, Jiangjin, and Chongqing) and by fitting the data obtained, we obtain the empirical formulas for the spreading time, flow rate amplitude percent decrease, and the flow rate amplitude-water level amplitude relationship at important stations.
As shown in Fig. 17, at different flow rates, unsteady flow can reach Zhutuo within 18-40 h (corresponding to a flow rate of approximately 11,000-2250 m 3 /s). The average flow rate amplitude is 34.7% of the flow rate amplitude of Yibin. When the amplitudes are the same, there is essentially no difference in the changes during flood, normal-water, and low-water periods; when the amplitude is 0-750 m 3 /s, the water level amplitudes during flood, normal-water, and low-water periods are all 0-0.4 m. The changes in unsteady flow relationships are shown in Table 2.

Analysis of the incoming flow conditions calculated with the daily adjusted unsteady flow and designed water level
The designed minimum navigable water level (hereafter referred to as the "designed low water level") is the lowest water level used in designing waterway projects that allows standard ships or fleets to navigate the waterway. This value is the most basic design parameter in waterway projects, such as in channels and wharfs. Its value is directly related to project cost and project operation benefits. The daily adjustment of the Xiangjiaba Hydropower Station along the Jinsha River inevitably has some impact on the hydrological conditions of the downstream channel. This section mainly calculates and analyzes the designed water level given the daily adjusted unsteady flow.

Determination of the operating conditions for calculation
When calculating the designed water level of the Shuifu-Chongqing river reach, the natural conditions and daily adjustment conditions of the Xiangjiaba Hydropower Station should be considered. Thus, four types of operating conditions are provided to identify the operating conditions of the designed water level.

Operating condition 1
Natural daily average incoming flow (21 years from 1991 to 2011). Before the Xiangjiaba Reservoir was built, according to the data from the hydrological station in the 21 years from 1991 to 2011, with a 98% guaranteed rate of flow, we obtained the designed minimum navigable flow at each main station along the Xiangjiaba-Chongqing river reach. That of the Pingshan Station was 1245 m 3 /s (Table 3).

Operating condition 2
Operating conditions for calculating the minimum outflow rate of the Xiangjiaba Hydropower Station along the Jinsha River.
According to the scheduling method of the Xiangjiaba Hydropower Station, the minimum discharge flow rate after the Xiangjiaba Reservoir was built reaches 1518 m 3 /s, and the duration is more than 17-20 h. Therefore, this flow rate is selected for the calculation of the designed water level under steady flow mode (Table 4).

Operating conditions 3 and 4
Operating conditions for calculating the daily adjusted unsteady flow after the reservoirs along the Jinsha River and the Minjiang River were built.
After the Xiangjiaba Reservoir was built, the minimum discharge flow rate was 1518 m 3 /s. Although the duration exceeds 17-20 h, the unsteady flows discharged from the Xiangjiaba Hydropower Station are flattened during the spreading process. The minimum flow rate increases minimally when the flows reach the Luzhou and Zhutuo stations. Therefore, operating condition 3 uses this flow rate to calculate the designed minimum water level under unsteady flow mode, and operating condition 4 calculates the daily average water level under unsteady flow mode.

Numerical results
As shown in Fig. 18, for operating condition 2, when the Xiangjiaba Hydropower Station discharges at the base flow rate (1518 m 3 /s), the designed water level of the Yibin-Chongqing river reach is approximately the same as that under natural conditions, although there is a slight increase of 0.03-0.07 m. Under operating condition 3, at the minimum unsteady discharge flow rate, the designed water level of the Yibin-Chongqing river reach increases by 0.11-0.25 m, approximately 0.19 m on average. Under operating condition 4, at the daily average unsteady flow rate, the designed water level of the Shuifu-Chongqing river reach is 0.17-0.32 m higher than that under natural conditions, averaging approximately 0.25 m.

Analysis of the changes in the water level of the Chongqing-Yibin river reach before and after sand excavation
Since the flow rate of the gauging stations along the path changes extensively during the year, the water level varies with the upstream water level scheduling process, river reach, and different flow rates during low-water, normalwater, and flood periods. Therefore, three flow rates are selected for calculation at Zhutuo Station during the The flood surface lines during low-water, normal-water, and flood periods after sand excavation are shown in Fig. 21. At the three flow rates, the changes in the Xiangjiaba-Fuling water surface line are essentially the same. The upstream water face line has a larger gradient, while the downstream water face line features a smaller gradient.
As shown in Fig. 22, the changes in the water level reduction amplitudes at the three flow rates are similar, and the water level reduction amplitudes at different flow levels are as follows: low-water period > normal-water period > flood period. However, for some river reaches, the water level reduction values in the normal-water period are larger than those in the flood period. This difference is mainly the result of the presence of a large number of sand pits in the point bar where the water cannot flow out during the low-water period, and during normal-water and flood periods when the high water level is high, the water level reduction amplitude increases due to the impact of the sand pits. Furthermore, the upstream side has more river reaches, with a water level reduction value more than 1 m greater than that of the downstream side. Table 3 shows that the overall water level at the three flow rates decreases by 0-2 m. In particular, the decreases at 700-780 km and 900-1050 km are relatively large at 1-2 m. From 552 to 600 km, which is the Three Gorges Reservoir region, there is a small water level reduction value, i.e., 0.2-0.5 m. The results of comparing the water level reduction value and the amount of sand excavation during the sand excavation along the Fuling-Xiangjiaba reach in 2016 are shown in Fig. 13.
The sand excavation change characteristics along the path are generally similar to the water level reduction amplitude. The amount of sand excavation gradually increases at 550-670 km, decreases at 670-850 km, suddenly rises at 900 km, and then declines at 900-1076 km. Likewise, the water level reduction amplitude increases abruptly at 550-700 km, gradually decreases at 700-860 km, suddenly increases at 860-900 km, and then gradually decreases. The change trend is almost the same as that of sand excavation. Due to improper sand excavation, the distance between each sand pit and the amount of sand excavation differs extensively. Therefore, the average value of the sand excavation volume and water level reduction amplitude within every 50 km is summarized to analyze the effect of sand excavation on the reduction in the water level. The sand excavation conditions and water level changes along the path are shown in Figs. 23 and 24.
According to Table 5, the ratio mainly falls within the range of 1-3 m, and the ranking is still as follows: flood < normal water < low water, which indicates that during the low-water period, sand excavation has the greatest impact on the water level. Additionally, the upstream sand excavation has a larger influence on the water level reduction amplitude. In particular, at 826-676 km, 10 million m 3 of sand is removed every 50 km, which causes the water level to drop by 0.2-0.3 m during the low-water period.
According to the relationship between the amount of sand excavation per 50 km and the change in the water level (Fig. 23), with the increase in sand excavation, the water level reduction amplitude grows. In particular, when the Xiangjiaba-Zhutuo and Zhutuo-Fuling river reaches have 10 million m 3 of sand excavated per 50 km, the water level decreases by 0.6 m and 0.33, respectively, mainly because the river channel becomes increasingly narrow on the upstream side. A narrower river reach results in a smaller flow area at the same flow rate and thus a larger impact of sand excavation on the change in the water level.

Results and discussion
The upstream portion of the Yangtze River along the Yibin-Chongqing river reach will undergo a waterway regulation project. The waterway level in this section will be upgraded from Class III to Class II. The daily regulation of unsteady flow and sand mining in this upstream section will have a great impact on the water level of this reach, Fig. 24 The relationship between the amount of sand extracted and the water level reduction every 50 km and changes in the water flow conditions will directly affect the implementation of the remediation project. During a low-water period, we reviewed the typical water gauges (at Lizhuang, Yujiawan, Shenbeizui, Zhutuo, Zhengjialiang, and Cuntan) in the section of the Yangtze River from Chongqing to Yibin (Figs. 25, 26). According to the principle of an approximately equal gradient in the low-water period, the water level was determined when the design flow rate was Qzhutuo≈2230 m 3 /s and Qcun-tan≈2900 m 3 /s. An overview of the results is shown in Table 6.
The model calculation results given the design flow rate show that the average decreases in the water level corresponding to the 20-km sand extraction near the   [31], Xiao [20] and the Water Resources Committee of the Yangtze River, the water levels at the Cuntan and Zhutuo stations dropped by approximately 0.15 m and 0.30 m, respectively. The propagation characteristics of unsteady flow and the influence of unsteady flow and sand mining on the design water level calculated by the model in this paper can provide technical support for the waterway regulation of this section of the river. The relationship between the water level and discharge at the Zhutuo hydrological station after sand mining is shown in Fig. 27.
In this paper, only unsteady flow and sand mining were calculated and analyzed under the current channel boundary conditions. After the cascade hub is completed and jointly scheduled upstream of the Jinsha River, the hydrological conditions will change further. The evolution and regulation of this section are a subject that we will focus on in subsequent work.

Conclusion
(1) The new implicit iterative method proposed in this study is much better than the chasing method. Suitable for the calculation of the spreading of multiwavelength and long-distance unsteady flow, this method can be used to calculate the multiwavelength unsteady flow of the Yibin-Chongqing river reach. The verification with measured data shows good results. (2) Through fitting analysis of the spreading time, flow rate amplitude percent decrease, and flow rate amplitude-water level amplitude relationship at important stations (Luzhou, Hejiang, Zhutuo, Jiangjin, and Chongqing), we obtained empirical formulas for the spreading time, flow rate amplitude percent decrease, and flow rate amplitude-water level amplitude relationship at these important stations.

Conflict of interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.