Coupled modeling of bed deformation and stability analysis on a typical vulnerable zone on the Ningxia reach of the Yellow River during flood season

The Ningxia reach of the Yellow River (NRYR) is wide and shallow, and instability and damages are likely to occur at its vulnerable zone. In this study, based on the one‐dimensional hydro‐sediment dynamics, two‐dimensional hydro‐sediment dynamics and vulnerable zone stability analysis model, a coupled model was established to analyze the scouring risk of the vulnerable zones in the NRYR during flood season. Considering the morphological change of bank slope at the vulnerable zones caused by flow erosion, the seepage stability of a typical spur dyke and protection embankment under the water level fluctuation conditions was analyzed. To fully reflect the coupled scouring effect of flow and sediment, the water flow scouring intensity index α is proposed. A significant increase of water flow scouring intensity index α will result in bank toe scouring, while its significant decrease will lead to silting of the bank toe. The probability of instability failure of the bank for the spur dyke slope is greater during the flood recession period. Compared with the influence of the change of river water level in the flood season, the slip stability of the bank slope is more sensitive to the change of bank toe erosion.


| INTRODUCTION
The Ningxia reach of the Yellow River (NRYR) has long been dynamic in channel morphology. Damming and diversion are important human activities that affect the water and sediment conditions of the NRYR (Fan, Shi, Zhou, & Du, 2012;Fan, Wang, & Ran, 2010). The river channel possesses complex morphology, large mainstream swing, and erodible riverbanks in the NRYR (Tian, Ma, Yuan, Wang, & Yue, 2019;Yuan, Tian, Wang, Liu, & Chen, 2018). Flood disasters have occurred frequently throughout its history. Since 1904, there have been five catastrophic floods with a peak discharge exceeding 5,000 m 3 /s. In recent years, the construction and operation of many large-scale water conservancy projects, as well as the protection and restoration of vegetation in the NRYR, have greatly altered the regional mechanism of sediment yield, thereby resulting in significant change in the water and sediment inflow process and its situation of bed deformation in the NRYR (Zhang et al., 2019). A section of dike that may be dangerous under actual conditions affected by warping river, soft soil, or violent flow is called the vulnerable zone. During the flood season, the near-shore part of the vulnerable zone is prone to scouring, thus resulting in steep ridges along the bank slope. For example, the maximum scouring depth of the vulnerable zone during the flood period of 2012 in NRYR is 4-5 m (Wang & Yue, 2015). This in turn threatens the stability of the bank slope, and seriously affects the flood control safety of the main dike of the Yellow River and its infrastructure on the river bank (Ta, Xiao, & Dong, 2008). River regulation projects such as a spur dyke and protection embankment were built on the vulnerable zone to protect the dykes, followed by river regime control and diversion (Papanicolaou, Elhakeem, & Wardman, 2011).
It is important to study the impact of flow scour on the stability of vulnerable zone bank slopes (Chiang, Tsai, & Yang, 2011;Stecca, Measures, & Hicks, 2017). The bank slope stability state under the scouring of river flow mainly depends on two factors: one is the bank slope section and filling material, while the other is the response of the bank slope under the load effect, such as the scour force and seepage pressure (Enlow et al., 2018;Nicholas, 2013).
Bed deformation is typically a key factor affecting river bank slope stability (Karmaker & Dutta, 2015). In recent decades, the morphodynamics model, in particular the 2D water and sediment mathematical model, has been used to study the evolution of alluvial channels (Ahilan, Guan, Sleigh, Wright, & Chang, 2018;Darby, Rinaldi, & Dapporto, 2007;Ham & Church, 2012;Lai et al., 2015;Langendoen et al., 2016;Nakagawa, Zhang, Baba, Kawaike, & Teraguchi, 2013;Rousseau, van de Wiel, & Biron, 2017;Tal & Paola, 2010). Quantitative research on the erosion rate of the bank slope faces difficulties in terms of two aspects: First, it is difficult to balance the solution time and accuracy of the complex motion model for large-scale river bed deformation (Pittaluga & Seminara, 2011); and second, it is difficult to determine the impact of river regulation projects such as the spur dyke and protection embankment on the erosion of the bank slope.
Seepage effect is one of the important factors affecting the stability of the bank slope. Seepage results in groundwater level change in the soil, which in turn affects the physical and mechanical properties of bank slope soil, such as pore water pressure (Akay, Ozer, Fox, & Wilson, 2018;Fox & Wilson, 2010;Karmaker & Dutta, 2015). For example, Fox and Felice (2014) studied the influence mechanism of groundwater seepage on bank stability. Nardi, Rinaldi, and Solari (2012) used model experiments to analyze the influence of seepage erosion and pore water pressure on bank slope stability. Scholars have studied the interaction among the river bed deformation, seepage change and bank slope instability. Rinaldi, Mengoni, Luppi, Darby, and Mosselman (2008) combined a model of fluvial bank erosion with groundwater flow and bank stability analyses to account for the influence of hydraulic erosion on mass failure processes. Rousseau et al. (2017) combined the genetic algorithm-based bank analysis model with the Telemac-Mascare two-dimensional water-sediment model, to analyze the phreatic line position in the bank slope and assess the bank slope stability. Deng, Xia, Zhou, and Lin (2019) proposed a coupled model for simulating bed deformation and bank erosion. The plane slide and collapse were calculated according to the method of different dual structural riverbanks. However, the existing mathematical models which can be applied for water and sediment transport are seldom able to simulate the coupled process of channel water and sediment transport, local bed deformation of vulnerable zones, lateral erosion of the bank toe, and bank slope instability. They also seldom consider the influence of different factors, such as riverbank soil properties and seepage changes, on the bank slope instability.
Based on the one-dimensional river water and sediment dynamics model, the two-dimensional river water and sediment dynamics model, and the vulnerable zone stability analysis model, the coupled model is established to analyze the scour safety of the vulnerable zones in NRYR during flood season. This study considers the influence of lateral scour of the bank toe at a vulnerable zone, near-shore bed deformation and seepage on the bank slope stability of the vulnerable zone, and greatly improves the calculation efficiency while meeting the simulation accuracy conditions. The research results provide technical support for the fine simulation of the vulnerable zone and river flood evaluation in the flood season.

| Brief introduction to the research reaches
The NRYR has a total length of 397 km and an average water depth of $8-10 m. The water and sediment have different sources, the amount of water and sediment is unevenly distributed during the year, and varies greatly from year to year. The NRYR can be divided into the plain section, reservoir section, and canyon section, of which only the canyon section consists of rocky river, with the remainder of the river section being sandygravely river or sandy river (Yuan, Yue, Tian, Cao, & Song, 2020). The reach from Qingtongxia to Mianjiawan was selected as the study area to be covered in this article. The reach is located in the Qingshi section of the NRYR, and has a total length of 10.2 km. The water surface is wide and shallow, and the reach has a main channel width of 400-700 m, average gradient of 0.62‰, and sinuousity of 1.16. This reach is a sandy-gravel braided river with poor scour resistance, and the riverbed deformations have changed significantly over the years. The study area includes two hydrological stations, respectively, in Qingtongxia and Shizuishan, and five stage gauging stations in Yesheng, Shiba, Gaorenzhen, Wuduizi, and Hongyazi. The location of the study area is shown in Figure 1. A large number of bank revetment projects and river regulation projects have been built in the NRYR, and the change of the channel morphology is limited to some extent. However, bank collapse caused by floods still occurs frequently, thus resulting in the destruction of embankments, ditches, roads and diversion and drainage outlets. The one-dimensional water and sediment model modeling area was from Qingtongxia to Shizuishan, and the two-dimensional water and sediment model modeling area was from Qingtongxia to Mianjiawan.

| Structural composition and physical and mechanical properties of the vulnerable zone
After the construction of the second-phase flood control project in the NRYR, a vulnerable zone project is covered in the modeling area: the vulnerable zone project of the Caijiahekou, which includes 13 protection embankments and two spur dykes, with each spur dyke having a length of 90 m and average spacing of 100 m. According to the results of an engineering geological investigation (Zhou & Wang, 2009), the vulnerable zone bank slope soil in the study area can be divided into two layers: the upper layer consists of artificial fill (mainly sandy loam soil), with a thickness of 2 m, while the lower layer consists of a fine sand and gravel mixed layer (mainly fine sand), showing a medium dense to dense state. The physical and mechanical parameters of the soil layers are shown in Table 1.

| Mathematical model
The calculation process of the flood scour stability assessment model for the vulnerable zone is shown in Figure 2. First, the water and sediment transport and bed deformation are calculated using the one-or two-dimensional hydro-sediment dynamics model. Second, the lateral erosion distance of the vulnerable zone bank toe is analyzed and calculated. Third, the groundwater level variation in the typical vulnerable zone soil is calculated by taking the water level in the river as the boundary condition. Finally, based on the calculation results of the river level, the lateral erosion distance of the bank toe and groundwater level change, the bank slope stability degree of typical vulnerable zone is analyzed, and whether the vulnerable zone will undergo a dangerous situation is determined.

| Selecting the stream channel model
The one-and two-dimensional calculations of the water area near the vulnerable zone are carried out in order to F I G U R E 1 Study area with hydrological station and water level station realize the combination of large-scale calculation and small-scale accurate simulation, thereby reducing computing time and improving the response speed of the model. The water flow, sediment transport and bed deformation are calculated by the one-dimensional water and sediment dynamics model. The governing equation of river flow is the Saint-Venant system of equations. In order to reflect the characteristics of sediment transport of the NRYR, the sediment transport governing equation adopted in this article is the sediment convectiondiffusion equation (Papanicolaou, Elhakeem, Krallis, Prakash, & Edinger, 2008;Runkel, 1998), and the bed load transport rate is calculated by using the Engelund and Hansensediment transport rate equation (Brignoli, Espa, & Batalla, 2017). The upstream boundary conditions of the one-dimensional water and sediment mathematical model are the cross-section flood discharge hydrograph and sediment concentration hydrograph, and the downstream boundary is the cross-section water level-discharge relationship. Suspended sediment transport accounts for more than 90% of the river's total sediment transport in the NRYR (Long & Zhang, 2002), thus a suspended sediment concentration hydrograph was used to represent the sediment concentration hydrograph.
The local bed deformation near the vulnerable zone is calculated using a two-dimensional water and sediment dynamics model (Wu, Jiang, & Wang, 2004). The basic governing equations of the hydrodynamic module are as follows: where t is the time, x and y are the horizontal Cartesian coordinates, u and v are the components of the flow velocity in the x and y direction, h is the flow depth, z s is still water surface elevation, g is the acceleration of gravity, ρ 0 is water density, I 0 is the magnitude of discharge due to point source, u s and v s are the velocities by which the water is discharged into the ambient water, u and v are the average flow velocities along the depth of the water, T xx , T xy , T yx , and T yy are the depth averaged turbulent stresses, Smagorinsky's equation is used for the lateral the turbulent stresses (Smagorinsky, 1963). τ bx and τ by are the x and y components of the surface and bottom stresses, and the Manning's equation is used to express bottom stresses (Wu et al., 2004).
T A B L E 1 Soil layer physical and mechanical parameters of the bank slope in the vulnerable zone Note: ω indicates water content, ρ 0 indicates natural density, e indicates porosity ratio, c indicates cohesion of soil, φ indicates internal friction angle, and K indicates permeability coefficient.
The calculation process of the flood scour stability assessment model for the vulnerable zone The basic governing equation of the sediment module is the sediment transport and diffusion equation in the following form (Xia, Lin, Falconer, & Wang, 2010): where S is suspended sediment concentration (SSC), p mod and q mod are the components of the flux correction value in the x and y direction, ε x and ε y are the turbulent diffusion coefficients of sediment in the x, y direction, α * is the coefficient of saturation recovery, ω * is the sediment settling velocity, and S * is the suspended sediment transport capacity.
The sediment transport capacity is calculated using the van-Rijn formula (van Rijn, 1984a, van Rijn, 1984b, and the equation of suspended sediment transport capacity is as follows: where f s represents the suspended sediment transport correction factor, c a represents the suspended sediment concentration nearby the bed surface, and V represents the water flow velocity. The bed load sediment transport capacity is shown in previous relevant research (van Rijn, 1984b). The upstream boundary conditions of the twodimensional water and sediment mathematical model are the flood discharge hydrograph and sediment concentration hydrograph of the Qingtongxia Hydrological Station, which is the nearest hydrological station upstream of the vulnerable zone. In addition, the downstream boundary conditions are the downstream water level and sediment concentration hydrograph calculated by the one-dimensional water and sediment model. Due to the fact that one-dimensional models cannot provide detailed lateral distributions of bed deformation, the bed deformation computed by the one-and two-dimensional models were not linked or compared in the vulnerable zone. This article focuses instead on the bed deformations computed by the two-dimensional water and sediment model of the riverbed surrounding the vulnerable zone project. A Finite volume high resolution scheme is used to solve the two-dimensional mathematical model of the hydrodynamics and sediment movement, based on the Roe Approximate Riemann solver (Bladé et al., 2012;Boscheri, Balsara, & Dumbser, 2014;Roe, 1981). Triangular unstructured meshes are used to partition the research area, and local refinement is performed for the meshes as the attachment to the spur dyke and protection embankment. The edge length of the main channel mesh is 20-40 m, that of the flood plain mesh is 50-70 m, and that of the mesh around the spur dyke and protection embankment is 5-10 m. Compared with the incoming sediment in the upper reaches of the river, the proportion of the sediment generated from bank toe erosion was very small. Therefore, sediment generated from the bank toe erosion was not incorporated in the two-dimensional model as a source-sink term.

| Bank toe erosion model
The improved Osman lateral erosion model is used to calculate the lateral erosion distance of the bank toe. The Osman lateral erosion calculation model is suitable for the clay bank slope (Osman Akode & Thorne Colin, 1988). In order to calculate the lateral erosion characteristics of fine sand under the action of flow scour, the influence factor m is added to the Osman formula to modify the lateral erosion rate of the bank toe under the water flow scour effect of the bank slope (Jiang, Zhu, Shen, & Wang, 2015), and the lateral erosion distance of the bank toe in Δt time is as follows: where ΔB is the erosion distance of bank slope soil caused by water flow scouring within the time Δt, m is the influence factor m = 2, f 1 is the lateral scour factor, τ sw is the scour force of water flow, τ c is the critical scour resistance of bank slope soil, γ is the unit weight of bank slope soil, and J is the hydraulic slope. A Shields curve is widely used in calculating the critical shear force of a noncohesive soil bank slope (Paphitis, 2001). The Shields curve is the relationship curve between the dimensionless critical relative drag force and particle Reynolds number drawn by Shields according to the experimental results. The formula for calculating the impact resistance force of the soil body is as follows: where θ cr indicates the Shields number, ρ s indicates the sediment density, and D indicates the soil representative diameter.

| Bank slope seepage-stability analysis model of vulnerable zone
In this study, the seepage process in the bank slope soil was simulated using the saturated-unsaturated seepage finite element method. The basic equation of twodimensional seepage is as follows (Lam, Fredlund, & Barbour, 1987): where h w (x, y, z) indicates the function of water head, k x and k y indicate the permeability coefficient for the x and y directions, and C indicates the water capacity, C = ∂θ ∂h , is the volumetric water content.
Next, the effective stress method is used to calculate the shear strength of the bank slope soil. The relationship between the permeability coefficient and matrix suction can be obtained using the method proposed by Fredlund, Xing, and Huang (1994), in combination with the saturated permeability coefficient of the soil. The equation is as follows: where k w is the unsaturated permeability coefficient of the soil, k s is the saturated permeability coefficient of the soil, ψ is the matric suction of the soil, θ ' is the derivative of θ, ψ aev is the air-entry value of the soil under consideration, b = ln (1000000), e is the natural number, 2.71828, and y is a dummy variable of integration representing the logarithm of suction. Based on the stream channel model and bank toe erosion model, the morphological change of the vulnerable zone bank slope is determined. Next, the soil moisture distribution, an important parameter affecting bank slope stability, is determined using the saturated-unsaturated seepage finite element method. On the basis of seepage finite element simulation, the Morgenstern-Price method is used to perform the bank slope stability analysis of the vulnerable zone, and the variation of the typical bank slope stability with time in the flood season is studied (Zhu, Lee, Qian, & Chen, 2005). In the Morgenstern-Price method, the optional inter-slice force functions include constant function, semi-sinusoidal function, peak-clipping sinusoidal function, trapezoidal function, and discrete point data. The relationship between shear force and normal force between soil slices can be expressed as follows: where X represents the shear force among soil slices, λ represents the weight of the inter-slice force function, N represents the normal force among soil slices, and f (x) represents the inter-slice force function. Two types of safety factors are defined according to the moment balance and force balance. Then, based on Equation (11), the relationship between shear force and normal force among soil strips is established. After the iterative process has been completed, the ratio between the shear force and the normal force among soil strips is changed until the two safety factors are the same, and the safety factor which satisfies both the moment balance and force balance conditions is the bank slope stability safety factor K s . The torque balance safety factor and force balance safety factor are shown as follows (Fredlund & Krahn, 1977): where, F m represents the safety factor of moment balance, F f represents the safety factor of force balance, c 0 represents the effective cohesion, ϕ 0 represents the effective friction angle, u represents the pore water pressure, W represents the soil weight, D p represents the concentrated point load; β 0 , R 0 , x 0 , f 0 , and d 0 represent the geometric parameters, and α 0 represents the inclination of the soil bottom.

| Scour intensity index
The results of previous studies have shown that the water-sediment relationship is the main external factor affecting the bed deformation (Purvis & Fox, 2016;Xia et al., 2016). In the present article in order to study the effects of flow velocity and sediment concentration on the local scour at the bank toe of the vulnerable zone in the flood season, the flow velocity is used to represent the incoming water condition, while the sediment concentration is used to represent the incoming sediment condition. Based on the relationship between water and sediment, physical properties of water and sediment, and boundary conditions, a scouring intensity index α of water flow is proposed in this article as follows: where ω is the average sediment settling velocity, a is dimensionless coefficient, the calculated value is 10 −5 . In an alluvial river, the factors that affect the sediment transport capacity of water include flow conditions, physical properties of water and sediment, and boundary conditions (Cao, Pan, Wang, Li, & Zhang, 2018). In order to study the effects of flow velocity on the local scour at the bank toe, V 2 /gh or Froude number is used to approximately represent flow intensity (Stelling & Duinmeijer, 2003), V =ω is used to approximately represent relative gravity settlement (Cao et al., 2018), aρ 0 is used to approximately represent the physical property of water and sediment, aρ 0 V 3 =ghω approximately represents the sediment transport capacity at a given river regulation project, and α approximately represents the ratio of sediment transport capacity to incoming sediment concentration under a given velocity.

| Model calibration
Based on the water level and suspended sediment concentration hydrograph of the Qingshi reach during the typical flood period in 2012, the parameters of the model are calibrated. The time step of the one-dimensional water and sediment model is 10s. The time step of the two-dimensional water and sediment model is 15 s, which is chosen to satisfy the Courant-Friedrichs-Lewy (CFL) condition for numerical stability (Hou, Liang, Simons, & Hinkelmann, 2013). The time step of the twodimensional seepage model is 1 hr. The Manning's n roughness values of the main channel are 0.020-0.030, and those of the flood plain are 0.030-0.038, which are fixed with the amount of sediment transport and the evolution of bed conditions. The diffusion coefficients in the x and y directions are both 0.028, the coefficient of saturation recovery under siltation conditions is 0.25, and the coefficient of saturation recovery under scouring conditions is 1.0. There are a total of 44 fixed sections in the calculation reach, two hydrological stations in Qingtongxia and Shizuishan, and five gauging stations in Yesheng, Shiba, Gaorenzhen, Wuduizi, and Hongyazi. The topographic data of the one-dimensional watersediment model of the river channel are taken from the cross-section topography measured in June, 2011, and the topographic data of the two-dimensional watersediment model of the river are the DEM data measured in June, 2011. The inlet boundary conditions are the measured discharge and the sediment concentration hydrograph of the Qingtongxia Hydrological Station for 2012. The lateral inflow conditions are the measured incoming water and sediment data of the Kushui River, and the outlet boundary conditions are the actual measured water level data of the Shizuishan Hydrological Station.

| Calculation scheme
The model calculation scheme is shown in Figure 3. Recently, the years of major floods in NRYR were 1981 and 2012 according to statistics. Although the measured maximum peak flow of each hydrological station in 1981 was greater than that in 2012, the channel morphology. Damming and diversion have altered greatly (Yuan, 2014). The actual measured flood of the Qingtongxia Hydrological Station in 2012 is selected as the typical flood, while the discharge hydrograph corresponding to the design flood with a 20-year recurrence period is deduced by using the homogeneous multiple enlargement method. Based on the actual measured data of the flood season in 1989, 1990, 2007, 2010, 2011, and 2012, the power function empirical relationship is then selected to establish the relationship between the discharge and sediment transport rate of the Qingtongxia Hydrological Station in the flood season, and to deduce the sediment concentration hydrograph corresponding to the design flood with a 20-year recurrence period.

| Model calibration
A comparison of the measured values and calculated values of the highest water level in the NRYR in 2012 is shown in Figure 4. The verification results of the water level and suspended sediment concentration are shown in Figure 5. The maximum difference between the F I G U R E 3 The design flood with a 20-year recurrence period corresponds to the discharge and sediment concentration hydrograph measured water level and the water level calculated from the model at the five water level stations, including Yesheng, Shiba, Gaorenzhen, Wuduizi, and Hongyaze is <20 cm, and the maximum difference between the measured water level and the water level calculated from the model in each cross-section is <20 cm as well. Furthermore, the maximum ratio of errors and mean flow depths in each cross-section is <5%. The sediment concentration hydrograph measured at the Shizuishan Hydrological Station is consistent with the calculated sediment concentration hydrograph, and the maximum difference between the two hydrographes is <10%.
To verify the calculation accuracy of the model in the NRYR, the riverbank inundation scope, which was determined according to the aerial images ( Figure 6) after the flooding in 2012, was analyzed based on a comparison with the calculation results. The established hydrosediment dynamics model can more accurately simulate the process of water and sediment transport in the NRYR during the flood season.

| Analysis of velocity and bed deformation
Based on the conditions of the designed flood, the maximum flow velocity distribution and the change in the river bed deformation in the reach of the vulnerable zone  Figures 7 and 8, respectively. The river shows the complex bed deformation state after the flood season, and the main channel mainly shows scour, whereas the flood plain mainly shows silt. By analyzing the overall change in the bed deformation, it was observed that the riverbed reached a maximum velocity of 3.25 m/s, a scour depth of $3.6 m, and a maximum silt thickness of $1.2 m. By analyzing the flow velocity distribution and bed deformation of the vulnerable zone at Caijiahekou, it was discovered that the river bed of the #12, #13, and #15 spur dykes and the protection embankments at Caijiahekou were affected by a greater flow velocity and scour depth, with the maximum flow velocity at the bank toe being >2.25 m/s, and the scour depth at the bank toe >1.0 m. The #12, #13, and #15 spur dykes and protection embankments at Caijiahekou were selected as the research objects in this study, by which changes in scour depth and the lateral erosion distance of bank toe in a vulnerable zone were observed before and after the flood.
The changes in scour depth at the bank toe of the spur dyke and protection embankment during the flood season are as shown in Figure 9. Based on the characteristics of water level fluctuation, the flood season can be categorized into the water-rising, flood peak and waterfalling periods, which spanned 1-28, 29-37, and 38-58 days, respectively. The change rule of the scour depth of the spur dyke and protection embankment over time is as follows: the scour depth of the spur dyke and protection embankment in the water-rising period increases gradually; and the change amplitude of the scour depth decreases gradually. The change amplitude of the scour depth of the spur dyke and protection embankment is smaller in the flood peak and waterfalling periods, whereas the scour depth of the spur dyke and protection embankment after the flood season is balanced. In addition, the #15 spur dyke at Caijiahekou has the greatest scour depth, that is, at up to 2.98 m. Based on the designed flood, the bank toe of the #15 spur dyke at Caijiahekou shows obvious alternating scouring and silting. In terms of the overall change trend, the scour depth of the bank toe varied logarithmically, and conformed to the change rule of the local scour depth of the spur dyke.
Next, based on the #15 spur dyke at Caijiahekou as the research object, the causes of scouring and the silting alternation phenomenon were analyzed. The change in the near-shore water flow scour intensity of the #15 spur dyke at Caijiahekou over time was calculated, as shown in Figure 10. When analyzing the relationship between the bank toe deformation and water-sediment conditions, the hysteresis effect between the two must be considered. The hysteresis effect refers to the fact that the bank toe deformation of an alluvial river can be constantly regulated by the change in external conditions such as the inflow water and sediment, and the regulation of bed deformation involves a certain hysteresis time. Furthermore, the response and regulation times differed owing to the complexity of the change in external conditions (Wu, Xia, Fu, Zhang, & Wang, 2008). The analysis  Figure 10 indicates that a significant increase in water flow scour intensity index α will result in scours at the bank toe, whereas a significant decrease will result in silts at the bank toe. In the 30-32 days range, α decreased, silt occurred at days 32-33, and the hysteresis time after silt regulation was 1-2 days. In the 32-34 days range, α increased, the resulting scour occurred at days 33-38, and the hysteresis time after scour regulation was 1-4 days. In the 34-38 days range, the α decreased, silt occurred at days 38-42, and the hysteresis time after silt regulation was 4 days. In the 40-42 days range, the α increased, scour occurred at days 42-48, and the hysteresis time after scour regulation was 2-6 days. After 42 days, α indicated fluctuations without significant changes, and the bank toe deformation basically remained stable.
Compared with other calculation methods applied to the hydrodynamic force and riverbed deformation (Deng et al., 2019;Nicholas, 2013;Rousseau et al., 2017), the model described herein nests the local two-dimensional water-sediment dynamic model inside a one-dimensional river channel water-sediment dynamic model and focuses on the hydrodynamic force and sediment movement near the river section of a vulnerable zone. This one-and twodimensional nested model can significantly improve the calculation efficiency greatly, provided that the simulation accuracy conditions are satisfied.

| Analysis of lateral adjustments
Under the condition of a flood designed with a 20-year recurrence period, the calculation results of the nearshore average velocity, average shear stress and lateral erosion distance of the #12, #13, and #15 spur dykes and protection embankments at Caijiahekou are as shown in Table 2. The shear stresses of the average near-shore water flow of the #12, #13, and #15 spur dykes and protection embankments at Caijiahekou during the flood season are 10.3, 12.9, and 13.1 pa, respectively. The average flow velocity of each spur dyke and protection embankment in the flood season was >2 m/s, and the lateral erosion distance was >1 m. The average flow velocity of the #15 spur dyke at Caijiahekou was 2.50 m/s during the flood season, and the lateral erosion distance at the bank toe was 1.78 m. The average flow velocity and lateral erosion distance of the #15 spur dyke at Caijiahekou in the flood season were the greatest.
The change process of flow velocity and shear stress of the #15 spur dyke at Caijiahekou during the flood season is shown in Figure 11. The change scope of shear stress of the near-shore water flow during the flood season was within the range of 10.2-19.5 pa, and that of the near-shore water flow was always greater than the critical shear stress of the soil body (4.47 pa). In addition, F I G U R E 9 Changes in scour depth at bank toe of spur dyke and protection embankment during the flood season F I G U R E 1 0 Changes in near-shore water flow scour intensity of the #15 spur dyke at Caijiahekou over time. α indicates water flow scour intensity index the bank toe was always in an erosion state, and the lateral erosion reached its maximum value at the end of the flood. In the early stage of the water-rising period (0-5 days), the near-shore velocity was smaller and the shear stress was greater, owing to the lower water level and the smaller riverbed scour depth. With an increase in the channel water level and riverbed scour depth, the shear stress stabilized gradually. During the flood peak period, the near-shore flow velocity increased significantly, the shear stress reached a higher peak value, and the bank toe erosion rate was significantly higher than those in the water-rising and water-falling periods. With the decrease in flow velocity, the nearshore shear stress in the water-falling period decreased significantly, and reached a minimum value at the end of the flood.

| Stability analysis of typical Bank slopes in flood season
The bank slope of the #15 spur dyke at Caijiahekou was affected significantly by scour erosion. The computational cross-section and coordinate system are shown in Figure 12. It is regarded as the research object to calculate the stability safety factors K s1 and K s2 , where K s1 indicate the stability safety factors considering and not considering scour erosion, respectively, as shown in Figure 13. When the scour erosion of the bank toe was considered, the river bank stability safety factor K s1 exhibited a certain fluctuation with the change in the river channel water level but decreased gradually under the scour effect of the near-shore water flow. The change scope of K s1 in the water-rising period was 1.067-1.629, T A B L E 2 Calculation results of near-shore average velocity, average shear stress and lateral erosion distance of #12, #13, and #15 spur dyke and protection embankments at Caijiahekou F I G U R E 1 1 Change process of flow velocity and shear stress of #15 spur dyke at Caijiahekou duringflood season F I G U R E 1 2 Computational crosssection of seepage and stability model the average value was 1.297, and the bank slope stability was relatively high. K s1 varied from 1.087 to 1.152 in the flood peak period, with an average value of 1.118, which is lower than that in the water-rising period. The K s1 varies from 0.934 to 1.129 in the water-falling period, with an average value of 1.011. During the water-falling period, the water level in the river channel decreased gradually, whereas the lateral water pressure decreased, the change in the groundwater level slightly lagged behind the water level in the river channel, and the pore water pressure in the bank slope was greater. During the water-falling period, the scour depth and lateral erosion distances at the bank toe remained stable, whereas the stability of the bank slope decreased significantly with the decrease in the river water level, phreatic water level and pore water pressure, and the possibility of bank slope instability increased. This result confirms Rinaldi's conclusion that bank slope instability mainly occurs during the water-falling period of floods, and at the early stages of the subsequent flood peak and water-rising period (Rinaldi et al., 2008). When the scour erosion of the bank toe was not considered, K s2 increased and decreased with the river water level during the flood season, owing to seepage in the slope. The ratio of K s2 to K s1 is expressed by β, which reflects the degree of influence of scour erosion on bank slope stability. During the water-rising period, with the increase in the scour depth and lateral erosion distance at the bank toe, β increased gradually from 1.08 to 1.85. During the flood peak and water falling periods, the scour depth and lateral erosion distance at the bank toe remained stable, and β fluctuates in the range of 1.73-1.89. These results show that the slip stability of the bank slope was more sensitive to the change in the bank toe erosion than to the effect of the river water level or bank pore water pressure change, particularly when the water level in the river channel increased.
With the exceptions of the dry and flood seasons in the summer and autumn, flood disasters occurred frequently in the NRYR during the ice flood seasons in the winter, and the damage to the vulnerable zone during the ice flood season was very serious. The method outlined herein does not consider the formation of ice jams and dams, or the possibility of ice flood hazards, nor has it considered the bank slope stability characteristics of the vulnerable zone under different ice hydrodynamic conditions. In the future, a bank slope stability calculation model considering scour on the bank slope and the change in water level during the ice flood season will be established, and the mechanism of damage to the vulnerable zone under the ice flood effect will be discussed. In addition, the formation and development of damage in the ice flood season under sandy environments will be investigated, and the evolution mechanism of bank collapse will be revealed.

| CONCLUSIONS
A coupled model integrating a hydrosediment dynamics module with a vulnerable zone stability analysis module was established in this study, with its performance assessed in the NRYR. The local scour depth of the bank toe in the vulnerable zone exhibited a logarithmic variation, which increased gradually in the water-rising period and fluctuated in the flood peak and water-falling periods. The bank toe deformation in the vulnerable zone during the flood season was closely related to the water and sediment inflow conditions. A significant increase in the scour intensity index α resulted in bank toe scours in the vulnerable zone; conversely, siltation occurred.
The lateral erosion of the spur dyke increased gradually during the flood season because of the higher velocity and shear stress of the near-shore water flow. The probability of instability failure on the bank slope in the vulnerable zone was greater during the flood waterfalling period, with a decrease in phreatic water level and pore water pressure. Compared with the effect of river water level or bank pore water pressure in the flood season, the slip stability of the bank slope was more sensitive to the erosion of the bank toe, which was the main factor affecting the bank slope stability of the vulnerable zone.
The limitations of the proposed model are as follows: (1) the channel model was insufficient for characterizing F I G U R E 1 3 Change in near-shore water level and bank slope stability safety factor K s of #15 spur dyke at Caijiahekou under designed flood condition the complex ice flood progress in an alluvial river; and (2) bank slope stability characteristics under different ice hydrodynamic conditions were not considered. In the future, a coupled representation of sediment transport and ice flood routing in addition to improved modules for bank slope seepage-stability analysis under freezing conditions can be adopted.