Coupled two-dimensional modeling of bed evolution and bank erosion in the Upper JingJiang Reach of Middle Yangtze River

Geomorphology 344 (2019) 10–24

Contents lists available at ScienceDirect

Geomorphology journal homepage: www.elsevier.com/locate/geomorph

Coupled two-dimensional modeling of bed evolution and bank erosion in the Upper JingJiang Reach of Middle Yangtze River Shanshan Deng, Junqiang Xia ⁎, Meirong Zhou State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China

a r t i c l e

i n f o

Article history: Received 19 March 2019 Received in revised form 17 July 2019 Accepted 17 July 2019 Available online 19 July 2019 Keywords: Coupled modeling Bank erosion Bed deformation Groundwater flow Upper Jingjiang Reach

a b s t r a c t At the initial stage of dam operation in an alluvial river, significant vertical bed incision and transverse bank erosion processes always occur in the river reach downstream of the dam, with the latter being usually an integrated product of several interacting processes. A coupled two-dimensional (2D) model for the processes of bed evolution and bank erosion has been proposed in this study, which integrates a 2D morphodynamic module with a bank-erosion module and a groundwater flow module. In addition to the action of fluvial erosion, the effect of groundwater flow on bank erosion was considered, with the changes in pore water pressure and matrix suction being calculated. The proposed model was calibrated and validated through simulating channel evolution processes in a 30.54 km sub-reach of the Upper Jingjiang Reach of the Middle Yangtze River, with the effects of infiltration recharge on bank erosion and incoming flow on the central bar evolution being investigated. The results show that: (i) the proposed model had a good performance in reproducing flow and sediment factors, in terms of water depth, depth-averaged velocity and suspended sediment concentration, daily averaged discharge and river stage, and flow and sediment diversion ratios in local anabranching subreaches; (ii) the proposed model correctly calculated the bank erosion processes in the study reach, with the calculated bank retreat width being close to the measured value, however, it was not enough to fully reproduce the evolution of central bars, indicating that a specialized and improved module may be necessary for central bar evolution and its interaction with bank erosion; (iii) the groundwater level change lagged behind the river stage variation, and it showed a more obvious phase lag and a higher value at the same river stage (corresponding to a higher pore water pressure), when encountering with a more rapidly changing river stage; and (iv) the rate of infiltration recharge caused by rainfall process, would generally increase the groundwater level and thus the bank erosion degree in the study reach. In addition, the shifting of incoming main flow might take an obvious impact on the calculation of nearby central bar evolution. © 2019 Elsevier B.V. All rights reserved.

1. Introduction During the initial stage of dam operation, channel evolution in an alluvial river below the dam is usually characterized by vertical bed incision and transverse bank erosion processes, with the latter being an integrated product of several interacting processes such as soil hydrology, groundwater flow, fluvial erosion, and mass failure. Multidimensional numerical models have been proposed to predict river morphodynamic processes in recent decades (Darby et al., 2002; Cantelli et al., 2007; Jia et al., 2010; Abderrezzak et al., 2016; Langendoen et al., 2016; Morianou et al., 2018). In particular, coupled two-dimensional (2D) models are developed based on the shallow water and sediment transport equations, the bed-deformation equation, and the process-based bank stability analysis (Darby et al., 2002; Abderrezzak et al., 2016; Langendoen et al., 2016; Morianou et al., ⁎ Corresponding author. E-mail address: [email protected] (J. Xia).

https://doi.org/10.1016/j.geomorph.2019.07.010 0169-555X/© 2019 Elsevier B.V. All rights reserved.

2018). These models are widely used to simulate channel evolution processes in experimental and natural cases. Unfortunately, the roles of some key factors influencing bank erosion remain unclear, which are also case-dependent (Xia et al., 2014). The complexity of bank erosion mechanisms leads to a big difficulty in accurate modeling of bank erosion process (Rinaldi and Nardi, 2013). Many efforts have been made in exploring the mechanisms of bank erosion, with some important progress being achieved (Rinaldi and Nardi, 2013). Fluvial erosion is always regarded as the most crucial factor controlling the occurrence of bank erosion (Nardi et al., 2013; Xia et al., 2014). Its quantification has been improved through the development of specific techniques for in situ measurements and through nearbank flow modeling (Rinaldi and Nardi, 2013). The influences of riparian vegetation and subaerial processes on bank erosion gradually have become an important part of recent studies (Simon et al., 2006; Camporeale et al., 2013; Rinaldi and Nardi, 2013; Kimiaghalam et al., 2015). The interaction between bank stability and groundwater flow was also investigated by some researchers (Darby et al., 2007; Rinaldi

S. Deng et al. / Geomorphology 344 (2019) 10–24

et al., 2008; Lindow et al., 2009; Fox and Wilson, 2010; Karmaker and Dutta, 2013; Rinaldi and Nardi, 2013; Masoodi et al., 2018). The changes in pore water content and pressure, associated with groundwater flow, are now recognized as one of the most important factors controlling the onset and timing of bank failure (Simon et al., 2000; Darby et al., 2007; Rinaldi et al., 2008; Rinaldi and Nardi, 2013). Rinaldi and Nardi (2013) pointed out that a crucial point is the transient characteristics driven by dynamic hydrological variables when accounting for the action of pore water pressure on bank erosion. Quantification of its temporal change requires an integrated analysis of various hydrological processes involved (such as rainfall, river hydrographs, and groundwater flow). Bank erosion models were improved by combining groundwater flow calculations with fluvial erosion calculations and bank stability analysis to illustrate the interactions between these processes (Darby et al., 2007). These models were then integrated with a 2D hydrodynamic model to simulate bank erosion process in alluvial rivers (Rinaldi et al., 2008; Luppi et al., 2009). However, the manual procedure for mesh adaptation and the ignorance of bed evolution may limit the application of these models in simulating the whole process of channel evolution. Some previous models implemented the function of simultaneous calculation of both bed evolution and bank erosion processes (Jia et al., 2010; Xia et al., 2013; Abderrezzak et al., 2016; Langendoen et al., 2016). However, the mechanisms of bank erosion were simplified in order to reduce the model complexity in these studies, with a few incorporating the influences of some non-fluvial factors such as groundwater flow. The objectives of the current study are to: (i) develop a coupled 2D model for the channel evolution process, through integrating a 2D morphodynamic module with a groundwater flow module and a bank erosion module (covering fluvial erosion calculation and bank stability analysis); (ii) apply the proposed model to simulate channel evolution in the Upper Jingjiang Reach in order to assess model performance; and (iii)

quantitatively investigate the influences of some factors on channel evolution. 2. Study area The Jingjiang Reach is located in the Middle Yangtze River, 102 km downstream of the Three Gorges Dam (TGD), as shown in Fig. 1a. It is conventionally divided into the Upper and Lower Jingjiang Reach (UJR and LJR), owing to the difference in geographical characteristics. The study area is an anabranching and slightly curved sub-reach (30.54 km) of the UJR, which comprises three central bars (Fig. 1b). Riverbanks in the study reach have a two-layer soil structure, with a thick cohesive upper layer (8–16 m) and a usually submerged thin noncohesive lower layer (Xia et al., 2016). Bed material in the study reach mainly comprises medium and fine sand, with a median diameter of about 0.20 mm (Xia et al., 2016). Under the impact of dramatically decreased incoming sediment load after the operation of the Three Gorges Project (TGP) in 2003, the study reach has been subjected to significant channel adjustments, covering obvious bed incision, bank erosion and central bar evolution. The maximum bank erosion rate in the study reach reached 37 m/a in 2003–2016, occurring in the sub-reach near Lalingzhou. The three central bars composed of sand and gravel in the study reach were significantly eroded, with the total area decreasing by around 70% in 2003–2016. Continuous erosion occurred at the central bars of Sanbatan and Jingchengzhou (Fig. 1b), leading to a 91.2% and 70.1% reduction in their areas, respectively. The evolution process of the central bar of Taipingkou (Fig. 1b), close to Lalinzhou, was different than the other two bars. Its area was firstly increased to almost 2 times the initial area during 2003–2010 and then decreased by 62.6% during 2011–2016. The transfer from aggradation to erosion at the Taipingkou central bar was speculated to be related to the bank erosion process at Lainlinzhou (Li et al., 2018), which coincidently ceased after 2011 because of the construction of bank revetment works. In addition,

Upper Yangtze River

Middle Yangtze River

(4500km)

Qinghai-Tibet Plateau

11

(955km)

Lower Yangtze River (938km)

Three Gorges Dam Yichang

(a)

Shanghai

Yangtze River Basin Zhicheng

Hukou Chenglingji

Jingjiang Reach Dongting Lake

(b) 2 1 4

3

1 Taipingkou 4 Lalingzhou 2 Sanbatan 3 Jinchengzhou

Fig. 1. Sketched maps of the Yangtze River and Jingjiang Reach.

12

S. Deng et al. / Geomorphology 344 (2019) 10–24

bedload transport in the study reach accounted for b4% of the suspended sediment transport (Xia et al., 2016), and thus was assumed negligible in the current study.

The flow moment equations in the ξ and η directions can be respectively given by ∂U U ∂U V ∂U UV ∂C ξ V 2 ∂C η − ¼ þ þ þ ∂t C ξ ∂ξ C η ∂η C ξ Cη ∂η C ξ C η ∂ξ

3. Model framework In order to simulate the full process of channel evolution, covering flow and sediment transport, bed deformation and bank erosion, a coupled model has been proposed herein. This model integrates a 2D morphodynamic module with two modules respectively for bank erosion and groundwater level change. Fig. 2 shows the flow chart of the whole simulation process. As indicated in Fig. 2, the 2D morphodynamic module is first used to calculate the detailed flow and sediment parameters and bed deformation. With the obtained flow conditions, the groundwater flow module is then used to calculate the variation of groundwater level inside the bank, which is associated with the changes in pore water pressure and matrix suction. The bank erosion module is parameterized by the outputs of the former two modules, and is adopted to determine the occurrences of bank-toe erosion and mass failure. 3.1. Morphodynamic module The depth-averaged 2D governing equations for shallow water in curvilinear coordinates are adopted to calculate the flow conditions (Xia et al., 2013). The flow continuity equation can be written as:   ∂Z 1 ∂
1 ∂
þ UhC η þ VhC ξ ¼ 0 C ξ C η ∂η ∂t C ξ C η ∂ξ

ð1Þ

Start

and ∂V U ∂V V ∂V UV ∂C η U 2 ∂C ξ − ¼ þ þ þ ∂t C ξ ∂ξ C η ∂η C ξ Cη ∂ξ C ξ C η ∂η g ∂Z ‐gn2 ‐ C η ∂η

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 U2 þ V 2 νt ∂ V ∂ V V þ ð þ Þ 4=3 C ξ C η ∂ξ2 ∂η2 h

Change in bank profile

1D groudwater flow

Mass failure (including fluvial erosion caclulation and stability analysis)

Fs
ð4Þ

where εξ and εη are the diffusion coefficients in the ξ and η directions; Sk is the suspended sediment concentration for the kth fraction (kg/m3); Sbk is the sediment source for the kth fraction, contributed by bank erosion (kg/(m2∙s)); Ek is the entrainment rate of bed material for the kth fraction (kg/(m2∙s)); Dk is the deposition rate of suspended sediment for the kth fraction (kg/(m2∙s)). The source term (Ek-Dk) is calculated by Ek − Dk = ωkαk(S∗k ‐ Sk), where ωk is the sediment settling velocity for the kth fraction (m/s); and αk is the corresponding recovery coefficient, representing the channel self-adjustment capacity to reach an equilibrium state. The suspended sediment transport capacity S∗k (kg/ m3) is calculated using the formula proposed by Zhang (1998), which is written as: 3

Planar failure occurs

ð3Þ

where t is the time (s); Z is the water level (m); h is the water depth (m); U and V are respectively velocity components in the ξ and η directions (m/s); Cξ and Cη are Lamé coefficients; n is the Manning roughness coefficient; g is the gravitational acceleration (9.81 m2/s); and νt is the horizontal turbulence viscosity coefficient (m2/s), and νt = φu∗h, where φ is an empirical parameter ranging from 0.0 to 1.0; u∗ is the fricpffiffiffi −1=6 tion velocity (m/s), calculated by u ¼ n g Uh ; U is the total velocity pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 (m/s), and U ¼ U þ V . Alternative direction implicit method (ADI) is adopted to solve the above governing equations on a fixed and orthogonal curvilinear mesh. Without considering the bedload transport, the 2D governing equation for suspended sediment transport is given by (Xia et al., 2013):

Cξ ∂ ∂ þ ½hεη ðS Þg þ ðEk −Dk Þ þ Sbk C η ∂η k ∂η

2D morphodynamic module Lateral sediment source

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 U2 þ V 2 νt ∂ U ∂ U Uþ ð 2 þ 2Þ 4=3 C ξ C η ∂ξ ∂η h

 ∂
 Cη ∂ ∂ 1 ∂
1 ∂ ðhSk Þ þ ½ C η UhSk þ f ½hε ðS Þ C VhSk  ¼ C ξ C η ∂ξ C ξ C η ∂ξ ξ C ξ ∂ξ k ∂t ∂η ξ

Initial conditions (Nt=0)

Yes

g ∂Z ‐gn2 ‐ C ξ ∂ξ

ð2Þ

Sk ¼ ΔP k λ

U ghωs

!m ð5Þ

Change in bank profile

Nt ≤ N

Yes

No

End Fig. 2. Flow chart of the simulation process.

where ΔP⁎k is the percentage of sediment transport capacity for the kth fraction; ωs is the mean settling velocity of non-uniform suspended sediment (m/s); λ and m are empirical coefficients (Zhang, 1998). ΔP⁎k is a weighted average of the percentage of suspended sediment ΔPsk and the percentage ΔPlk estimated by the formula of Li (1987), which can be written as:
 ΔP k ¼ θp ΔP sk þ 1−θp ΔP lk

ð6Þ

S. Deng et al. / Geomorphology 344 (2019) 10–24

and −6ωk 1−Ak ð1−e κu Þ ωk ΔP lk ¼ ΔP bek −6ωk M X 1−Ak ΔP bek ð1−e κu Þ ωk k¼1

ð7Þ

where θp is the weight coefficient, ranging between 0.0 and 1.0; M is the number of fractions; ΔPbek is the percentage of bed material for the kth fraction; κ is Kármán constant; Ak is a coefficient that relates to the vertical turbulence intensity δv and the normal distribution function Φ (ωk/ δv). More details in calculating ΔPlk can be referred to Li (1987). ΔPsk is related to both the incoming sediment condition, and the suspended sediment transport process described by Eq. (4). Bed evolution would lead to a change in ΔPbek, and consequently cause a variation in ΔP⁎k. The detailed calculation procedure of ΔPbek can be referred to the study of Wang et al. (2008). The deformation rate of the river bed, caused by the non-equilibrium transport of suspended sediment, can be given by (Xia et al., 2013): ρ0

M M X ∂Z b X ðDk ‐Ek Þ ωk α k ðSk −Sk Þ ¼ ¼ ∂t k¼1 k¼1

ð8Þ

3.2. Module of groundwater table variation The groundwater table inside riverbank usually shows a seasonal variation, influenced by dynamic hydrological variables, such as the varying river stage and the rainfall and evaporation (Darby et al., 2007). The variation in the groundwater table is greatly associated with the changes in pore water pressure and soil properties, which may play an important role in bank erosion processes in some alluvial rivers. In the study reach, the riverbanks consist of a thick cohesive soil layer of low permeability, which would lead to a slow decrease of the groundwater table (or pore water pressure) in the recession period and may greatly influence the bank stability. Note that groundwater flow may also induce seepage erosion to bank soil particles, whereas it is seldom reported in the study reach and thus is not considered herein. The one-dimensional (1D) governing equation for groundwater flow with a free surface can be written as (Chiang et al., 2011): μ

 ∂Z g ∂Z g ∂ ¼ kc Þþq Zg ∂t ∂y ∂y

ð9Þ

where Zg is the groundwater level (m); kc is the hydraulic conductivity pffiffiffiffiffi (m/s); μ is the specific yield, estimated by μ ¼ 0:117 7 kc (Gu, 2000); q is the recharge rate from infiltration (m/s); and y is the spatial coordinate (m).

3.3. Bank erosion module Bank erosion in alluvial rivers usually consists of two processes: fluvial erosion by flowing water and mass failure by gravity and pore water pressure (Thorne et al., 1998; Darby et al., 2007; Langendoen et al., 2016). The calculations of these processes are presented in the following sections, with relevant parameterization also illustrated. 3.3.1. Fluvial erosion The fluvial erosion rate relies on both flow intensity and the resistance-to-erosion capacity of bank soil, and is usually expressed by an empirical relation based on the excess flow shear stress, which can be written as (Darby et al., 2007): Ef ¼


a kd τ f ‐τc 0

3.3.2. Mass failure Based on field observations and laboratory experiments, planar and rotational failures are the most common modes for mass failure in the UJR (Xia and Zong, 2015). However, only planar failure was considered in the current study, with the bank stability analysis approach proposed by Osman and Thorne (1988) being improved to determine the occurrence of mass failure. The influences of pore water pressure and matrix suction variations, resulting from the varying groundwater table, are considered in this study, along with the effect of the hydrostatic confining force of in-channel water.

i0

i

1

j toe-1

G

River stage

ð10Þ

where Ef is the fluvial erosion rate (m/s); τf is the flow shear stress (N/ m2); τc is the critical shear stress of bank soil (N/m2); kd is the erodibility coefficient (m3/(N∙s)); and a is an empirically-derived exponent, generally assumed to be 1.0. τf is obtained from the 2D morphodynamic module, with the value at the computational node near the bank-toe being adopted. According to the study of Xia and Zong (2015), τc ranges from 0.17 N/m2 to 1.18 N/m2 for cohesive bank soils in the Jingjiang reach, and increases with the dry density, whereas it varies between 0.14 N/m2 and 0.52 N/m2 for non-cohesive bank soils. Based on the previous studies (Hanson and Simon, 2001; Karmaker and Dutta, 2011), kd can be empir0 ically calculated bykd ¼ λ0 τc m , where λ′ was calibrated in the current research to be around 2.0 times the value in the research of Hanson and Simon (2001), and m′ directly used the value in their research.

Groundwater table

H2

if τ f Nτc if τ f ≤τc

j toe-2

Ht

H1

Both pore water pressure and matrix suction vary with the groundwater table, with the former acting as a destabilizing force and the latter contributing to the apparent cohesion of bank soil. A simple hypothesis is made in this study that the pore water pressure and matrix suction hydrostatically vary with the distance below and above the groundwater table (Rinaldi and Casagli, 1999). In addition, the average groundwater level, within a 5 m distance from the bank surface, is adopted to calculate the pore water pressure and matrix suction, which are then used to parameterize the bank erosion module.



where Zb is the bed level (m); and ρ’ is the dry density of bed material (kg/m3).

13

PU

P

j toe

2

sbk

3

j top Fig. 3. Sketched bank profile at the mode of planar failure.

Fig. 4. Scheme for lateral sediment source distribution.

14

S. Deng et al. / Geomorphology 344 (2019) 10–24

(i, j te-1) (i, j te)

(i, j tp)

(i, jte ) New

Computational node

(i, jtp +1)

(i, jtp ) New

(i, j te-2) (i, j te-1)

Bank-top

(i, j te)

(i, jtp )

Bank-top

Computational node

(b)

(a)

Bed incision

Bank-toe

Bed incision

Bank-toe

Fig. 5. Bank profile representation and adjustment.

Mass failure can be divided into initial failure and subsequent parallel failures, based on the approach of Osman and Thorne (1988). The sketched bank profile is shown in Fig. 3. The failure surface is assumed to pass through the bank-toe, with the presence of tensional crack on the bank-top. With these assumptions, the failure surface angle and bank stability degree are required to be calculated, in order to determine the occurrence of mass failure and the bank geometry after failure. Bank stability degree: The bank stability degree is normally represented by the safety factor Fs, namely the ratio of the maximum resistant force to the driving force exerted on the potential failure surface. Including the contributions of pore water pressure and matrix suction, Fs can be given by

Thorne, 1988), and can be calculated by:

β¼

1 2

( tan−1

"

H1 H2

# ) 2  1:0−K 2 tani0 þ ϕ

ð12Þ

where K is the ratio of Ht to H1, and Ht = (2c/γ) tan (45° + ϕ/2) (Wang and Kuang, 2007); γ is the natural specific weight of bank soil (kN/m3); c is the cohesion of bank soil (kN/m2); ϕ is the friction angle of bank soil (°); i0 is the initial bank angle (°); and H2 is the bank height (m) shown in Fig. 3. The surface angle for parallel failure equals the bank angle (β) for the initial failure. 3.4. Lateral sediment source


 0 0 c L þ Sm tanϕb þ N p −U l tanϕ Fs ¼ G sinβ þ P V cosβ−P sinθ

ð11Þ

where G is the unit-length weight of the potential failure block (kN/m), with the natural and saturated specific weight of bank soil being used for the portions above and under the groundwater table, respectively; PV is the pore water pressure on the tension crack (kN/m); P is the hydrostatic confining force (kN/m); θ is the angle between the direction of P and the normal direction of failure surface (°); c′ is the effective cohesion (kN/m2); ϕ′ is the effective friction angle (°); L is the length of failure surface (m), with L = (H1-Ht)/sinβ; H1 is the total bank height (m); Ht is the tension crack depth (m); β is the bank angle (°); tanϕb represents the rate of increase in shear strength relative to the matric suction; Sm is the total matric suction (kN/m) on the unsaturated portion of failure surface; Np is the total pressure on the failure surface (kN/m), and Np = Gcosβ + Pcosθ; Ul is the total uplift force (kN/m), with Ul = PU + PV sinβ, where PU is the pore water pressure on the saturated portion of failure surface (kN/m). Failure surface angle: The failure surface angle β for the initial failure depends on bank soil properties and bank geometry (Osman and

According to the research of Xia and Zong (2015), almost 50% of the failed bank material would be flushed away immediately by the flowing water after the occurrence of a mass failure event, with the remaining part depositing at the bank-toe region. Therefore, the bank erosion volume VsE, i that is transformed into suspended sediment, is f m f calculated by VsE, i = VE, i + RtsVE, i, where VE, i is the fluvial erosion volume (m3/m), Rts is the transform proportion of failed material into suspended sediment concentration and Rts = 50%; and th 3 m VE, i is the mass failure volume at the i section (m /m). Furthermore, two assumptions were made to calculate the lateral sediment source Sbk in Eq. (4). First, the length of eroded bank was assumed to be the product of Reb and the total length (Δξi) within the sub-reach between the (i-1)th and the ith sections. Second, Sbk was assumed to be uniformly distributed between the bank-toe nodes (i, jtoe) and (i, jtoe-1), and linearly distributed between the nodes (i, jtoe-2) and (i, jtoe-1) and between the nodes (i, jtoe) and (i, jtop), with a value of zero at nodes (i, jtoe-2) and

Jing34 Sha4

Soil properties

Jing34

Jing45

Jing55

Natural specific weight γ (kN/m3) Saturated specific weight γs (kN/m3) Hydraulic conductivity kc (cm/s) Measured cohesion c (kN∙m−2) Measured friction angle φ (°) *Effective cohesion c΄ (kN∙m−2) *Effective friction angle φ΄ (°) Dry density (t/m3)

18.46 18.85 8.95 × 10−5 8.80 27.90 6.20 25.10 1.30

17.89 18.09 5.42 × 10−7 22.45 17.45 15.72 15.71 1.34

18.19 19.16 2.08 × 10−6 15.30 27.80 10.71 25.02 1.52

* c΄ equals 70% of c, and φ΄ accounts for 90% of φ (Xia and Zong, 2015) .

Jing42

Jing45

Lalinzhou

Wan15

Table 1 Bank soil properties measured at three sections in or close to the study reach.

Jing36

Jing32

Jing33 Jing35 Jing29

Fixed section

Jing51 Jing52

Fig. 6. Computational mesh adopted in the model calibration.

S. Deng et al. / Geomorphology 344 (2019) 10–24

1.5 1.0

10.0

0.0 20.0 15.0

Depth (m)

1.5 1.0

10.0 5.0

800

1200

0

1600 20.0

(g)

2.0

Depth (m)

1.5 1.0

1200

400

800

5.0 800

0.0 0

600

1200

600

1200

Depth (m)

(m) 1.5 1.0

400

800

5.0

25.0

Depth (m)

1.5 1.0 0.5

400

800

400

800

15.0 10.0

1200

Distance from the left bank (m)

1200

400

800

1800

(o)

0.10

0 0.30

(q)

5.0 0

600

0.20

1200

20.0

0.0

0.10

0.00 0

(p)

1200

0.20

(n)

10.0

1200

800

(l)

0.30

0

0.0 0

400

0.30

15.0

0.5

0.10

1800

20.0

2.0

0.20

0.00 0

1800

1600

(i)

0

SSC (kg/m3)

(k)

5.0

0.5

1200

0.30

1200

10.0

1.0

800

0.40

15.0

Depth (m)

Velocity (m/s)

(j)

400

0.00 400

20.0

1.5

0.20

0

(h)

0

2.0

(f) 0.40

0.40

10.0

1200

0.10

1600

0.0 0

Velocity (m/s)

800

15.0

0.5

Velocity (m/s)

400

SSC (kg/m3)

Velocity (m/s)

2.5

400

SSC (kg/m3)

0

0.20

0.00

0.0

0.5

2.0

(e) SSC (kg/m3)

(d)

2.0

0.30

0.00 0.60

SSC (kg/m3)

Velocity (m/s)

15.0

5.0

0.5 2.5

(c)

0.40

20.0

Depth (m)

Velocity (m/s)

2.0

(b)

25.0

Cal. Meas.

SSC (kg/m3)

(a)

2.5

15

1200

(r)

0.20 0.10 0.00

0

400

800

1200

Distance from the left bank (m)

0

400

800

1200

Distance from the left bank (m)

Fig. 7. Comparisons between the calculated and measured depth-averaged velocities, depths and SSCs during the 2004 flood period at six fixed sections: (a)–(c) at Wan15; (d)–(f) at Sha4; (g)–(i) at Jing35; (j)–(l) at Jing38; (m)–(o) at Jing45; (p)–(r) at Jing51.

(i, jtop), as shown in Fig. 4. Sbk can be thus calculated by: Sbk ¼

V sE;i−1;i ρb 2
 ΔP bk Δt Δξi Δη1 þ 2Δη2 þ Δη3

ð13Þ

where Δη1, Δη2 and Δη3 are near-bank transverse intervals (m), as shown in Fig. 4; Δt is the time step (s); VsE, i−1, i is the bank erosion volume (m3), which is transformed into suspended sediment within Δξi, and calculated by VsE, i−1, i = 0.5Reb(VsE, i−1 + VsE, i)Δξi; ρb is the dry density of bank soil (kg/m 3); and ΔP bk is the percentage of bank soil for the kth fraction. In addition, the remaining part of failed bank soil is assumed to be evenly deposited within the bank-toe region, without considering its influence on near-bank bed material

gradation. The corresponding increase in bed level (ΔZb′) is evalum ated by ΔZ b′ = 0.5Reb (1 ‐ Rts )(Vm E, i + V E, i−1 )/L d , where L d is the width of bank-toe region (m), and is set to 5% - 10% of the channel width in this study. 3.5. Mesh adaption A fixed mesh is adopted in this study, with the bank-top (i, jtp) and bank-toe nodes (i, jte) being recorded at each section. If the bank-top passes through the node (i, jtp) after a bank failure (Fig. 5a), the nodes (i, jtp) and (i, jtp + 1) would be respectively recorded as the new bank-toe node (i, jte)New and bank-top node (i, jtp)New, with the level of node (i, jtp) being set to the level of node (i, jte). These bed level changes

S. Deng et al. / Geomorphology 344 (2019) 10–24

(a)

40

42 Cal. Meas.

30 20 10 30/Jun. 30/Jul. 29/Aug. 28/Sep. 28/Oct. Date

(b)

2.0 SSC (kg/m3)

50

River stage (m)

Discharge (103 m3/s)

16

40 38 36

(c)

1.6 1.2 0.8

34

0.4

32 30/Jun. 30/Jul. 29/Aug. 28/Sep. 28/Oct. Date

0.0 30/Jun. 30/Jul. 29/Aug. 28/Sep. 28/Oct. Date

Fig. 8. Comparisons between the calculated and measured hydrographs in the 2004 flood season at Shashi: (a) daily averaged discharge; (b) river stage; and (c) SSC.

would provide a feedback on the subsequent hydrodynamic simulation. However, if the bank-top does not pass through the adjacent node (i, jtp) as shown in Fig. 5b, the bank profile change would not have any influence on the subsequent hydrodynamic simulation. In this case, the feedback between bank erosion and near-bank flow transport is reasonably ignored, given the fact that the magnitude of bank retreat width is small when scaled to the channel width in the study reach such that bank erosion would not cause a significant change in the depth-averaged flow conditions over a relatively short period. 4. Model application Model calibration and validation processes were conducted herein, by simulating the channel evolution processes in, respectively, the 2004 and 2008 flood periods and comparing the calculations with the corresponding measurements. 4.1. Model calibration 4.1.1. Model setup The proposed model was adopted to simulate the channel evolution process during the 2004 flood period (from July to October) in the study reach. The topography measured in July 2004 was collected from Changjiang Water Resources Commission (CWRC), along with the bed material gradation measured in 2003 and the daily averaged hydrological data at the Shashi station. The depth-averaged flow velocities and suspended sediment concentrations measured at six fixed sections were also collected from the CWRC to assess the model performance. In addition, cross-sectional profiles at 28 fixed sections in the study reach were annually measured in October by CWRC. Bank soil samples at three sections (Fig. 1b), in or close to the study reach, were collected by Xia and Zong (2015), with indoor geotechnical tests conducted to

Jing33: Wm=13.5 m Wc=10.7 m Jing34: Wm=11.9 m Wc=9.0 m Jing35: Wm=3.1 m Wc=5.1 m

measure the physical and mechanical soil properties, as shown in Table 1. Bank soil properties adopted in the proposed model were interpolated from the measured soil properties at these three sections. Due to the lack of hydrometric stations at the inlet and outlet of the study reach, the flow and sediment hydrographs obtained from a 1D morphodynamic model was used to provide the boundary conditions of the 2D model in this study. The 1D model was first adopted to simulate the flow and sediment transport processes in the sub-reach between the hydrometric station of Zhicheng and the water gauge station of Haoxue, which was chosen to cover the study reach, as shown in Fig. 1. The measured daily averaged water discharge and sediment concentration at Zhicheng were used as the upstream boundary conditions for the 1D model, with the measured river stage at Haoxue as the downstream boundary condition. The calculated results of the 1D model were compared with the corresponding measurements at the hydrometric stations of Zhicheng and Shashi and the water gauge station of Chengjiawan, with close agreement obtained. The rootmean-squared-errors (RMSEs) for discharge and sediment concentration were one to two orders of magnitude less than the corresponding mean values, with the RMSE for river stage being b2.5% of the variation range. Therefore, the calculated results of the 1D model were believed to be able to correctly characterize the boundary conditions for the 2D model in the current study. Based on the calculation of the 1D model, the discharge at the inlet section (Jing29) of 2D modeling varied from 10,945 m3/s to 49,613 m3/s during the 2004 flood period, carrying a sediment concentration of 0.091–1.778 kg/m3, and the river stage changed from 34.05 m to 41.65 m. The computational mesh used in the 2D simulation was locally refined and consisted of 205 × 40 nodes, as shown in Fig. 6. The minimum mesh size was 51.1 m in the streamwise direction and 14.5 m in the transverse direction. The locally refined area was chosen to cover the right bank-toe region at Lalinzhou where serious bank erosion usually occurred.

Zb 5m 10m 15m 20m 25m 30m 35m 40m

2m/s

IBT

Lalinzhou FBT

Fig. 9. Comparison between the initial and final bank-top lines in the 2004 flood period.

S. Deng et al. / Geomorphology 344 (2019) 10–24

17

Jing34

Jing42 Jing33

Jing35

Sanbatan

Final central bar Initial central bar

Zb

50 40 30 20 10 0m

Jinchengzhou

Fig. 10. Calculated evolution processes of central bars over the 2004 flood period, with the boundary of central bar characterized by the 30 m contour.

4.1.2. Comparisons of flow and sediment conditions Fig. 7 shows the comparisons between the calculated results using the proposed 2D model and the measurements at six fixed sections, in terms of depth-averaged velocity, water depth and suspended sediment concentration (SSC). As indicated in Fig. 7, the calculated results generally agreed well with the corresponding measurements. The calculated cross-section averaged velocities ranged from 1.15 to 1.44 m/s at these sections, close to the measured values of 1.29–1.55 m/s, with the absolute errors b0.40 m/s. The difference between the calculated and measured cross-section averaged depths ranged between 0.14 m and 0.62 m at sections of Jing35 - Jing51, with a relative error (RE) of b5%; whereas it was 1.60–1.82 m at sections Wan15 and Sha4, with an RE of 12–20%. The difference between the calculated and measured values of cross-section averaged SSC had a maximum value of 0.08 kg/m3 at Jing35 and a minimum value of 0.01 kg/m3 at Sha4. In addition, the calculated transverse distributions of depth-averaged velocity and water depth also agreed relatively well with the measurements. For example, the calculated maximum velocity Vmax (1.69 m/s) at Jing38 occurred at a distance of 359 m from the left bank, comparing well with the measured Vmax of 1.65 m/s at a distance of 224 m. The calculated maximum depth of 17.54 m at Jing38 occurred at a distance of 881 m from the left bank, comparing well with the measured value of 15.30 m at a distance of 973 m. The calculated transverse distribution of depth-averaged SSC was relatively less satisfying, but still reflected the general characteristics, especially at sections of Jing38 and Jing51, as shown in Fig. 7. Fig. 8 shows the comparisons between the calculated and measured daily mean hydrographs of water discharge, river stage and SSC at Shashi. As indicated in Fig. 8, the calculations also agreed well with the measurements. The RMSEs for discharge and SSC were 1507 m3/s and

(a)

Shashi sation

Taipingkou

Sanbatan Songzikou Reference station˖ Shashi River stage: 29.57m Time: 2004/3/15 Coordinate system˖ WGS84 49N

Jinchengzhou

0.08 kg/m3, respectively, which accounted for 7.81% and 30.40% of the corresponding average values. The RMSE for river stage was 0.33 m, accounting for 4.09% of the variation range (8.01 m). Additionally, the calculated maximum discharge was just 3.85% less than the measured value, with the relative error of the maximum SSC being 20.43%. The difference between the calculated and measured maximum river stages was 0.21 m, accounting for 2.57% of the variation range. Therefore, it can be concluded that the proposed model reproduced well the flow and sediment transport processes in the study reach during the 2004 flood period. 4.1.3. Bank erosion Fig. 9 shows the comparison between the initial bank-top line (IBT) and the final bank-top line (FBT) during the 2004 flood period, and the flow field (at the calculated maximum discharge of 45,186 m 3 /s at the Shashi station) near Lalinzhou where obvious bank erosion occurred. Considering that the amount of bank retreat (measured as a width) was small when scaled to the channel dimensions, it was magnified by a factor of 10 in Fig. 9 for a better view. As can be seen from Fig. 9, there were high velocities near the right bank at Lalinzhou, indicating that intensive fluvial erosion occurred in the right bank-toe region. Bank erosion was calculated to mainly occur along the right banks between sections of Jing33 - Jing35, which was consistent with the findings based on the annually surveyed profiles at the 28 fixed sections. In addition, the calculated retreat widths (Wc) at sections of Jing33 - Jing35 were respectively 10.7 m, 9.0 m and 5.1 m, which were close to the measured values (Wm) of 13.5 m, 11.9 m and 3.1 m (Fig. 9). It should be noted that the measured retreat width was obtained from the annually surveyed cross-sectional profiles, and thus it represented the bank erosion

(b)

Boundary of central bar in 2004

Reference station˖ Shashi River stage: 29.82 m Time: 2005/4/3 Coordinate system˖ WGS84 49N

Fig. 11. Remote sensing images taken on 15 Mar. 2004 and 3 Apr. 2005 when the river stages were close to 30 m at Shashi.

18

S. Deng et al. / Geomorphology 344 (2019) 10–24

(a)

Depth (m)

15.0

1.0

10.0

2.5

0 20.0

(d)

1.5

0.10

(e)

10.0

1.0

5.0

0.5

2.5

0

200 400 600 800 1000 1200 1400 1600

(g)

Depth (m)

0.5 0.0 200

400

600

5.0

800 1000 1200 1400 15.0

200

400

600

0.9

0 0.10

(k)

5.0

0.7 0.5

0.05

800 1000 1200 1400

10.0

1.1

0.0 0

1.5

300

600

900

1200

1500

1800

300

600

900

1200

1500

SSC (kg/m3)

Depth (m)

5.0

0

200

400

600

(p)

2.0

200

400

600

15.0 10.0

0.0 0

200

400

600

800 1000 1200 1400

Distance from the left bank (m)

900

1200

1500

1800

0.04

200

400

600

800 1000 1200 1400

(r)

0.08 0.06 0.04 0.02

5.0

0.0

600

0.06

0

SSC (kg/m3)

Depth (m)

0.5

300

(o)

0.08

0.10

20.0

1.0

0.02

800 1000 1200 1400

(q)

25.0

1.5

0.04

0.00 0

800 1000 1200 1400

0.06

0.02

0.0

0.5

800 1000 1200 1400

(l)

0

(n)

10.0

0.7

600

0.08

1800

15.0

0.9

400

0.10

20.0

1.1

200

0.00 0

(m)

1.3

(i)

0.10

SSC (kg/m3)

Depth (m)

1.3

200 400 600 800 1000 1200 1400 1600

0.00 0

(j)

1.5

Velocity (m/s)

0

0.0 0

Velocity (m/s)

0.02

0.15

10.0

1.0

0.04

200 400 600 800 1000 1200 1400 1600

15.0

1.5

0.06

(h)

20.0

2.0

(f)

0.08

SSC (kg/m3)

0

200 400 600 800 1000 1200 1400 1600

0.00

0.0

0.0

Velocity (m/s)

0

200 400 600 800 1000 1200 1400 1600

15.0

Depth (m)

Velocity (m/s)

2.0

0.05 0.00

5.0 200 400 600 800 1000 1200 1400 1600

0.10

SSC (kg/m3)

Velocity (m/s)

1.5

0

(c)

0.15

Meas.

0.5

Velocity (m/s)

(b)

20.0 Cal.

SSC (kg/m3)

2.0

0.00 0

200

400

600

800 1000 1200 1400

Distance from the left bank (m)

0

200

400

600

800 1000 1200 1400

Distance from the left bank (m)

Fig. 12. Comparisons between the calculated and measured depth-averaged velocities, depths and SSCs during the 2008 flood period at six fixed sections: (a)–(c) at Wan15; (d)–(f) at Sha4; (g)–(i) at Jing35; (j)–(l) at Jing38; (m)–(o) at Jing45; (p)–(r) at Jing51.

process during the whole hydrological year of 2004 (October 2003– October 2004), whereas the current simulation presented the calculated width during the 2004 flood period of July–October. However, based on the field investigations by CWRC, about 72% of bank failures occurred during the period of July–October in the Middle and Lower Yangtze River (CWRC, 2014). The proposed model was thus considered to be able to reproduce the major bank erosion process during the whole hydrological year.

4.1.4. Central bar evolution The boundary of the central bar in the study reach was characterized by the 30 m contour of bed level. Fig. 10 gives a comparison between the calculated initial and final boundaries, in order to demonstrate the evolution processes of central bars during the study period. Calculated results (Fig. 10) show that: (i) three central bars were obviously eroded, with a decrease in their areas; (ii) the central bar of Taipingkou shifted downstream due to a degradation at the upstream end and an

(a)

Cal. Meas.

30.0 20.0 10.0 0.0 30/Jun.

30/Jul.

29/Aug.

28/Sep.

40

0.80

(b)

(c)

0.60

36

0.40

32 28 30/Jun.

28/Oct.

19

SSC (kg/m3)

40.0

River stage (m)

Discharge (103 m3/s)

S. Deng et al. / Geomorphology 344 (2019) 10–24

0.20

30/Jul.

29/Aug.

28/Sep.

28/Oct.

0.00 30/Jun.

Date

Date

30/Jul.

29/Aug.

28/Sep.

28/Oct.

Date

Fig. 13. Comparisons between the calculated and measured hydrographs in the 2008 flood season at Shashi: (a) daily averaged discharge; (b) river stage; and (c) SSC.

aggradation at the downstream end; and (iii) obvious erosion occurred along the right margin of Sanbatan and at the upstream end of Jinchengzhou, with a slight aggradation at the downstream end of Sanbatan. Fig. 11 shows the remote sensing images of the study reach on 15 March 2004 and on 3 April 2005, when the river stages at Shashi station were close to 30 m. The comparison between these two images indicates that: (i) a dramatic erosion occurred along the right margin of Sanbatan, with a slight aggradation at the downstream end; (ii) the central bar of Jinchengzhou was significantly eroded, with the upstream end shifting downstream and the downstream end being cut off. In addition, the erosion along the right margin of Jinchengzhou seemed to promote the development of the right branch; and (iii) the upstream aggradation and downstream degradation ultimately led to an upstream shifting and an expansion of the central bar of Taipingkou. Comparing the calculated results with the observations from remote sensing images illustrates that the calculation reflected the major evolution characteristics of the bars of Sanbatan and Jinchengzhou, but incorrectly reproduced the evolution trend of Taipingkou in the study period. This disagreement may be partly caused by the insufficient characterization of the change in incoming flow regime and the ignorance of flow diversion of the tributary of Songzikou (Fig. 11).

4.2. Model validation 4.2.1. Model setup In order to carry out the validation for the 2D model, flow and sediment transport processes in the sub-reach between Zhicheng and Shashi during the 2008 flood period were first calculated by the 1D morphodynamic model. The results obtained at the specified sections were then used as the boundary conditions of the proposed 2D model. The computational mesh used in the 2D simulation consisted of 205 × 43 nodes, with a minimum mesh size of 50.5 m in the streamwise direction and of 15.1 m in the transverse direction.

Jing33: Jing34: Jing35: Jing36:

4.2.2. Comparisons of flow and sediment conditions The calculated and measured flow and sediment conditions during the 2008 flood period were compared at six fixed sections, as shown in Fig. 12. The calculations agreed well with the measurements, in terms of depth-averaged velocity, depth and SSC. For example, the calculated mean velocity and depth of the sections ranged from 0.98 m/s to 1.19 m/s and from 8.1 m to 11.85 m, respectively, at these sections, which were close to the measured values of 0.96–1.20 m/s and of 7.90–14.35 m. In addition, the calculated transverse distributions of depth-averaged velocity and water depth agreed well with the measurements. The agreement between the calculated and measured SSCs was relatively less satisfying, with the relative error of cross-section averaged SSC b10%, except for the values of 23.7% at the section of Sha4 and 33.3% at the section of Jing51. Fig. 13 indicates relatively close agreement between the calculated and measured hydrographs of daily mean water discharge and river stage at Shashi during the 2008 flood season, and indicates a slight underestimation of daily mean SSC. 4.2.3. Bank erosion Fig. 14 shows the initial (IBT) and final bank-top lines (FBT) and the flow field (at the calculated maximum discharge of 34,343 m3/s at the Shashi station) during the 2008 flood period. The bank retreat width was magnified by a factor of 10 for a clearer view, which led to the predicted FBT to exceed the channel boundary in Fig. 14. Calculated results indicate that obvious bank erosion occurred along the right bank between the sections of Sha4 and Jing36, where the main flow in the right branch approached, as shown in Fig. 14. The calculation slightly overestimated the area of bank erosion, as compared with the observations from the annually measured cross-sectional profiles, which indicated that bank erosion occurred along the right bank between the sections of Jing34 and Jing36. The calculated retreat widths (Wc) at Jing34 - Jing36 ranged from 21.4 m to 39.0 m during the study period, agreeing with the measured values (Wm) of 23.7–41.6 m in October 2007–October 2008.

Wm=0.0 m Wc=39.0 m Wm=41.6 m Wc=21.4 m Z : 5m 10m 15m 20m 25m 30m 35m 40m b Wm=23.7 m Wc=26.9 m Wm=28.8 m Wc=20.7 m

2m/s

Jing32

IBT FBT Sha4

Lanlinzhou

Jing33

Fig. 14. Comparison between the initial and final bank-top lines during the 2008 flood period.

20

S. Deng et al. / Geomorphology 344 (2019) 10–24

Jing34

Jing36

Jing42 Jing33 Jing35

Sanbatan Final central bar Initial central bar

Jinchengzhou

Zb: 0m 5 10 15 20 25 30 35 40 45 50

Fig. 15. Calculated evolution processes of central bars during the 2008 flood period, with the boundary of central bar characterized by the 30 m contour.

4.2.4. Evolution of central bars Calculated results also indicate that the central bar of Sanbatan was slightly eroded, and the central bars of Taipingkou and Jinchengzhou were slightly changed during the 2008 flood period, as shown in Fig. 15. However, observations from remote sensing images show that the central bar of Sanbatan was slightly eroded from March 2008 to February 2009, whereas obvious erosion occurred at Jinchengzhou and at the downstream end of Taipingkou, and significant aggradation occurred at the upstream end of Taipingkou (Fig. 16). Therefore, the proposed model seems to incorrectly predict the evolution trends of Taipingkou and Jinchengzhou during the 2008 flood period. The disagreement infers that the bed deformation equation (Eq. (8)) does not correctly characterize the evolution processes of central bars in the study reach, and thus it may be necessary to develop an improved module for the central bar evolution in the future.

the left branch, and thus led to high average values of 87.2% and 85.1% for RQL and of 93.7% and 90.2% for RQSL during these two flood seasons. According to the data by CWRC, the measured values of RQL and RQSL were 41.5% and 49.8%, respectively, in August 2004 near Sanbatan, and were 81.3% and 90.0% in June 2005 near Jinchengzhou. Furthermore, the measured annually mean values of RQL and RQSL during the flood seasons from 2003 to 2009, were 41.6% and 47.9% near Sanbatan and 91.0% and 91.1% near Jinchengzhou, respectively. Therefore, it can be concluded that the calculated flow and sediment diversion ratios were in good agreement with the measurements. In addition, calculated results indicate that: (i) the average flow and sediment diversion ratios for the left branch were decreased slightly in the 2008 flood period, as compared with the results during the 2004 flood period; and (ii) the right branch near Jinchengzhou was nearly dry (with RQL approaching 100%) when the total water discharge was b10,000 m3/s in October, as shown in Fig. 17e and f.

4.3. Water and sediment diversion ratios in the anabranching sub-reaches

4.4. Mass failure and groundwater level change

Fig. 17 shows the calculated water and sediment diversion ratios in the anabranching sub-reaches near the three central bars during the 2004 and 2008 flood periods. Fig. 17a and b indicate that: (i) near the central bar of Taipingkou, the calculated water diversion ratio for the left branch (RQL = QL/Q, where QL = water discharge entering the left branch; and Q = total water discharge) had an average value of 56.0% in 2004 and 53.3% in 2008, along with a sediment diversion ratio (RQSL = QsL/Qs, where QsL = the sediment discharge in the left branch; and Qs = the total sediment discharge) of 59.1% and 54.3%; (ii) near the central bar of Sanbatan, the calculated RQL was on average about 38.6% in both years, with RQSL at 43.9% in 2004 and 37.2% in 2008; and (iii) near the central bar of Jingchengzhou, the main flow went through

Fig. 18a and b, respectively, show the variations in the safety factor (Fs), average groundwater level (AGWL), river stage, and near-bank flow shear stress at the right bank of Jing34 in the 2004 and 2008 flood seasons. Bank failure occurred more frequently in 2008, mainly caused by the larger flow shear stresses. There was an obvious difference between the river stage hydrographs in these two flood seasons. In 2004, the rive stage hydrograph had a peak of about 41.8 m, with the flow overtopping the bank for seven days. However, the river stage hydrograph in 2008 had several peaks, with the highest peak reaching 40.16 m, indicating a more frequent and rapid change in river stage. Due to the difference in river stage, the AGWL in these two years changed differently. Generally, a more rapid change in river

(a)

Shashi sation

(b)

Sanbatan Taipingkou

Reference station˖ Shashi River stage: 29.98 m Time: 2008/3/26 Coordinate system˖ WGS84 49N

Boundary of central bar in 2008

Jinchengzhou

Reference station˖ Shashi River stage: 29.38 m Time: 2009/2/9 Coordinate system˖ WGS84 49N

Fig. 16. Remote sensing images taken on 26 Mar. 2008 and 9 Feb. 2009 with the river stages at Shashi of around 30 m.

Left diversion ratio (%)

Left diversion ratio (%)

Diversion ratio (%)

S. Deng et al. / Geomorphology 344 (2019) 10–24

80

Water

(a)

Sediment

60

70

21

(b)

55

60 50 50 40 55

45

(c)

55

50

50

45

45

40

40

35

35

30 100

(d)

30

(e)

100

95

95

90

90

85

85

80

80

75 1/Jul/04

31/Jul/04

30/Aug/04

29/Sep/04

(f)

75 1/Jul/08.

29/Oct/04

31/Jul/08. 30/Aug/08. 29/Sep/08. 29/Oct/08.

38

2.0

36 1.5 34 1.0

32

bank-toe

0.5

30 28 1/Jan. 42

Water level (m)

2.5

0.0 21/Jan.

10/Feb.

1/Mar.

21/Mar.

10/Apr.

(b) bankfull level

6.0

40

5.0

38

4.0

36

3.0

34

2.0

32 30 1/Jan.

1.0

bank-toe 21/Jan.

Failure 0.0 10/Feb.

1/Mar.

21/Mar.

10/Apr.

Shear stress (N/m2)

bankfull level

3.0

Safety factor

40

Water level (m)

River stage AGWL Safety factor Shear stress

(a)

Shear stress (N/m2)

42

Safety factor

Fig. 17. Calculated water and sediment diversion ratios in the anabranching sub-reaches near the central bars of: (a) Taipingkou in 2004; (b) Taipingkou in 2008; (c) Sanbatan in 2004; (d) Sanbatan in 2008; (e) Jinchengzhou in 2004; (f) Jinchengzhou in 2008.

Date Fig. 18. Changes in the safety factor (Fs), average groundwater level (AGWL), river stage, and near-bank flow shear stress at Jing34: (a) during the 2004 flood period; and (b) during the 2008 flood period.

22

S. Deng et al. / Geomorphology 344 (2019) 10–24

41

R1(Rq=0%)

AGWL (m)

40

R2(Rq=5%)

R3(Rq=10%)

R4(Rq=20%)

(a)

39 38 37 36

35 01/Jul. 100 (b)

21/Jul.

R2(Rq=5%)

70 Pw (%)

10/Aug.

30/Aug. R3(Rq=10%)

Jing 32 Jing 33

40

19/Sep.

09/Oct.

29/Oct.

R4(Rq=20%)

Jing 35

Jing 34

Jing 36

10 -20 143.3 143.8 144.7 145.3 145.5 145.8 146.1 146.4 146.8 147.1 147.5 147.8 Distance from the TGD (km) Fig. 19. Changes in the AWGL and bank retreat width in the cases of R1–R4: (a) AWGL; (b) retreat width.

stage would lead to a more delayed variation in AGWL. Therefore, a more obvious phase lag and a larger difference were observed between the river stage and the AGWL in 2008, as compared with the differences in 2004 (Fig. 18a). The delayed change in AWGL caused a higher pore water pressure relative to the hydrostatic confining force during the drawdown period, which was unfavorable to bank stability. Therefore, bank failures mostly occurred during the drawdown period (Fig. 18b), when there was a relatively higher pore water pressure and a rapidly decreasing hydrostatic confining force. 5. Discussion The effect of the infiltration recharge rate on bank erosion was investigated herein through numerical tests, along with the analysis of the impact of the incoming main flow shifting on the calculated result of central bar evolution. 5.1. Impact of infiltration recharge on bank erosion The aforementioned simulation in the 2008 flood period was referred to as R1, with three additional test runs (R2 - R4) conducted to investigate the impact of the infiltration recharge rate q in Eq. (9) on bank erosion. Rq (=q/kc) was set to be 0%, 5%, 10% and 20%, respectively, in the cases of R1 - R4. Fig. 19a shows the temporal changes in the AWGL at the right bank of Jing 34 in the cases of R1–R4. As indicated in Fig. 19a, the AWGL obviously increased with an increase in the infiltration

recharge rate, resulting in a higher pore water pressure, which was unfavorable to bank stability. Fig. 19b shows the percentage change (Pw) of bank retreat width in the cases of R2 - R4, relative to the retreat width in R1. It can be seen that: when considering the effect of infiltration recharge (q N 0), there was generally a larger bank retreat width (a positive value of Pw), with the average retreat width at sections between Jing32 - Jing36 increasing by 3.0% in R4 (Rq = 20.0%). However, Pw did not monotonically increase with Rq, and the retreat width even decreased (a negative Pw) with q N 0 at some sections. These phenomena indicate that even though the infiltration recharge, maybe due to the rainfall process, would generally lead to severe bank erosion, its exact impacts on bank erosion may be site specific. Furthermore, the decrease of bank retreat width may be explained as follows: (i) the increased infiltration recharge expedited the occurrence of bank failure, and thus might cause a larger number of bank failure events at some sections, but it led to a smaller volume of individual failed block material because of a lower fluvial erosion degree before failure; and (ii) the changes in the timing and angle of the initial bank failure, due to the increased infiltration recharge, influenced subsequent parallel failures. 5.2. Effect of the incoming flow condition on calculated central bar evolution Fig. 20 shows the re-calculated central bar evolution in the 2004 flood period, with a rightward (Fig. 20a) and leftward (Fig. 20b) shifting of the incoming main flow. Note that the discharges were the same in Sanbatan

Sanbatan Taipingkou Jing34

Velocity distribution of incoming flow on 7 th Oct.

Final central bar Initial central bar

Velovity (m/s)

Sha4 Jing33

(b) Jing42

Jing35

3.0

Taipingkou Jing34

Sha4 Jing33 Right-shifted Original

2.0

Velocity distribution of incoming flow on 7 th Oct.

1.0 0.0 300

Zb: 0m 5 10 15 20 25 30 35 40 45 50

700 1100 1500 Distance from the left bank (m)

Final central bar Initial central bar

Velovity (m/s)

(a)

Jing35 Left-shifted

2.0

Jing42 Original

1.0 0.0 300

700 1100 1500 Distance from the left bank (m)

Zb: 0m 5 10 15 20 25 30 35 40 45 50

Fig. 20. Calculated central bar evolution processes in the 2004 flood period, with the incoming main flow shifting toward: (a) the right; and (b) the left.

S. Deng et al. / Geomorphology 344 (2019) 10–24

(a)

(b)

3.0

Original

Velovity (m/s)

Right-shifted Left-shifted

Velovity (m/s)

3.0

23

Right-shifted Left-shifted

Original

2.0

2.0

1.0

1.0

0.0

0.0 0

400 800 1200 Distance from the left bank (m)

1600

0

400 800 1200 Distance from the left bank (m)

1600

Fig. 21. Transverse distributions of depth-averaged velocity on 7 October 2004, with shifted incoming main flow, at sections of: (a) Sha4; and (b) Jing35.

these cases. Comparing Fig. 20a with Fig. 10 demonstrates that the erosion was reduced at the upstream end and the aggradation was increased at the downstream end of Taipingkou, when the incoming flow was shifted to the right. However, a comparison between Fig. 20b and Fig. 10 shows the calculated erosion of Taipingkou was just slightly reduced with a leftward shifting of the incoming flow. In addition, the shifting of incoming flow did not have an obvious impact on the calculated adjustments of Sanbatan and Jinchengzhou. Fig. 21 shows the transverse distributions of depth-averaged velocity at sections of Sha4 and Jing35 on 7 October 2004, with a rightward and leftward shifting of incoming main flow. As illustrated in Fig. 21, the shifting of incoming flow exerted a greater impact on the velocity distribution at Sha4 (close to the inlet boundary). However, it had a relatively insignificant impact at Jing35. This indicates that the flow was enough developed in the sub-reach near Jing35 during the simulation period, which was thus slightly related to the incoming velocity distribution. Furthermore, because of the insignificant change in the velocity near Lalinzhou, the calculated bank erosion was also slightly influenced by the shifting of incoming main flow. Therefore, it can be concluded that the incoming main flow shifting would just have an obvious impact on the calculated channel evolution in the sub-reach near the inlet.

5.3. Effects of other factors on the calculated results Numerical tests were conducted in order to demonstrate the impacts of the transform proportion Rts and Reb in Eq. (13) on the calculated suspended sediment concentration (SSC) hydrograph. Results indicated that the effects of Rt and Reb were insignificant in the study reach. The daily averaged suspended sediment concentration at Shashi station was increased by b5%, as the value of Rt or Reb increased from 0.5 to 1.0. Based on the measurements, we evaluate potential explanations for the relative unsatisfying agreement between the calculated and measured SSCs. We believe that this disagreement may be partly caused by the relatively insufficient characterizations of sediment entrainment and deposition processes in the study reach. Table 2 shows the relationship between the measured depth-averaged SSC and the term U3/ghws in Eq. (5), which was calculated on the basis of relevant measurements in 2004 and 2008. As indicated in Table 2, SSC generally has a power re3

lationship with U =ghws . Assuming that SSC is linearly related to the sediment transport capacity S⁎ in the study reach undergoing a degradation process, it can be inferred that Eq. (5) is applicable to characterize S⁎

in the study reach. However, the exponent of the power relationship changed at different sections, implying that the parameter group (λ, m) in Eq. (5) probably also varied along the study reach, due to the changes in topographic and hydraulic conditions. However, in 2D modeling the variations in λ and m were difficult to be calibrated and were usually treated as constants. 6. Conclusion A coupled 2D model for channel evolution process was proposed by integrating the 2D morphodynamic module with the bank erosion module and the groundwater flow module. The variation of the groundwater table induced by changes in river stage was considered, which would lead to changes in the pore water pressure and the matrix suction for unsaturated bank soil. Model calibration and validation processes were conducted by comparing the calculated results with corresponding measurements. The effects of infiltration recharge on bank erosion and of incoming flow on central bar evolution were also investigated through numerical tests, along with the discussion of the reasons for the relatively unsatisfying disagreement between the calculated and measured suspended sediment concentration. Following conclusions can be draw from this study: (i) The proposed model shows a good performance in reproducing the flow and sediment transport processes, with the calculations agreeing well with the measurements in terms of depth, depthaveraged velocity and suspended sediment concentration, daily averaged hydrographs of discharge and river stage, and water and sediment diversion ratios in the local anabranching subreaches. In addition, the calculated degree of bank erosion was generally in agreement with the measurements. However, the central bar evolution was not correctly reproduced, inferring that an improved module needs to be developed in the future for the central bar evolution and its interaction with bank erosion. (ii) The groundwater level change was delayed as compared with the river stage change, and a more rapid variation in river stage would lead to a more obvious phase lag and a higher value (at the same river stage) of groundwater level. An increase in the infiltration recharge rate, due to rainfall or other processes, would lead to a higher groundwater level and thus generally increase the bank erosion degree in the study reach. However, at a

Table 2 3

Relationship between the measured SSC and U =ghws at different sections. 3

Sections

Best relationship with U =ghws

Exponent

Correlation Coefficient

Wan 15 Sha4 Jing35 Jing38 Jing45 Jing51

Power Power Power Power Power Linear

0.53 0.76 0.74 0.59 1.05 1.00

0.71 0.46 0.68 0.24 0.81 0.63

24

S. Deng et al. / Geomorphology 344 (2019) 10–24

specified section, the change trend in the degree of bank erosion, with an increased rate in infiltration recharge, was also related to other factors such as the timing of bank failure and the cumulative fluvial erosion action before failure. (iii) The incoming flow in the study reach had an obvious effect on the calculated evolution of the central bar that was close to the inlet boundary, with the leftward and rightward shifting of the incoming main flow (under the same discharge) causing different degrees of erosion to the nearby central bar. Acknowledgements This work was supported mainly by the National Natural Science Foundation of China (Grant nos. 51725902 & 51579186), and it was supported partly by the Program of the National Key Research and Development Plan (Grant nos. 2016YFC0402303/05). The constructive suggestions of the anonymous reviewers and the editors are gratefully acknowledged. References Abderrezzak, K.E.K., Moran, A.D., Tassi, P., Ata, R., Hervouet, J.M., 2016. Modelling river bank erosion using a 2D depth-averaged numerical model of flow and noncohesive, non-uniform sediment transport. Adv. Water Resour. 93, 75–88. https:// doi.org/10.1016/j.advwatres.2015.11.004. Camporeale, C., Perucca, E., Ridolfi, L., Gurnell, A.M., 2013. Modeling the interactions between river morphodynamics and riparian vegetation. Rev. Geophys. 51 (3), 379–414. https://doi.org/10.1002/rog.20014. Cantelli, A., Wong, M., Parker, G., Paola, C., 2007. Numerical model linking bed and bank evolution of incisional channel created by dam removal. Water Resour. Res. 43, W07436. https://doi.org/10.1029/2006WR005621. Chiang, S., Tsai, T., Yang, J., 2011. Conjunction effect of stream water level and groundwater flow for riverbank stability analysis. Environ. Earth Sci. 62 (4), 707–715. https:// doi.org/10.1007/s12665-010-0557-8. CWRC (Changjiang Water Resources Commission), 2014. Bank Erosion Investigation and Research in Key Reaches of Middle Yangtze River. Scientific Report of CWRC, Wuhan (in Chinese). Darby, S.E., Alabyan, A.M., Van de Wiel, M.J., 2002. Numerical simulation of bank erosion and channel migration in meandering rivers. Water Resour. Res. 38 (9), 1–2. https:// doi.org/10.1029/2001WR000602. Darby, S.E., Rinaldi, M., Dapporto, S., 2007. Coupled simulations of fluvial erosion and mass wasting for cohesive river banks. J. Geophys. Res. 112 (F3). https://doi.org/ 10.1029/2006JF000722. Fox, G.A., Wilson, G.V., 2010. The Role of subsurface flow in hillslope and stream bank erosion: a review. Soil Sci. Soc. Am. J. 74 (3), 717–733. https://doi.org/10.2136/ sssaj2009.0319. Gu, W.C., 2000. Principles and Applications of Seepage Modeling. China Building Material Press, Beijing (in Chinese). Hanson, G.J., Simon, A., 2001. Erodibility of cohesive streambeds in the loess area of the Midwestern USA. Hydrol. Process. 15 (1), 23–38. https://doi.org/10.1002/hyp.149. Jia, D.D., Shao, X.J., Wang, H., Zhou, G., 2010. Three-dimensional modeling of bank erosion and morphological changes in the Shishou bend of the middle Yangtze River. Adv. Water Resour. 33 (3), 348–360. https://doi.org/10.1016/j.advwatres.2010.01.002. Karmaker, T., Dutta, S., 2011. Erodibility of fine soil from the composite river bank of Brahmaputra in India. Hydrol. Process. 25 (1), 104–111. https://doi.org/10.1002/hyp.7826. Karmaker, T., Dutta, S., 2013. Modeling seepage erosion and bank retreat in a composite river bank. J. Hydrol. 476, 178–187. https://doi.org/10.1016/j.jhydrol.2012.10.032. Kimiaghalam, N., Goharrokhi, M., Clark, S.P., Ahmari, H., 2015. A comprehensive fluvial geomorphology study of riverbank erosion on the Red River in Winnipeg, Manitoba,

Canada. Journal of Hydrology 529 (3), 1488–1498. https://doi.org/10.1016/j. jhydrol.2015.08.033. Langendoen, E.J., Mendoza, A., Abad, J.D., Tassi, P., Wang, D., Ata, R., El kadi Abderrezzak, K., Hervouet, J.M., 2016. Improved numerical modeling of morphodynamics of rivers with steep banks. Adv. Water Resour. 93, 4–14. https://doi.org/10.1016/j. advwatres.2015.04.002. Li, Y.T., 1987. Preliminary study on the gradation of bed material load in equilibrium. Journal of Sediment Research 1, 82–87 (In Chinese). Li, Y.W., Xia, J.Q., Zhou, M.R., Deng, S.S., 2018. Analysis on evolution of mid-channel sandbars in Shashi reach after the Three Gorges Project operation. Journal of Hydroelectric Engineering 37 (10), 76–85 (In Chinese). Lindow, N., Fox, G.A., Evans, R.O., 2009. Seepage erosion in layered stream bank material. Earth Surf. Process. Landf. 34 (12), 1693–1701. https://doi.org/10.1002/esp.1874. Luppi, L., Rinaldi, M., Teruggi, L.B., Darby, S.E., Nardi, L., 2009. Monitoring and numerical modelling of riverbank erosion processes: a case study along the Cecina River (central Italy). Earth Surf. Process. Landf. 34 (4), 530–546. https://doi.org/10.1002/esp.1754. Masoodi, A., Noorzad, A., Tabatabai, M.M., Samadi, A., 2018. Application of short-range photogrammetry for monitoring seepage erosion of riverbank by laboratory experiments. J. Hydrol. 558, 380–391. https://doi.org/10.1016/j.jhydrol.2018.01.051. Morianou, G.G., Kourgialas, N.N., Karatzas, G.P., Nikolaidis, N.P., 2018. Assessing hydromorphological changes in Mediterranean stream using curvilinear grid modeling approach-climate change impacts. Earth Sci. Inf. 11 (2), 205–216. https://doi.org/ 10.1007/s12145-017-0326-2. Nardi, L., Campo, L., Rinaldi, M., 2013. Quantification of riverbank erosion and application in risk analysis. Nat. Hazards 69 (1), 869–887. https://doi.org/10.1007/s11069-0130741-8. Osman, A.M., Thorne, C.R., 1988. Riverbank stability analysis. I: Theory. J. Hydraul. Eng. 114 (2), 134–150. https://doi.org/10.1061/(ASCE)0733-9429(1988)114:2(134). Rinaldi, M., Casagli, N., 1999. Stability of streambanks formed in partially saturated soils and effect of negative pore water pressures: the Sieve River (Italy). Geomorphology 26 (4), 253–277. https://doi.org/10.1016/S0169-555X(98)00069-5. Rinaldi, M., Nardi, L., 2013. Modeling interactions between riverbank hydrology and mass failures. J. Hydrol. Eng. 18 (10), 1231–1240. https://doi.org/10.1061/(ASCE)HE.19435584.0000716. Rinaldi, M., Mengoni, B., Luppi, L., Darby, S.E., Mosselman, E., 2008. Numerical simulation of hydrodynamics and bank erosion in a river bend. Water Resour. Res. 44, W09428. https://doi.org/10.1029/2008WR007008. Simon, A., Curini, A., Darby, S.E., Langendoen, E.J., 2000. Bank and near-bank processes in an incised channel. Geomorphology 35 (3–4), 193–217. https://doi.org/10.1016/ S0169-555X(00)00036-2. Simon, A., Pollen, N., Langendoen, E., 2006. Influence of two woody riparian species on critical conditions for streambank stability: Upper Truckee River, California. JAWRA Journal of the American Water Resources Association 42 (1), 99–113. https://doi. org/10.1111/j.1752-1688.2006.tb03826.x. Thorne, C.R., Alonso, C., Borah, D., Darby, S., Diplas, P., Julien, P., Knight, D., Pizzuto, J., Quick, M., Simon, A., 1998. River width adjustment. I: Processes and mechanisms. J. Hydraul. Eng. 124 (9), 881–902. https://doi.org/10.1061/(ASCE)0733-9429(1998) 124:9(881). Wang, Y.G., Kuang, S.F., 2007. Critical height of bank collapse. J. Hydraul. Eng. 38 (10), 1158–1165 (in Chinese). Wang, G., Xia, J., Wu, B., 2008. Numerical simulation of longitudinal and lateral channel deformations in the braided reach of the lower Yellow River. J. Hydraul. Eng. 134 (8), 1064–1078. https://doi.org/10.1061/(ASCE)0733-9429(2008)134:8(1064. Xia, J.Q., Zong, Q.L., 2015. Mechanisms and Numerical Modeling of Bank Erosion in the Jingjiang Reach of the Middle Yangtze River. Science Press, Beijing (in Chinese). Xia, J.Q., Wang, Z.B., Wang, Y., Yu, X., 2013. Comparison of morphodynamic models for the Lower Yellow River. JAWRA Journal of the American Water Resources Association 49 (1), 114–131. https://doi.org/10.1111/jawr.12002. Xia, J.Q., Zong, Q.L., Deng, S.S., Xu, Q.X., Lu, J.Y., 2014. Seasonal variations in composite riverbank stability in the Lower Jingjiang Reach, China. J. Hydrol. 519, 3664–3673. https://doi.org/10.1016/j.jhydrol.2014.10.061. Xia, J.Q., Deng, S.S., Lu, J.Y., Xu, Q.X., Zong, Q.L., Tan, G.M., 2016. Dynamic channel adjustments in the Jingjiang Reach of the Middle Yangtze River. Sci. Rep. 6, 22802. https:// doi.org/10.1038/srep22802. Zhang, R.J., 1998. River and Sediment Dynamics. China Water and Power Press, Beijing (in Chinese).