# Evolution Mechanism of Meandering River Downstream Gigantic Hydraulic Project I: Hydrodynamic Models and Verification

Authors: Linshuang Liu; Haiqing Zhu; Chengtao Huang; Li Zheng
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The evolution mechanism of meandering river is one of essential references to predict the evolution disciplines of meandering river. Jingjiang curved reaches are the typical meandering river; more specifically, they are located downstream the gigantic hydraulic project named Three Gorge Project (TGP). Because the incoming water and sediment condition have been changed by the gigantic project, the evolution behavior of Jingjiang curved reaches changes a lot, making the evolution behavior unpredictable. However, traditional two-dimensional (2D) hydrodynamic model could not simulate the transportation characteristics of unbalanced and suspended sediment, leading the predict results of 2D model are far from the measured data. This paper presents a theoretical and numerical approach that explores the evolution mechanism of meandering river downstream gigantic hydraulic project. Firstly, the evolution behavior and evolution disciplines of Jingjiang curved reaches were classified before and after the hydraulic project implement respectively. Secondly, a 2D hydrodynamic model was set up and verified according to the measured data. And then a superior three-dimensional (3D) numerical model, whose boundary conditions were simulated by the results of the 2D numerical model, considering the unbalanced water and sediment transportation properties, was developed and verified by the measured data. Research results show that the 2D model displayed a reasonable accuracy in predicting the water level, branch diversion ratio, and flow velocity; the 3D model displayed a better accuracy in predicting the water lever, vertical flow velocity, longitudinal flow velocity, sediment concentration, and sediment variable quantity. Both 2D and 3D models could be applied to study the evolution mechanism of meandering river; especially the proposed 3D model considering the sediment transport in longitudinal, transverse, and vertical directions will improve the accuracy of behavior prediction and will help decision-making for the river regulation.


---

## Body

## 1. Introduction

TGP is the largest water conservation project in human history. It started on December 14, 1994; test impoundment began on June 1, 2003, and impounded to the normal water level (means above the sea level 175 m) on October 26, 2010. It provides comprehensive benefits including flood control, hydropower, shipping, and water supply [1]. At the meanwhile, the great project changes the intrinsic circumstance of water and sediment, which has a substantial influence on the evolution laws of the downstream reaches, especially reaches near TGP [2].Jingjiang reach, in middle Yangtze River, is a typical meandering river, and it is the nearest meandering reach downstream the TGP. Jingjiang reach consists of 16 local curved reaches, shown as Figure1, ① Guanzhou curved reach (CR), ② Jiangkou CR, ③ Yuanshi CR, ④ Shashi CR, ⑤ Maijiazhai CR, ⑥ Haoxue CR, ⑦ Ouchikou CR, ⑧ Nianziwan CR, ⑨ Tiaoguan CR, ⑩ Laijiapu CR, ⑪ Jianli CR, ⑫ Zhuanqiao CR, ⑬ Fanzui CR, ⑭ Xiongjiazhou CR, ⑮ Qigongling CR, and⑯ Guanyingzhou CR.Figure 1
Schematic of the covered reach in the downstream of the TGD.Under the natural condition, curved reach follows the traditional evolution disciplines including the following: (1) in the low water period, the main stream will go along the river bank; (2) in the high water period, the main stream will go downstream directly; (3) in the long term, concave bank will be scoured while convex bank will be deposited [3–7] (Schumm, 1960). Jingjiang reach followed the above disciplines until the TGD was established.Present researches [8–10] and measured data (shown in Figure 2) showed that, after Three Gorge Reservoir impoundment, water and sediment go downstream changed and brought about new evolution disciplines such like: (1) in the flood period, flood peak decreased; (2) in the dry period, total water 65 flow increased; (3) the duration of normal water period increased; (4) sediment went to the downstream driver decreased dramatically, usually decreased more than 80%. In recent years, we have observed that in Jingjiang reach, some concave banks begin to deposit while some convex banks begin to scour which present an opposite disciplines to the traditional disciplines. For example, in ⑤ Maijiazhai CR and ⑧ Nianziwan CR, the convex bank beach has been scoured while the concave bank bed has been deposited obviously. The concave bank beds of ⑨ Tiaoguan CR and ⑬ Fanzui CR have been deposited severely, forming a new diara near the concave bank. Conversely, there ought to be a deposit at the convex of ⑮ Qigongling CR, but it has been scoured severely, and the convex bank bed has been cut into two parts by a small ditch at the end of the convex band bed. Therefore, river beach and groove distribution tend to be unpredictable and would impact the flood control, irrigation, and shipping of the downstream of Yangtze River.Figure 2
Evolution history of Jingjiang reach downstream of the TGD.In brief, water and sediment transport of curved reach in longitudinal, transverse, and vertical directions are strong. However, the traditional two-dimensional (2D) numerical model could not work very well on the complicated dynamic circumstance.This paper presents a theoretical and numerical study on the evolution disciplines of the downstream Yangtze River. Firstly, the evolution properties were summarized and classified according to history measured data. Then a two-dimensional (3D) numerical model was set up and verified according to the measured data. At last a superior three-dimensional (3D) numerical model, whose boundary conditions were simulated by the results of the 2D numerical model, considering the unbalanced water and sediment transportation properties, was developed and verified by the measured data. Results show that the proposed 3D model displayed a reasonable agreement on measured data and could be used to study the evolution mechanism of meandering river.


## 2. Evolution Properties of Meandering River

### 2.1. Before the TGD Impounded Water

Jingjing reach originated from Yunmeng settlement in the Tertiary period, and then the channel changed frequently under the natural circumstance with complicated mechanism. As recorded literature [11–13], it experienced four periods: (1) formation of the delta bed; (2) attenuation of the delta bed; (3) formation of the meandering reach; and (4) propagation of the meandering reach.Jingjiang reach can be divided into two parts at Ouchikou: (1) upper-Jingjiang reach, between Zhicheng and Ouchikou, and (2) lower-Jingjiang reach, between Ouchikou to Chenglingji [14]. In the last three hundred years, upper-Jingjiang reach evolved lightly while lower-Jingjiang reach evolved severely, shown as Figure 2. It shows that upper-Jingjiang reach wiggled its river bed gradually, but the river channel kept its general configuration. Lower-Jingjiang reach wiggled its river bed frequently, forming and attenuating meanders, making the river alternately go downstream curvedly and directly. In order to control the original propagation, a series of drainage valves, levee construction, and bank protection projects had been set up along Jingjiang reach. Till 1970s, the artificial facilities have made great difference for slowing down the wigging of the river bed [11, 15]. For example, because of man-made cut-bend projects at Zhongzhouzi and Shangchewan (shown in Figure 2, in the picture 1973 AD), along with the effect of natural cut-bend at Shatanzi, wigging of lower-Jingjiang reach has been controlled generally.Before impounding of TGR (from 1970s to year 2003), in the low water period, the main stream in curved reaches went along the channel bank; in the high water period, the main stream went downstream directly [11, 14, 16]. It followed the traditional evolution disciplines: concave bank scoured while convex bank deposited, shown as Figures 3(a) and 3(b). In general, the evolution disciplines of Jingjiang reach presented two types:Figure 3
Cross section change of curving reaches in Jingjiang river.
(a)
Section position map (b)
Cross section change before impoundment of the Three Gorges Reservoir (c)
Cross section change after impoundment of the Three Gorges Reservoir(1) Type A, concave bank was scoured while convex bank was deposited. For most curving reaches, concave bank was scoured while convex bank was deposited; as a result meandering radius of the river will decline gradually, such as ④ Shishou CR, ⑨ Tiaoguan CR, ⑩ Laijiapu CR, ⑪ Jianli CR, and ⑬ Fanzui CR.(2) Type B, convex bank was scoured and will even be cut-bend. For minority meandering reaches, under the effects of natural incoming water and river bank boundaries, the convex was scoured and even a little groove could propagate to be a natural cut-bend channel, making the meandering reach into straight reach, such as ⑭ Xiongjiazhou CR, ⑮ Qigongling CR, and⑯ Guanyingzhou CR.


### 2.2. After the TGD Impounded Water

The incoming water and sediment changed dramatically after impoundment of TGR [15, 17–19]; consequently the evolution disciplines of the downstream curved reach changed as follows.(1) Type A has changed into an opposite appearance, which means concave bank was deposited while convex bank scoured, such as Shashi CR, Hekou CR, Laijiapu CR, Faizui CR, and Xingjiazhou CR, shown as Figures3(a) and 3(c).Taking Guanzhou CR (a typical Type A reach) for instance, bed near convex bank had been scoured obviously and bed near concave bank had been deposited obviously after impounding of TGR, for the average depth of scouring grooves were 5 m, and 3 m for depositing area near concave bank, shown as Figure4. The distributary ratios of convex at Guanzhou CR changed from 19.07% (at year 2003) to 34.08% (at year 2012).Figure 4
Erosion and siltation distribution of ①Guanzhou CR (2003.02~2012.02).As well as ⑨Tiaoguan CR, its convex bank (named Jijiazui beach) was scoured badly; the average scouring depth was more than 5 m, shown as Figure5. Its maximum distance of zero-meter isobath wentting to the convex bank were 230m (2003~2006) and 420m (2006~2009). The zero-meter isobath area decreased from 1.88 (km)2 (2002) to 0.56 (km)2 (2009), shrinking almost two thirds. Because of the scouring effect, the area of Jijiazui beach became smaller and smaller. Therefore, water flow in ⑨Tiaoguan CR got diminishing constraint, tended to go straight along the river channel, which made the nearby concave bank (Xiaojiaguai) became deposited. There had been an obvious sediment deposition in nearby concave bank in August, 2009, and it became a diara in August, 2014. In the following years, the diara is steady and divides the main river channel into two channels in low water period, shown as Figure 6.Figure 5
Erosion and siltation distribution of ⑨Tiaoguan CR (2002.10~2009.09).Figure 6
Isobath variation in ⑨Tiaoguan CR. (a) zero-meter and (b) four-meter.
(a) (b)(2) Type B has kept its appearance but accelerates its speed or aggravate its level, such as Wakouzi CR, Majiazui CR and Qigongling CR; Distributary ratios near convex bank at Wakouzi CR and Majiazui CR has been increased; Little channeling groove near convex bank at Qigongling CR and Guanyinzhou CR has been developed, even separates the beach into two parts.Taking Qigongling CR (a typical Type B reach) for instance, the convex bank already be scoured and resulted in little channeling grooves before impounding of TGR; and then the beach became shrunken rapidly and the grooves became extended rapidly after the TGR impounding water, shown as Figure7. From Tables 1 and 2, beach area at the convex bank was 1.7 (km)2 (at year 2004) but dropped to 0.28 (km)2 (at year 2012). The river island area at the concave bank was 0.07 (km)2 (at year 2004) but soared to 1.92 (km)2 (at year 2012).Table 1
The geometric appearance of the convex beach at Qigongling CR (the length/width is the maximum value of the beach).
Measured time Area (km2) Length (m) Width (m) 2006.04.23 1.70 3616 662 2007.04.24 0.75 2333 394 2009.03.29 0.48 1236 454 2009.09.06 0.43 1181 385 2009.11.13 0.28 892 377 2012.02.15 0.28 912 184Table 2
The geometric appearance of the river island at Qigongling CR (the length/width is the maximum value of the river island).
Measured time Area(km2) Length (m) Width (m) 2004.03.12 0.11 1300 120 2006.04.23 1.19 5858 354 2007.04.24 1.79 6482 533 2009.03.29 2.31 3481 866 2009.09.06 2.16 3297 873 2009.11.13 1.69 2355 1139 2012.02.15 1.92 3520 850Figure 7
Erosion and siltation distribution of ⑭ Xiongjiazhou CR, ⑮ Qigongling CR, and⑯ Guanyingzhou CR (2002.02~2010.02).


## 2.1. Before the TGD Impounded Water

Jingjing reach originated from Yunmeng settlement in the Tertiary period, and then the channel changed frequently under the natural circumstance with complicated mechanism. As recorded literature [11–13], it experienced four periods: (1) formation of the delta bed; (2) attenuation of the delta bed; (3) formation of the meandering reach; and (4) propagation of the meandering reach.Jingjiang reach can be divided into two parts at Ouchikou: (1) upper-Jingjiang reach, between Zhicheng and Ouchikou, and (2) lower-Jingjiang reach, between Ouchikou to Chenglingji [14]. In the last three hundred years, upper-Jingjiang reach evolved lightly while lower-Jingjiang reach evolved severely, shown as Figure 2. It shows that upper-Jingjiang reach wiggled its river bed gradually, but the river channel kept its general configuration. Lower-Jingjiang reach wiggled its river bed frequently, forming and attenuating meanders, making the river alternately go downstream curvedly and directly. In order to control the original propagation, a series of drainage valves, levee construction, and bank protection projects had been set up along Jingjiang reach. Till 1970s, the artificial facilities have made great difference for slowing down the wigging of the river bed [11, 15]. For example, because of man-made cut-bend projects at Zhongzhouzi and Shangchewan (shown in Figure 2, in the picture 1973 AD), along with the effect of natural cut-bend at Shatanzi, wigging of lower-Jingjiang reach has been controlled generally.Before impounding of TGR (from 1970s to year 2003), in the low water period, the main stream in curved reaches went along the channel bank; in the high water period, the main stream went downstream directly [11, 14, 16]. It followed the traditional evolution disciplines: concave bank scoured while convex bank deposited, shown as Figures 3(a) and 3(b). In general, the evolution disciplines of Jingjiang reach presented two types:Figure 3
Cross section change of curving reaches in Jingjiang river.
(a)
Section position map (b)
Cross section change before impoundment of the Three Gorges Reservoir (c)
Cross section change after impoundment of the Three Gorges Reservoir(1) Type A, concave bank was scoured while convex bank was deposited. For most curving reaches, concave bank was scoured while convex bank was deposited; as a result meandering radius of the river will decline gradually, such as ④ Shishou CR, ⑨ Tiaoguan CR, ⑩ Laijiapu CR, ⑪ Jianli CR, and ⑬ Fanzui CR.(2) Type B, convex bank was scoured and will even be cut-bend. For minority meandering reaches, under the effects of natural incoming water and river bank boundaries, the convex was scoured and even a little groove could propagate to be a natural cut-bend channel, making the meandering reach into straight reach, such as ⑭ Xiongjiazhou CR, ⑮ Qigongling CR, and⑯ Guanyingzhou CR.


## 2.2. After the TGD Impounded Water

The incoming water and sediment changed dramatically after impoundment of TGR [15, 17–19]; consequently the evolution disciplines of the downstream curved reach changed as follows.(1) Type A has changed into an opposite appearance, which means concave bank was deposited while convex bank scoured, such as Shashi CR, Hekou CR, Laijiapu CR, Faizui CR, and Xingjiazhou CR, shown as Figures3(a) and 3(c).Taking Guanzhou CR (a typical Type A reach) for instance, bed near convex bank had been scoured obviously and bed near concave bank had been deposited obviously after impounding of TGR, for the average depth of scouring grooves were 5 m, and 3 m for depositing area near concave bank, shown as Figure4. The distributary ratios of convex at Guanzhou CR changed from 19.07% (at year 2003) to 34.08% (at year 2012).Figure 4
Erosion and siltation distribution of ①Guanzhou CR (2003.02~2012.02).As well as ⑨Tiaoguan CR, its convex bank (named Jijiazui beach) was scoured badly; the average scouring depth was more than 5 m, shown as Figure5. Its maximum distance of zero-meter isobath wentting to the convex bank were 230m (2003~2006) and 420m (2006~2009). The zero-meter isobath area decreased from 1.88 (km)2 (2002) to 0.56 (km)2 (2009), shrinking almost two thirds. Because of the scouring effect, the area of Jijiazui beach became smaller and smaller. Therefore, water flow in ⑨Tiaoguan CR got diminishing constraint, tended to go straight along the river channel, which made the nearby concave bank (Xiaojiaguai) became deposited. There had been an obvious sediment deposition in nearby concave bank in August, 2009, and it became a diara in August, 2014. In the following years, the diara is steady and divides the main river channel into two channels in low water period, shown as Figure 6.Figure 5
Erosion and siltation distribution of ⑨Tiaoguan CR (2002.10~2009.09).Figure 6
Isobath variation in ⑨Tiaoguan CR. (a) zero-meter and (b) four-meter.
(a) (b)(2) Type B has kept its appearance but accelerates its speed or aggravate its level, such as Wakouzi CR, Majiazui CR and Qigongling CR; Distributary ratios near convex bank at Wakouzi CR and Majiazui CR has been increased; Little channeling groove near convex bank at Qigongling CR and Guanyinzhou CR has been developed, even separates the beach into two parts.Taking Qigongling CR (a typical Type B reach) for instance, the convex bank already be scoured and resulted in little channeling grooves before impounding of TGR; and then the beach became shrunken rapidly and the grooves became extended rapidly after the TGR impounding water, shown as Figure7. From Tables 1 and 2, beach area at the convex bank was 1.7 (km)2 (at year 2004) but dropped to 0.28 (km)2 (at year 2012). The river island area at the concave bank was 0.07 (km)2 (at year 2004) but soared to 1.92 (km)2 (at year 2012).Table 1
The geometric appearance of the convex beach at Qigongling CR (the length/width is the maximum value of the beach).
Measured time Area (km2) Length (m) Width (m) 2006.04.23 1.70 3616 662 2007.04.24 0.75 2333 394 2009.03.29 0.48 1236 454 2009.09.06 0.43 1181 385 2009.11.13 0.28 892 377 2012.02.15 0.28 912 184Table 2
The geometric appearance of the river island at Qigongling CR (the length/width is the maximum value of the river island).
Measured time Area(km2) Length (m) Width (m) 2004.03.12 0.11 1300 120 2006.04.23 1.19 5858 354 2007.04.24 1.79 6482 533 2009.03.29 2.31 3481 866 2009.09.06 2.16 3297 873 2009.11.13 1.69 2355 1139 2012.02.15 1.92 3520 850Figure 7
Erosion and siltation distribution of ⑭ Xiongjiazhou CR, ⑮ Qigongling CR, and⑯ Guanyingzhou CR (2002.02~2010.02).


## 3. Two-Dimensional Model

The evolution disciplines of curved reaches were complex. To improve the analysis accuracy, a two-dimensional (2D) numerical model was built. The model was verified by the measured data and used to study the evolution disciplines of curved reaches.


### 3.1. Fundamental Equations

A pane two-dimensional mathematical model was use to simulated the river [19, 20]. The transformation of normal curve equations was shown as follows from (1) to (3).Flow continuous equation is(1)J∂Z∂t+∂HU∂ζ+∂HV∂η=0where J was water surface slope;Z was water level, m;H is water depth, m;U was velocity in X direction, m/s; V was velocity in X direction, m/s; ξ was direction of river mainstream; η was direction of river width.Flow motion equation is(2)J∂M∂t+∂MU∂ζ+∂MV∂η=-ghJξx∂Z∂ζ+ηx∂Z∂η+∂∂ζDJq11Mζ+q12Mη+∂∂ηDJq12Mζ+q22Mη-gn2Mu2+v2h4/3JwhereM andN were discharged per unit width inx andy direction separately,M=uhandN=vh;n was Manning Roughness Coefficient;g was gravity acceleration, m/s2;D was Turbulent Viscosity Coefficient,L2T−1; ξx,ηx,ξy, and ηy were the conversion coefficient of X-ξ, X-η, Y-ξ, and Y-η. q11=ξx2+ξy2, q12=ξxηx+ξyηy, and q22=ηx2+ηy2; Mx, My,Nx, andNy were partial derivative, asMx wasM’s partial derivative in X direction.Motion equation of the water goes vertically along the river channel:(3)J∂N∂t+∂NU∂ζ∂NV∂η=-ghJζy∂Z∂ζ+ηy∂Z∂η+∂∂ζDJq11Nζ+q12Nη+∂∂ηDJq12Nζ+q22Nη-gn2Nu2+v2h4/3J


### 3.2. Roughness Factor

The channel bed resistance includes sand-grain resistance, sand-wave resistance, bed configuration resistance, and settled deposit resistance. Sections of Jingjiang channel were wide and deep, and channel bed was made of sand. An empirical equation for roughness factor is chosen as(4)n=C1·1+1hn1·fA,r,uwhere n is roughness factor of the river channel; C1 is unidimensional roughness factor; h is water depth; n1 is roughness factor of the across section; A is water flow area; r is radius curvature; u is current velocity; f(A,r,u) is the plane morphological function of the river; the value is generally between 0.5 and 1.5 [21].


### 3.3. Computational Meshes

The computational results of 2D model will be used to simulate the boundary conditions of 3D model, while the boundary conditions of 2D model are obtained from the observed data and its interpolation results. Taking Nandiguai-Tashiyi reach (including Hekou CR, Tiaoguan CR and Laijiapu CR) and Tashiyi-Chenglingji reach (including Jianli CR, Tianxingge CR, Fanzui CR, Xiongjiazhou CR, Qigongling CR, Guanyingzhou CR) for instances, the mesh details are as follows:(  1) Tianguan-Tashiyi Reach. Nandiguai-Tashiyi reach begins at Nandiguai and ends at Tashiyi, while the total length is 37 km. The meshed model is shown as Figure 8, and there are a total of 520×120 grid points. The type of the grid is calculated from Elliptic differential equations [22]; along the river channel there are 520 grid points and the distance between every two grid points is 40~100m; furthermore, perpendicular to the river channel there are 100 grid points and the distance between every two grid points is 10~20m.Figure 8
Mesh details of Nandiguai-Tashiyi reach.(  2) Tashiyi-Chenglingji Reach.Tashiyi-Chenglingji reach begins at Tashiyi and ends at Chenglingji, while the total length is 96 km. The meshed model is shown as Figure 9, and there are a total of 900×80 grid points. The type of the grid is the orthogonal curvilinear coordinate mesh [23]; along the river channel there are 900 grid points and the distance between every two grid points is 67~144 m; furthermore, perpendicular to the river channel there are 80 grid points and the distance between every two grid points is 5~40 m.Figure 9
Mesh details of Tashiyi-Chenglingji reach.


### 3.4. Verification

#### 3.4.1. Verification of Fixed Bed Model

River toughness factor and sand transport capacity are calibrated according to the measured data, including water level, distributary ratio of branch, current velocity distribution, and river scouring and depositing amount.(  1) Water Level. Water levels at different locations are calculated by the 2D model at flood period, median water period, and drought period, respectively. With the measured data, the calculating results are shown in Tables 3, 4, and 5 accordingly. All of the water level errors are less than 0.11m. Predicted water level displays a reasonable agreement with measured data.Table 3
Water level at flood period (Q is flow rate, hm is measured value, hp is predicted value, and γh is error, γh=hp-hm).
Sectionlocation Jianli CR @ July, 2007 Fanzui CR @ Aug, 2010 Tiaoguan CR @ Sep, 2010 Q =28900 m3/s Q =20444 m3/s Q =19735 m3/s h m h c γ h h m h c γ h h m h c γ h (m) (m) (m) (m) (m) (m) (m) (m) (m) 1#L 31.46 31.5 0.04 29.02 29.00 -0.02 32.15 32.16 0.00 2#L 31.26 31.37 0.11 28.88 28.86 -0.02 31.92 31.94 0.02 3#L 31.1 31.1 0.00 28.70 28.67 -0.03 31.77 31.78 0.01 4#L 31.1 31.1 0.00 28.50 28.51 0.01 31.62 31.62 0.00 5#L 30.99 30.97 -0.02 28.35 28.40 0.05 31.46 31.43 -0.03 6#L / / / 28.15 28.16 0.01 31.33 31.29 -0.04 7#L / / / / / / 31.04 31.02 -0.02Table 4
Water level at median water period.
Sectionlocation Jianli CR @ Oct, 2008 Fanzui CR @ Sep, 2009 Q =16322 m3/s Q =14232 m3/s h m h c γ h h m h c γ h (m) (m) (m) (m) (m) (m) 1#L 29.45 29.48 0.03 27.12 27.18 0.06 2#L 29.33 29.33 0.00 26.99 27.05 0.06 3#L / / / 26.84 26.88 0.04 4#L 29.21 29.22 0.01 26.69 26.66 -0.03 5#L 28.92 28.87 -0.05 26.62 26.58 -0.04 7#L    26.42 26.37 -0.05Table 5
Water level at drought period.
Sectionlocation Jianli CR @ Jun, 2006 Qigongling CR @ Feb, 2012 Tiaoguan CR @ Feb, 2014 Q =5231 m3/s Q =6262 m3/s Q =6518 m3/s h m h c γ h h m h c γ h h m h c γ h (m) (m) (m) (m) (m) (m) (m) (m) (m) 1#L 23.11 23.17 0.06 29.02 29.00 -0.02 24.40 24.41 -0.01 2#L 22.94 22.96 0.02 28.88 28.86 -0.02 24.19 24.20 -0.01 3#L 22.57 22.52 -0.05 28.70 28.67 -0.03 24.06 24.08 -0.02 4#L 22.22 22.25 0.04 28.50 28.51 0.01 23.99 23.98 0.01 5#L 22.14 22.14 0 28.35 28.40 0.05 23.84 23.83 0.01 6#L 23.11 23.17 0.06 28.15 28.16 0.01 23.68 23.71 -0.03 7#L / / / / / / 23.51 23.56 -0.05(  2) Distributary Ratio of Branch.Jianli CR as Figure 10, a typical branching curve, is chosen to validate distributary ratio of the model. Distributary ratios in Jianli CR were calculated by the 2D model at flood period (July 2007, Q=28900 m3/s), median water period (Oct, 2008, Q=16322m3/s; Sep, 2010, Q=19175m3/s), and drought period, respectively (Jan 2006, 5231m3/s). From the results, the calculated values of each branch diversion ratio are in good agreement with the measured values and errors are within 2%. Results show that the model can reflect variation of branch diversion and meet the requirement of calculation. The split ratio verification is shown in Table 6.Table 6
Distributary ratio at different flow rate (φm is measured value, φp is predicted value, and γφ is error, γφ=φp-φm).
Q (m3/s) Date Left branch Right branch φ m(%) φ p(%) γ φ(%) φ m(%) φ p(%) γ φ(%) Low water period 5231 2006.01 3.76 3.14 -0.62 96.24 96.86 0.62  Median water period 16322 2008.10 8.50 10.36 1.86 91.50 89.64 -1.86 19175 2010.09 10.77 10.22 -0.55 89.23 89.78 0.55  High water period 28900 2007.07 8.96 8.02 -0.94 91.04 91.98 0.94Figure 10
Sketch map of Jianli CR (December, 2014).(  3) Flow Velocity Distribution Verification. Calculated flow velocity values in different curved reach are vivificated by measured values. Typical sections in Tiaoguan CR and Laijiapu CR are chosen; comparison diagram with calculated values and measured values are drew in Figure 11. All of the errors between measured value and measured values are less than 0.2m/s, and the model meets the requirements of flow velocity simulation.Figure 11
Cross section velocity verification of Tiaoguan CR on December, 2014 ( —calculated value, ■measured data): (a) Q=7633 m3/s; (b) Q=19263m3/s,.
(a) (b)


#### 3.4.2. Verification of Moving Bed Model

According to the field testing points, reaches listed in Table6 are selected to compare the measured and calculated variable quantity of sediment. The geomorphic features in March, 2010 are set as the boundary conditions for the 2D numerical model; the incoming sediment and water are considered as extra loading. Table 7 shows that the 2D numerical model has a reasonable accuracy, for the maximum error is 15.2%. Therefore, the 2D numerical model can be used as a technique to study the dynamic evolution disciplines downstream of the TGD.Table 7
Comparison of depositing and scouring amount in selected reaches (from 2010.3 to 2012.2). (α is the measured data, β is the calculated value, positive values mean depositing, and minus values mean scouring. γ is the relative error, γ=β-α/α×100%).
Reach name α β γ (×104 m3) (×104 m3) (%) Tiaoguan CR 163.4 185.4 13.46 Laijiapu CR 351.9 303.8 -13.67 Tashiyi to Hongshuigang reach -549.58 -485.03 -11.7 Yanchuantao to Jingjiangmen reach -1030.35 -873.62 -15.2 Xiongjiazhou to Chenglingji reach 855.98 758.24 -11.4


## 3.1. Fundamental Equations

A pane two-dimensional mathematical model was use to simulated the river [19, 20]. The transformation of normal curve equations was shown as follows from (1) to (3).Flow continuous equation is(1)J∂Z∂t+∂HU∂ζ+∂HV∂η=0where J was water surface slope;Z was water level, m;H is water depth, m;U was velocity in X direction, m/s; V was velocity in X direction, m/s; ξ was direction of river mainstream; η was direction of river width.Flow motion equation is(2)J∂M∂t+∂MU∂ζ+∂MV∂η=-ghJξx∂Z∂ζ+ηx∂Z∂η+∂∂ζDJq11Mζ+q12Mη+∂∂ηDJq12Mζ+q22Mη-gn2Mu2+v2h4/3JwhereM andN were discharged per unit width inx andy direction separately,M=uhandN=vh;n was Manning Roughness Coefficient;g was gravity acceleration, m/s2;D was Turbulent Viscosity Coefficient,L2T−1; ξx,ηx,ξy, and ηy were the conversion coefficient of X-ξ, X-η, Y-ξ, and Y-η. q11=ξx2+ξy2, q12=ξxηx+ξyηy, and q22=ηx2+ηy2; Mx, My,Nx, andNy were partial derivative, asMx wasM’s partial derivative in X direction.Motion equation of the water goes vertically along the river channel:(3)J∂N∂t+∂NU∂ζ∂NV∂η=-ghJζy∂Z∂ζ+ηy∂Z∂η+∂∂ζDJq11Nζ+q12Nη+∂∂ηDJq12Nζ+q22Nη-gn2Nu2+v2h4/3J


## 3.2. Roughness Factor

The channel bed resistance includes sand-grain resistance, sand-wave resistance, bed configuration resistance, and settled deposit resistance. Sections of Jingjiang channel were wide and deep, and channel bed was made of sand. An empirical equation for roughness factor is chosen as(4)n=C1·1+1hn1·fA,r,uwhere n is roughness factor of the river channel; C1 is unidimensional roughness factor; h is water depth; n1 is roughness factor of the across section; A is water flow area; r is radius curvature; u is current velocity; f(A,r,u) is the plane morphological function of the river; the value is generally between 0.5 and 1.5 [21].


## 3.3. Computational Meshes

The computational results of 2D model will be used to simulate the boundary conditions of 3D model, while the boundary conditions of 2D model are obtained from the observed data and its interpolation results. Taking Nandiguai-Tashiyi reach (including Hekou CR, Tiaoguan CR and Laijiapu CR) and Tashiyi-Chenglingji reach (including Jianli CR, Tianxingge CR, Fanzui CR, Xiongjiazhou CR, Qigongling CR, Guanyingzhou CR) for instances, the mesh details are as follows:(  1) Tianguan-Tashiyi Reach. Nandiguai-Tashiyi reach begins at Nandiguai and ends at Tashiyi, while the total length is 37 km. The meshed model is shown as Figure 8, and there are a total of 520×120 grid points. The type of the grid is calculated from Elliptic differential equations [22]; along the river channel there are 520 grid points and the distance between every two grid points is 40~100m; furthermore, perpendicular to the river channel there are 100 grid points and the distance between every two grid points is 10~20m.Figure 8
Mesh details of Nandiguai-Tashiyi reach.(  2) Tashiyi-Chenglingji Reach.Tashiyi-Chenglingji reach begins at Tashiyi and ends at Chenglingji, while the total length is 96 km. The meshed model is shown as Figure 9, and there are a total of 900×80 grid points. The type of the grid is the orthogonal curvilinear coordinate mesh [23]; along the river channel there are 900 grid points and the distance between every two grid points is 67~144 m; furthermore, perpendicular to the river channel there are 80 grid points and the distance between every two grid points is 5~40 m.Figure 9
Mesh details of Tashiyi-Chenglingji reach.


## 3.4. Verification

### 3.4.1. Verification of Fixed Bed Model

River toughness factor and sand transport capacity are calibrated according to the measured data, including water level, distributary ratio of branch, current velocity distribution, and river scouring and depositing amount.(  1) Water Level. Water levels at different locations are calculated by the 2D model at flood period, median water period, and drought period, respectively. With the measured data, the calculating results are shown in Tables 3, 4, and 5 accordingly. All of the water level errors are less than 0.11m. Predicted water level displays a reasonable agreement with measured data.Table 3
Water level at flood period (Q is flow rate, hm is measured value, hp is predicted value, and γh is error, γh=hp-hm).
Sectionlocation Jianli CR @ July, 2007 Fanzui CR @ Aug, 2010 Tiaoguan CR @ Sep, 2010 Q =28900 m3/s Q =20444 m3/s Q =19735 m3/s h m h c γ h h m h c γ h h m h c γ h (m) (m) (m) (m) (m) (m) (m) (m) (m) 1#L 31.46 31.5 0.04 29.02 29.00 -0.02 32.15 32.16 0.00 2#L 31.26 31.37 0.11 28.88 28.86 -0.02 31.92 31.94 0.02 3#L 31.1 31.1 0.00 28.70 28.67 -0.03 31.77 31.78 0.01 4#L 31.1 31.1 0.00 28.50 28.51 0.01 31.62 31.62 0.00 5#L 30.99 30.97 -0.02 28.35 28.40 0.05 31.46 31.43 -0.03 6#L / / / 28.15 28.16 0.01 31.33 31.29 -0.04 7#L / / / / / / 31.04 31.02 -0.02Table 4
Water level at median water period.
Sectionlocation Jianli CR @ Oct, 2008 Fanzui CR @ Sep, 2009 Q =16322 m3/s Q =14232 m3/s h m h c γ h h m h c γ h (m) (m) (m) (m) (m) (m) 1#L 29.45 29.48 0.03 27.12 27.18 0.06 2#L 29.33 29.33 0.00 26.99 27.05 0.06 3#L / / / 26.84 26.88 0.04 4#L 29.21 29.22 0.01 26.69 26.66 -0.03 5#L 28.92 28.87 -0.05 26.62 26.58 -0.04 7#L    26.42 26.37 -0.05Table 5
Water level at drought period.
Sectionlocation Jianli CR @ Jun, 2006 Qigongling CR @ Feb, 2012 Tiaoguan CR @ Feb, 2014 Q =5231 m3/s Q =6262 m3/s Q =6518 m3/s h m h c γ h h m h c γ h h m h c γ h (m) (m) (m) (m) (m) (m) (m) (m) (m) 1#L 23.11 23.17 0.06 29.02 29.00 -0.02 24.40 24.41 -0.01 2#L 22.94 22.96 0.02 28.88 28.86 -0.02 24.19 24.20 -0.01 3#L 22.57 22.52 -0.05 28.70 28.67 -0.03 24.06 24.08 -0.02 4#L 22.22 22.25 0.04 28.50 28.51 0.01 23.99 23.98 0.01 5#L 22.14 22.14 0 28.35 28.40 0.05 23.84 23.83 0.01 6#L 23.11 23.17 0.06 28.15 28.16 0.01 23.68 23.71 -0.03 7#L / / / / / / 23.51 23.56 -0.05(  2) Distributary Ratio of Branch.Jianli CR as Figure 10, a typical branching curve, is chosen to validate distributary ratio of the model. Distributary ratios in Jianli CR were calculated by the 2D model at flood period (July 2007, Q=28900 m3/s), median water period (Oct, 2008, Q=16322m3/s; Sep, 2010, Q=19175m3/s), and drought period, respectively (Jan 2006, 5231m3/s). From the results, the calculated values of each branch diversion ratio are in good agreement with the measured values and errors are within 2%. Results show that the model can reflect variation of branch diversion and meet the requirement of calculation. The split ratio verification is shown in Table 6.Table 6
Distributary ratio at different flow rate (φm is measured value, φp is predicted value, and γφ is error, γφ=φp-φm).
Q (m3/s) Date Left branch Right branch φ m(%) φ p(%) γ φ(%) φ m(%) φ p(%) γ φ(%) Low water period 5231 2006.01 3.76 3.14 -0.62 96.24 96.86 0.62  Median water period 16322 2008.10 8.50 10.36 1.86 91.50 89.64 -1.86 19175 2010.09 10.77 10.22 -0.55 89.23 89.78 0.55  High water period 28900 2007.07 8.96 8.02 -0.94 91.04 91.98 0.94Figure 10
Sketch map of Jianli CR (December, 2014).(  3) Flow Velocity Distribution Verification. Calculated flow velocity values in different curved reach are vivificated by measured values. Typical sections in Tiaoguan CR and Laijiapu CR are chosen; comparison diagram with calculated values and measured values are drew in Figure 11. All of the errors between measured value and measured values are less than 0.2m/s, and the model meets the requirements of flow velocity simulation.Figure 11
Cross section velocity verification of Tiaoguan CR on December, 2014 ( —calculated value, ■measured data): (a) Q=7633 m3/s; (b) Q=19263m3/s,.
(a) (b)


### 3.4.2. Verification of Moving Bed Model

According to the field testing points, reaches listed in Table6 are selected to compare the measured and calculated variable quantity of sediment. The geomorphic features in March, 2010 are set as the boundary conditions for the 2D numerical model; the incoming sediment and water are considered as extra loading. Table 7 shows that the 2D numerical model has a reasonable accuracy, for the maximum error is 15.2%. Therefore, the 2D numerical model can be used as a technique to study the dynamic evolution disciplines downstream of the TGD.Table 7
Comparison of depositing and scouring amount in selected reaches (from 2010.3 to 2012.2). (α is the measured data, β is the calculated value, positive values mean depositing, and minus values mean scouring. γ is the relative error, γ=β-α/α×100%).
Reach name α β γ (×104 m3) (×104 m3) (%) Tiaoguan CR 163.4 185.4 13.46 Laijiapu CR 351.9 303.8 -13.67 Tashiyi to Hongshuigang reach -549.58 -485.03 -11.7 Yanchuantao to Jingjiangmen reach -1030.35 -873.62 -15.2 Xiongjiazhou to Chenglingji reach 855.98 758.24 -11.4


## 3.4.1. Verification of Fixed Bed Model

River toughness factor and sand transport capacity are calibrated according to the measured data, including water level, distributary ratio of branch, current velocity distribution, and river scouring and depositing amount.(  1) Water Level. Water levels at different locations are calculated by the 2D model at flood period, median water period, and drought period, respectively. With the measured data, the calculating results are shown in Tables 3, 4, and 5 accordingly. All of the water level errors are less than 0.11m. Predicted water level displays a reasonable agreement with measured data.Table 3
Water level at flood period (Q is flow rate, hm is measured value, hp is predicted value, and γh is error, γh=hp-hm).
Sectionlocation Jianli CR @ July, 2007 Fanzui CR @ Aug, 2010 Tiaoguan CR @ Sep, 2010 Q =28900 m3/s Q =20444 m3/s Q =19735 m3/s h m h c γ h h m h c γ h h m h c γ h (m) (m) (m) (m) (m) (m) (m) (m) (m) 1#L 31.46 31.5 0.04 29.02 29.00 -0.02 32.15 32.16 0.00 2#L 31.26 31.37 0.11 28.88 28.86 -0.02 31.92 31.94 0.02 3#L 31.1 31.1 0.00 28.70 28.67 -0.03 31.77 31.78 0.01 4#L 31.1 31.1 0.00 28.50 28.51 0.01 31.62 31.62 0.00 5#L 30.99 30.97 -0.02 28.35 28.40 0.05 31.46 31.43 -0.03 6#L / / / 28.15 28.16 0.01 31.33 31.29 -0.04 7#L / / / / / / 31.04 31.02 -0.02Table 4
Water level at median water period.
Sectionlocation Jianli CR @ Oct, 2008 Fanzui CR @ Sep, 2009 Q =16322 m3/s Q =14232 m3/s h m h c γ h h m h c γ h (m) (m) (m) (m) (m) (m) 1#L 29.45 29.48 0.03 27.12 27.18 0.06 2#L 29.33 29.33 0.00 26.99 27.05 0.06 3#L / / / 26.84 26.88 0.04 4#L 29.21 29.22 0.01 26.69 26.66 -0.03 5#L 28.92 28.87 -0.05 26.62 26.58 -0.04 7#L    26.42 26.37 -0.05Table 5
Water level at drought period.
Sectionlocation Jianli CR @ Jun, 2006 Qigongling CR @ Feb, 2012 Tiaoguan CR @ Feb, 2014 Q =5231 m3/s Q =6262 m3/s Q =6518 m3/s h m h c γ h h m h c γ h h m h c γ h (m) (m) (m) (m) (m) (m) (m) (m) (m) 1#L 23.11 23.17 0.06 29.02 29.00 -0.02 24.40 24.41 -0.01 2#L 22.94 22.96 0.02 28.88 28.86 -0.02 24.19 24.20 -0.01 3#L 22.57 22.52 -0.05 28.70 28.67 -0.03 24.06 24.08 -0.02 4#L 22.22 22.25 0.04 28.50 28.51 0.01 23.99 23.98 0.01 5#L 22.14 22.14 0 28.35 28.40 0.05 23.84 23.83 0.01 6#L 23.11 23.17 0.06 28.15 28.16 0.01 23.68 23.71 -0.03 7#L / / / / / / 23.51 23.56 -0.05(  2) Distributary Ratio of Branch.Jianli CR as Figure 10, a typical branching curve, is chosen to validate distributary ratio of the model. Distributary ratios in Jianli CR were calculated by the 2D model at flood period (July 2007, Q=28900 m3/s), median water period (Oct, 2008, Q=16322m3/s; Sep, 2010, Q=19175m3/s), and drought period, respectively (Jan 2006, 5231m3/s). From the results, the calculated values of each branch diversion ratio are in good agreement with the measured values and errors are within 2%. Results show that the model can reflect variation of branch diversion and meet the requirement of calculation. The split ratio verification is shown in Table 6.Table 6
Distributary ratio at different flow rate (φm is measured value, φp is predicted value, and γφ is error, γφ=φp-φm).
Q (m3/s) Date Left branch Right branch φ m(%) φ p(%) γ φ(%) φ m(%) φ p(%) γ φ(%) Low water period 5231 2006.01 3.76 3.14 -0.62 96.24 96.86 0.62  Median water period 16322 2008.10 8.50 10.36 1.86 91.50 89.64 -1.86 19175 2010.09 10.77 10.22 -0.55 89.23 89.78 0.55  High water period 28900 2007.07 8.96 8.02 -0.94 91.04 91.98 0.94Figure 10
Sketch map of Jianli CR (December, 2014).(  3) Flow Velocity Distribution Verification. Calculated flow velocity values in different curved reach are vivificated by measured values. Typical sections in Tiaoguan CR and Laijiapu CR are chosen; comparison diagram with calculated values and measured values are drew in Figure 11. All of the errors between measured value and measured values are less than 0.2m/s, and the model meets the requirements of flow velocity simulation.Figure 11
Cross section velocity verification of Tiaoguan CR on December, 2014 ( —calculated value, ■measured data): (a) Q=7633 m3/s; (b) Q=19263m3/s,.
(a) (b)


## 3.4.2. Verification of Moving Bed Model

According to the field testing points, reaches listed in Table6 are selected to compare the measured and calculated variable quantity of sediment. The geomorphic features in March, 2010 are set as the boundary conditions for the 2D numerical model; the incoming sediment and water are considered as extra loading. Table 7 shows that the 2D numerical model has a reasonable accuracy, for the maximum error is 15.2%. Therefore, the 2D numerical model can be used as a technique to study the dynamic evolution disciplines downstream of the TGD.Table 7
Comparison of depositing and scouring amount in selected reaches (from 2010.3 to 2012.2). (α is the measured data, β is the calculated value, positive values mean depositing, and minus values mean scouring. γ is the relative error, γ=β-α/α×100%).
Reach name α β γ (×104 m3) (×104 m3) (%) Tiaoguan CR 163.4 185.4 13.46 Laijiapu CR 351.9 303.8 -13.67 Tashiyi to Hongshuigang reach -549.58 -485.03 -11.7 Yanchuantao to Jingjiangmen reach -1030.35 -873.62 -15.2 Xiongjiazhou to Chenglingji reach 855.98 758.24 -11.4


## 4. Three-Dimensional Model

The 2D model could describe the plant evolution disciplines of curved reaches. For curved reaches with strong three-dimensional characteristics, 2D model could not describe the vertical characteristic of water and sediment. Thus, a superior three-dimensional (3D) numerical model was built to analyze the water and sediment movement characteristic. The 3D model’s boundary conditions were given by the results of the 2D numerical model. The two models were verified by the measured data, and part of results by 3D model was also verified by 2D model.


### 4.1. Fundamental Equations

#### 4.1.1. Water Flow Equation

The three-dimensional (3D) water flow equation in the Cartesian coordinate system is (5), and Navier-Stokes equation [24] is (6):(5)∂ui∂xi=0(6)ρdudt=-∂P∂xi+PXi+μΔuiwhere xi is the coordinate axis of the system; ρ is fluid density; P is water pressure; ui is water current velocity; μ is water dynamic viscosity coefficient (generally equal to 1×10-6m2sec-1).For the Large Eddy Simulation (LES) model [25], the physical variable should be factorized into large scale pulsation and small scale pulsation after the fundamental equation done the space-filtering operation. For example,f is a 3D physical variable; differential operator off can be space-filtering operated as (7)∂fx¯∂x=∂fx¯∂x(8)∂nfx¯∂xn=∂nfx¯∂xnEquations (7) and (8), Smagorinsky pattern [26, 27], and Space filtering operation (9) are used to solve (5) and (6); therefore (10a)–(10g) were achieved:(9)f¯¯≠f¯,f¯g¯¯≠f¯g¯,f¯g′¯≠0(10a)∂ui¯∂xi=0(10b)∂ui¯∂t=-1ρ∂P∂xi-∂ui¯uj¯∂xj+∂∂xj2νeSij¯+gi(10c)νe=ν+νt(10d)νt=CsΔ22Sij¯Sij¯1/2(10e)Sij¯=12∂ui¯∂xj+∂uj¯∂xi(10f)Δ=Δx1Δx2Δx31/3(10g)P=p¯+23qIn the equations, the upper line represents the space filter variable;νt is spectral space eddy viscosity coefficient; ν is motion viscosity coefficient, taking 1.0∗10−6 m2sec−1;Cs is Smagorinsky coefficient, taking 0.1; Δxi is grid width;p is pressure;q is spectral space eddy kinetic energy.


#### 4.1.2. Unbalanced Sediment Transport Equation

The unbalanced sediment transport equation can be deduced based on the diffusion theory and displayed as(11)∂S∂t+∂∂xiui-ωsδi3S=∂∂xiεs∂S∂xiwhere S is sediment concentration; δi3 is the Kronecker factor, (12)δi3=1,i=30,i≠3;ωs is deposition velocity of sediment; σs is the Schmidt factor, generally 0.5≤σs≤1.0. However, it [28] suggested σs=1.0 in their 3D numerical model.


#### 4.1.3. Suspended Sediment Transportation Equation

(13) ω s S + ε s ∂ S ∂ z = D b - E b = γ s ′ ∂ Z b ∂ t


#### 4.1.4. Bed Load Transport Equation

(  1) Bed Load Discharge.The simulation of bed load discharge is one of the pivotal hypotheses of the 3D numerical model. The bulk bed load discharge qb∗ is a function of the particle velocity Ub, the granule jumps height δb, and the bulk sediment concentration Sbv; the relationship shown as (14) and [30] (as (15)) has been one the most popularly used equations:(14)
q
b
∗
=
U
b
δ
b
S
b
v
(15)
q
b
∗
=
0.053
ρ
s
-
ρ
ρ
g
D
50
1.5
T
2.1
D
∗
0.3
T
<
3
0.100
ρ
s
-
ρ
ρ
g
D
50
1.5
T
1.5
D
∗
0.3
T
≥
3where D∗=D50ρs-ρ/ρg/ν21/3 is particle parameter and D50 is the corresponding particle size with a cumulative percentage distribution of 50%; T=τb′-τb,cr/τb,cr is the variable of phase shifting, τb′ is the effective shear stress of river bed, τb,cr is the critical effective shear stress of river bed, and τb′=abτb,cr, αb=(C/C′)2,C is the Chezy coefficient [31], C′=18log⁡(12Rb/3D90); Rb is the hydraulic radius and D90 is the corresponding particle size with a cumulative percentage distribution of 90%;Afterwards, after measured four sectional bed load discharge and suspended sediment discharge of the Nile, along with a series filed testing date, Abdel-Fattah [32] proposed a new factor λ to modify the bed load discharge. Hence, (14) could be rewritten as (16)
q
b
∗
=
0.053
ρ
s
-
ρ
ρ
g
D
50
1.5
T
2.1
D
∗
0.3
T
<
2.5
0.100
ρ
s
-
ρ
ρ
g
D
50
1.5
T
m
1.5
D
∗
0.3
T
>
2.5where Tm=λτb′-τb,cr/τb,cr, λ could be calculated by (17)λ=e0.45a+0.2βorλ=e1.8-2.4σg+0.6σg2where α=1-σg, β=d50/d90-d10/d50, and σg=d84/d50+d50/d16/2.The critical effective shear stress of river bedτb,cr can be deduced by the flat river bed experiment [33], shown as (18)τb,cr=ρs-ρgθcrD50Rijin [34] had modified the following by a series research, considering the longitudinal slope of river bed effects on the Shield critical shear stress:(19)τb,cr,s=kβkγτb,crwhere kβ is longitudinal slope factor, kβ=sin⁡(φ-β)/sin⁡φ, β is the longitudinal slope toe, it is positive when downward grade while it is negative when upward grade, and φ is slope of repose. kγ is horizontal slope factor, kγ=cos⁡γ(1-tan2γ/tan2φ), and γ is the horizontal slope toe.Furthermore, bed load discharge should be revised as the influence of the slope. It could be divided into longitudinal bed load sediment rate and transversal bed load sediment rate.Formula of modified longitudinal bed load sediment rate after being revised is(20)qb,s∗′=αsqb,s∗Transversal bed load sediment rate formula after being revised is(21)qb,n∗′=vbub+ετb,crτb,s0.5tan⁡γqb,s∗′where qb,s∗=ub/Ubqb∗ is transversal bed load sediment rate formula; qb,n∗=vb/Ubqb∗ is transversal bed load sediment rate formula; Ub was confluence velocity of transverse velocity (vb) and longitudinal velocity (ub); αs is the slope surface influence coefficient; τb,s is the longitudinal component of the shear stress for the bed surface; ε is a parameter (generally equal to 1.5).(  2) Unbalanced Sediment Transport Equation of Bed Load.The unbalanced sediment transport equation of bed load was proposed by Rijn [29]:(22)γs′ρs∂zb,b∂t+∂δbSbv′∂t+Db-Eb+∂qb,s∂x+∂qb,n∂y=0where Sbv′ is the bulk sediment concentration of the bottom sand bed; ∂Zb,b/∂t is the bed deformation caused by the bed load transportation.Wellington [35] had supposed that(23)γs′ρs∂zb,b∂t=1Lsqb-qb∗where Ls is the recovery distance of the bottom sand bed, Ls=3d50D∗0.5T0.9.When the effects of suspended sediment are very small and can be ignored, (22) can be rewritten as (24)1Lsqb-qb∗+∂δbSbv′∂t+∂qb,s∂x+∂qb,n∂y=0Becauseqb=UbδbSbv′, (24) can be written as(25)1Lsqb-qb∗+∂qb/Ub∂t+∂qb,s∂x+∂qb,n∂y=0Thereby, if the time derivative terms of suspended sediment are ignored, the deformation equation of bed load is(26)1Lsqb-qb∗+∂qb,s∂x+∂qb,n∂y=0Deformation of river bed is the sum of deformations caused by suspended sediment and bed load sediment.


### 4.2. Initial Conditions

Physical variables such as current velocity, water level, and turbulent kinetic energy are aggravated 3% as pulsating part from the original measured date, respectively. Initial value of velocity and water level of every mesh had be calculated by 2D model and average allocated to every mesh in 3D model. k coefficient andε coefficient value could be given by experience. Dynamic water pressure is zero.


### 4.3. Boundary Conditions

#### 4.3.1. Entrance Boundary

Suspended sediment concentration distribution in river entrance section used measured date, or vertical sediment concentration distribution formula. Bed load sediment concentration in entrance section is given to keep no deformation occurring in the section.


#### 4.3.2. Export Boundary

The normal gradient of suspended sediment concentration in export section is zero.


#### 4.3.3. Wall Surface and Bottom Surface Boundaries

(  1) Current Velocity at the Surfaces.The walls and bottom of the river are coarse surfaces; they would hinder the water flow away. In this 3D model, the “law of the wall” (Ban, 1995) is used to simulate the current velocity at the interfaces, which is the water and solid interface.(27)uu∗=1kln⁡zu∗ν+Bwhere u∗ is friction velocity; k=0.41; B=5.0.At the meanwhile, the shear stress at the water up surface is supposed to be zero, and the shear stress at the water and soil interface can be deduced from(28)τwρ=ν∂u∂z(  2) Stress at the Water-Soil Interface.The stress at the water-soil interface can be calculated by (29)τbx=ρCdubub2+vb2τby=ρCdvbub2+vb2In the moving bed model, current velocity can be calculated by (29) use logarithmic wall function:(30)ubu∗=1kln⁡zz0where friction velocity u∗=τb/ρ. Thus,Cd can be calculated by (31)~(34):(31)Cd=11/kln⁡zb/z02(32)z0=νEu∗(33)E=eκB-ΔB(34)ΔB=0Re<2.25B-8.5+1κln⁡Resin⁡0.428+ln⁡Re-0.8112.25≤Re<90B-8.5+1κln⁡ReRe≥90(35)Re=u∗ksνwhere z0 is the height when current velocity is equal to 0; E is rough parameters of bed surface; B=5.2; ks is rough height of river bottom or riverbank, ks=2.5d [30].(  3) Sediment Concentration at the Surfaces.These is no sediment exchange at the surfaces, and sediment concentration at the surfaces can be deduced from (36)ωsS+εs∂S∂zz=zfree=0(  4) Sediment Concentration at the Near-Bed Region.The near-bed region is the layer in the water but right near the river bed, the sediment concentration of the near-bed region Sb is a significant hypothesis for the 3D model, and it can be calculated by (37)S=Sb-Sb∗+ce-ω/εszAnd replacing the variables with the known equations: z=δb,S=Sb,z=zk-δb,and S=Sk, 40 can be rewritten as follows:(38)Sb=Sk+Sb∗1-e-ω/εszk-δb(  5) Sediment Carrying Capacity of the Near-Bed Region.As a more commonly used method at present, van Rijn is used to calculate sediment carrying capacity. Sediment carrying capacity of the near-bed region can be deduced from (39)Svb∗=0.015d50T1.5aD∗0.3(40)D∗=d50ρs-ρ/ρgν21/3(41)a=0.5ks(42)a≥0.01d50


### 4.4. Computational Meshes

Unstructured grid mesh method is used in the proposed 3D model, and three typical meandering reaches have been chosen to verify the mesh technique:(  1) Tiaoguan CR. Tiaoguan CR is about 11 km in length and 1.1~2km in width, the mesh size is 20~30 m in the horizontal direction, and along the water channel depth direction there are 10 layers grids, and there are a total of 5925 grids, shown as Figure 12(a).Figure 12
Mesh details of: (a) Tiaoguan CR, (b) Laijiapu CR, and (c) Qigongling CR.
(a) (b) (c)(  2) Laijiapu CR.Laijiapu CR is about 8 km in length and 0.7~1.3km in width, the mesh size is 17~30 m in the horizontal direction, and along the water channel depth direction there are 10 layers grids, and there are a total of 3604 grids, shown as Figure 12(b).(  3) Chibakou CR.Chibakou CR is about 9 km in length and 0.7~1.7km in width, the mesh size is 20~50 m in the horizontal direction, and along the water channel depth direction there are 10 layers grids, and there are a total of 4905 grids, shown as Figure 12(c).


### 4.5. Verification

#### 4.5.1. Verification by Flume Test

Onishi classic model test [36] was used to validate the mathematical model. The given water level at the model entrance was 0.1333 m, the average current velocity was 0.536 m/s, and the average diameter of the particles was 0.25 mm, shown as Figure 13. Verification results show that the average water lever errors between calculated value and testing value by total sediment transport model (TSTM) are less than 0.06 m in all sections, while errors are less than 0.04 m by suspended sediment transport model (SSTM). Thus, TSTM is chosen to simulate the river, shown as Figure 14. Vertical flow velocity errors are mainly produced by fluctuation of water level. The errors are less than 0.2m/s, shown as Figure 15. The predicted sediment concentration in the top layer was almost the same as the testing data, while the predicted sediment concentration of in bottom layer was vibrating within 100 g/m3 at the testing data, shown as Figure 16.Figure 13
Onishi classic model test plane layout.Figure 14
Water depth value contrast. ○: the first test set; △: the second test set; —: calculated value by TSTM. – –: calculated value by SSTM.Figure 15
Vertical flow velocity contrast. —: calculated value. ○: test value.Figure 16
Sediment concentration value contrast. —: calculated value. □: test value.Briefly, the parameters calibration results indicate that the proposed 3D model displayed a reasonable agreement with the testing data, and could be used to do comprehensive verification for fixed bed model and moving bed model.


#### 4.5.2. Verification by Measured Data

The fixed bed model is vivificated by measured data of February 2010. Figures17 and 18 show the velocity comparison of cross section and longitudinal section in chosen reaches. Velocity errors of cross section are within 0.2m/s. Velocity errors of longitudinal section are within 0.4m/s.Figure 17
Velocity contrast of cross section.Figure 18
Velocity contrast of longitudinal section: (a) Tiaoguan CR, (b) Laijiapu CR, and (c) Qigongling CR.
(a) (b) (c)The moving bed model used the geomorphic data of February 2012 as the original inputs and used the geomorphic data of October 2012 as the verification materials. Figure19 shows that the predicted scouring or depositing tendency of the reaches is the same as the measured data and the predicted variation of the sediment is slightly vibrating within 20% error near the measured data, which already meet the calculation accuracy of engineering projects.Figure 19
Scouring or depositing distribution diagram from February, 2012 to October, 2012. (a) Tiaoguan CR, calculated value. (b) Tiaoguan CR, measured data. (c) Laijiapu CR, calculated value. (d) Laijiapu CR, measured data. (e) Qigongling CR, calculated value. (f) Qigongling CR, measured data.
(a) (b) (c) (d) (e) (f)According to the verification of moving bed model, the proposed 3D model behaves reliably in predict the transportation properties of water, sediment, and even the deformation of the river bed. It can be used in further study on river evolution mechanism and can be used in broader prediction of river evolution disciplines.


## 4.1. Fundamental Equations

### 4.1.1. Water Flow Equation

The three-dimensional (3D) water flow equation in the Cartesian coordinate system is (5), and Navier-Stokes equation [24] is (6):(5)∂ui∂xi=0(6)ρdudt=-∂P∂xi+PXi+μΔuiwhere xi is the coordinate axis of the system; ρ is fluid density; P is water pressure; ui is water current velocity; μ is water dynamic viscosity coefficient (generally equal to 1×10-6m2sec-1).For the Large Eddy Simulation (LES) model [25], the physical variable should be factorized into large scale pulsation and small scale pulsation after the fundamental equation done the space-filtering operation. For example,f is a 3D physical variable; differential operator off can be space-filtering operated as (7)∂fx¯∂x=∂fx¯∂x(8)∂nfx¯∂xn=∂nfx¯∂xnEquations (7) and (8), Smagorinsky pattern [26, 27], and Space filtering operation (9) are used to solve (5) and (6); therefore (10a)–(10g) were achieved:(9)f¯¯≠f¯,f¯g¯¯≠f¯g¯,f¯g′¯≠0(10a)∂ui¯∂xi=0(10b)∂ui¯∂t=-1ρ∂P∂xi-∂ui¯uj¯∂xj+∂∂xj2νeSij¯+gi(10c)νe=ν+νt(10d)νt=CsΔ22Sij¯Sij¯1/2(10e)Sij¯=12∂ui¯∂xj+∂uj¯∂xi(10f)Δ=Δx1Δx2Δx31/3(10g)P=p¯+23qIn the equations, the upper line represents the space filter variable;νt is spectral space eddy viscosity coefficient; ν is motion viscosity coefficient, taking 1.0∗10−6 m2sec−1;Cs is Smagorinsky coefficient, taking 0.1; Δxi is grid width;p is pressure;q is spectral space eddy kinetic energy.


### 4.1.2. Unbalanced Sediment Transport Equation

The unbalanced sediment transport equation can be deduced based on the diffusion theory and displayed as(11)∂S∂t+∂∂xiui-ωsδi3S=∂∂xiεs∂S∂xiwhere S is sediment concentration; δi3 is the Kronecker factor, (12)δi3=1,i=30,i≠3;ωs is deposition velocity of sediment; σs is the Schmidt factor, generally 0.5≤σs≤1.0. However, it [28] suggested σs=1.0 in their 3D numerical model.


### 4.1.3. Suspended Sediment Transportation Equation

(13) ω s S + ε s ∂ S ∂ z = D b - E b = γ s ′ ∂ Z b ∂ t


### 4.1.4. Bed Load Transport Equation

(  1) Bed Load Discharge.The simulation of bed load discharge is one of the pivotal hypotheses of the 3D numerical model. The bulk bed load discharge qb∗ is a function of the particle velocity Ub, the granule jumps height δb, and the bulk sediment concentration Sbv; the relationship shown as (14) and [30] (as (15)) has been one the most popularly used equations:(14)
q
b
∗
=
U
b
δ
b
S
b
v
(15)
q
b
∗
=
0.053
ρ
s
-
ρ
ρ
g
D
50
1.5
T
2.1
D
∗
0.3
T
<
3
0.100
ρ
s
-
ρ
ρ
g
D
50
1.5
T
1.5
D
∗
0.3
T
≥
3where D∗=D50ρs-ρ/ρg/ν21/3 is particle parameter and D50 is the corresponding particle size with a cumulative percentage distribution of 50%; T=τb′-τb,cr/τb,cr is the variable of phase shifting, τb′ is the effective shear stress of river bed, τb,cr is the critical effective shear stress of river bed, and τb′=abτb,cr, αb=(C/C′)2,C is the Chezy coefficient [31], C′=18log⁡(12Rb/3D90); Rb is the hydraulic radius and D90 is the corresponding particle size with a cumulative percentage distribution of 90%;Afterwards, after measured four sectional bed load discharge and suspended sediment discharge of the Nile, along with a series filed testing date, Abdel-Fattah [32] proposed a new factor λ to modify the bed load discharge. Hence, (14) could be rewritten as (16)
q
b
∗
=
0.053
ρ
s
-
ρ
ρ
g
D
50
1.5
T
2.1
D
∗
0.3
T
<
2.5
0.100
ρ
s
-
ρ
ρ
g
D
50
1.5
T
m
1.5
D
∗
0.3
T
>
2.5where Tm=λτb′-τb,cr/τb,cr, λ could be calculated by (17)λ=e0.45a+0.2βorλ=e1.8-2.4σg+0.6σg2where α=1-σg, β=d50/d90-d10/d50, and σg=d84/d50+d50/d16/2.The critical effective shear stress of river bedτb,cr can be deduced by the flat river bed experiment [33], shown as (18)τb,cr=ρs-ρgθcrD50Rijin [34] had modified the following by a series research, considering the longitudinal slope of river bed effects on the Shield critical shear stress:(19)τb,cr,s=kβkγτb,crwhere kβ is longitudinal slope factor, kβ=sin⁡(φ-β)/sin⁡φ, β is the longitudinal slope toe, it is positive when downward grade while it is negative when upward grade, and φ is slope of repose. kγ is horizontal slope factor, kγ=cos⁡γ(1-tan2γ/tan2φ), and γ is the horizontal slope toe.Furthermore, bed load discharge should be revised as the influence of the slope. It could be divided into longitudinal bed load sediment rate and transversal bed load sediment rate.Formula of modified longitudinal bed load sediment rate after being revised is(20)qb,s∗′=αsqb,s∗Transversal bed load sediment rate formula after being revised is(21)qb,n∗′=vbub+ετb,crτb,s0.5tan⁡γqb,s∗′where qb,s∗=ub/Ubqb∗ is transversal bed load sediment rate formula; qb,n∗=vb/Ubqb∗ is transversal bed load sediment rate formula; Ub was confluence velocity of transverse velocity (vb) and longitudinal velocity (ub); αs is the slope surface influence coefficient; τb,s is the longitudinal component of the shear stress for the bed surface; ε is a parameter (generally equal to 1.5).(  2) Unbalanced Sediment Transport Equation of Bed Load.The unbalanced sediment transport equation of bed load was proposed by Rijn [29]:(22)γs′ρs∂zb,b∂t+∂δbSbv′∂t+Db-Eb+∂qb,s∂x+∂qb,n∂y=0where Sbv′ is the bulk sediment concentration of the bottom sand bed; ∂Zb,b/∂t is the bed deformation caused by the bed load transportation.Wellington [35] had supposed that(23)γs′ρs∂zb,b∂t=1Lsqb-qb∗where Ls is the recovery distance of the bottom sand bed, Ls=3d50D∗0.5T0.9.When the effects of suspended sediment are very small and can be ignored, (22) can be rewritten as (24)1Lsqb-qb∗+∂δbSbv′∂t+∂qb,s∂x+∂qb,n∂y=0Becauseqb=UbδbSbv′, (24) can be written as(25)1Lsqb-qb∗+∂qb/Ub∂t+∂qb,s∂x+∂qb,n∂y=0Thereby, if the time derivative terms of suspended sediment are ignored, the deformation equation of bed load is(26)1Lsqb-qb∗+∂qb,s∂x+∂qb,n∂y=0Deformation of river bed is the sum of deformations caused by suspended sediment and bed load sediment.


## 4.1.1. Water Flow Equation

The three-dimensional (3D) water flow equation in the Cartesian coordinate system is (5), and Navier-Stokes equation [24] is (6):(5)∂ui∂xi=0(6)ρdudt=-∂P∂xi+PXi+μΔuiwhere xi is the coordinate axis of the system; ρ is fluid density; P is water pressure; ui is water current velocity; μ is water dynamic viscosity coefficient (generally equal to 1×10-6m2sec-1).For the Large Eddy Simulation (LES) model [25], the physical variable should be factorized into large scale pulsation and small scale pulsation after the fundamental equation done the space-filtering operation. For example,f is a 3D physical variable; differential operator off can be space-filtering operated as (7)∂fx¯∂x=∂fx¯∂x(8)∂nfx¯∂xn=∂nfx¯∂xnEquations (7) and (8), Smagorinsky pattern [26, 27], and Space filtering operation (9) are used to solve (5) and (6); therefore (10a)–(10g) were achieved:(9)f¯¯≠f¯,f¯g¯¯≠f¯g¯,f¯g′¯≠0(10a)∂ui¯∂xi=0(10b)∂ui¯∂t=-1ρ∂P∂xi-∂ui¯uj¯∂xj+∂∂xj2νeSij¯+gi(10c)νe=ν+νt(10d)νt=CsΔ22Sij¯Sij¯1/2(10e)Sij¯=12∂ui¯∂xj+∂uj¯∂xi(10f)Δ=Δx1Δx2Δx31/3(10g)P=p¯+23qIn the equations, the upper line represents the space filter variable;νt is spectral space eddy viscosity coefficient; ν is motion viscosity coefficient, taking 1.0∗10−6 m2sec−1;Cs is Smagorinsky coefficient, taking 0.1; Δxi is grid width;p is pressure;q is spectral space eddy kinetic energy.


## 4.1.2. Unbalanced Sediment Transport Equation

The unbalanced sediment transport equation can be deduced based on the diffusion theory and displayed as(11)∂S∂t+∂∂xiui-ωsδi3S=∂∂xiεs∂S∂xiwhere S is sediment concentration; δi3 is the Kronecker factor, (12)δi3=1,i=30,i≠3;ωs is deposition velocity of sediment; σs is the Schmidt factor, generally 0.5≤σs≤1.0. However, it [28] suggested σs=1.0 in their 3D numerical model.


## 4.1.3. Suspended Sediment Transportation Equation

(13) ω s S + ε s ∂ S ∂ z = D b - E b = γ s ′ ∂ Z b ∂ t


## 4.1.4. Bed Load Transport Equation

(  1) Bed Load Discharge.The simulation of bed load discharge is one of the pivotal hypotheses of the 3D numerical model. The bulk bed load discharge qb∗ is a function of the particle velocity Ub, the granule jumps height δb, and the bulk sediment concentration Sbv; the relationship shown as (14) and [30] (as (15)) has been one the most popularly used equations:(14)
q
b
∗
=
U
b
δ
b
S
b
v
(15)
q
b
∗
=
0.053
ρ
s
-
ρ
ρ
g
D
50
1.5
T
2.1
D
∗
0.3
T
<
3
0.100
ρ
s
-
ρ
ρ
g
D
50
1.5
T
1.5
D
∗
0.3
T
≥
3where D∗=D50ρs-ρ/ρg/ν21/3 is particle parameter and D50 is the corresponding particle size with a cumulative percentage distribution of 50%; T=τb′-τb,cr/τb,cr is the variable of phase shifting, τb′ is the effective shear stress of river bed, τb,cr is the critical effective shear stress of river bed, and τb′=abτb,cr, αb=(C/C′)2,C is the Chezy coefficient [31], C′=18log⁡(12Rb/3D90); Rb is the hydraulic radius and D90 is the corresponding particle size with a cumulative percentage distribution of 90%;Afterwards, after measured four sectional bed load discharge and suspended sediment discharge of the Nile, along with a series filed testing date, Abdel-Fattah [32] proposed a new factor λ to modify the bed load discharge. Hence, (14) could be rewritten as (16)
q
b
∗
=
0.053
ρ
s
-
ρ
ρ
g
D
50
1.5
T
2.1
D
∗
0.3
T
<
2.5
0.100
ρ
s
-
ρ
ρ
g
D
50
1.5
T
m
1.5
D
∗
0.3
T
>
2.5where Tm=λτb′-τb,cr/τb,cr, λ could be calculated by (17)λ=e0.45a+0.2βorλ=e1.8-2.4σg+0.6σg2where α=1-σg, β=d50/d90-d10/d50, and σg=d84/d50+d50/d16/2.The critical effective shear stress of river bedτb,cr can be deduced by the flat river bed experiment [33], shown as (18)τb,cr=ρs-ρgθcrD50Rijin [34] had modified the following by a series research, considering the longitudinal slope of river bed effects on the Shield critical shear stress:(19)τb,cr,s=kβkγτb,crwhere kβ is longitudinal slope factor, kβ=sin⁡(φ-β)/sin⁡φ, β is the longitudinal slope toe, it is positive when downward grade while it is negative when upward grade, and φ is slope of repose. kγ is horizontal slope factor, kγ=cos⁡γ(1-tan2γ/tan2φ), and γ is the horizontal slope toe.Furthermore, bed load discharge should be revised as the influence of the slope. It could be divided into longitudinal bed load sediment rate and transversal bed load sediment rate.Formula of modified longitudinal bed load sediment rate after being revised is(20)qb,s∗′=αsqb,s∗Transversal bed load sediment rate formula after being revised is(21)qb,n∗′=vbub+ετb,crτb,s0.5tan⁡γqb,s∗′where qb,s∗=ub/Ubqb∗ is transversal bed load sediment rate formula; qb,n∗=vb/Ubqb∗ is transversal bed load sediment rate formula; Ub was confluence velocity of transverse velocity (vb) and longitudinal velocity (ub); αs is the slope surface influence coefficient; τb,s is the longitudinal component of the shear stress for the bed surface; ε is a parameter (generally equal to 1.5).(  2) Unbalanced Sediment Transport Equation of Bed Load.The unbalanced sediment transport equation of bed load was proposed by Rijn [29]:(22)γs′ρs∂zb,b∂t+∂δbSbv′∂t+Db-Eb+∂qb,s∂x+∂qb,n∂y=0where Sbv′ is the bulk sediment concentration of the bottom sand bed; ∂Zb,b/∂t is the bed deformation caused by the bed load transportation.Wellington [35] had supposed that(23)γs′ρs∂zb,b∂t=1Lsqb-qb∗where Ls is the recovery distance of the bottom sand bed, Ls=3d50D∗0.5T0.9.When the effects of suspended sediment are very small and can be ignored, (22) can be rewritten as (24)1Lsqb-qb∗+∂δbSbv′∂t+∂qb,s∂x+∂qb,n∂y=0Becauseqb=UbδbSbv′, (24) can be written as(25)1Lsqb-qb∗+∂qb/Ub∂t+∂qb,s∂x+∂qb,n∂y=0Thereby, if the time derivative terms of suspended sediment are ignored, the deformation equation of bed load is(26)1Lsqb-qb∗+∂qb,s∂x+∂qb,n∂y=0Deformation of river bed is the sum of deformations caused by suspended sediment and bed load sediment.


## 4.2. Initial Conditions

Physical variables such as current velocity, water level, and turbulent kinetic energy are aggravated 3% as pulsating part from the original measured date, respectively. Initial value of velocity and water level of every mesh had be calculated by 2D model and average allocated to every mesh in 3D model. k coefficient andε coefficient value could be given by experience. Dynamic water pressure is zero.


## 4.3. Boundary Conditions

### 4.3.1. Entrance Boundary

Suspended sediment concentration distribution in river entrance section used measured date, or vertical sediment concentration distribution formula. Bed load sediment concentration in entrance section is given to keep no deformation occurring in the section.


### 4.3.2. Export Boundary

The normal gradient of suspended sediment concentration in export section is zero.


### 4.3.3. Wall Surface and Bottom Surface Boundaries

(  1) Current Velocity at the Surfaces.The walls and bottom of the river are coarse surfaces; they would hinder the water flow away. In this 3D model, the “law of the wall” (Ban, 1995) is used to simulate the current velocity at the interfaces, which is the water and solid interface.(27)uu∗=1kln⁡zu∗ν+Bwhere u∗ is friction velocity; k=0.41; B=5.0.At the meanwhile, the shear stress at the water up surface is supposed to be zero, and the shear stress at the water and soil interface can be deduced from(28)τwρ=ν∂u∂z(  2) Stress at the Water-Soil Interface.The stress at the water-soil interface can be calculated by (29)τbx=ρCdubub2+vb2τby=ρCdvbub2+vb2In the moving bed model, current velocity can be calculated by (29) use logarithmic wall function:(30)ubu∗=1kln⁡zz0where friction velocity u∗=τb/ρ. Thus,Cd can be calculated by (31)~(34):(31)Cd=11/kln⁡zb/z02(32)z0=νEu∗(33)E=eκB-ΔB(34)ΔB=0Re<2.25B-8.5+1κln⁡Resin⁡0.428+ln⁡Re-0.8112.25≤Re<90B-8.5+1κln⁡ReRe≥90(35)Re=u∗ksνwhere z0 is the height when current velocity is equal to 0; E is rough parameters of bed surface; B=5.2; ks is rough height of river bottom or riverbank, ks=2.5d [30].(  3) Sediment Concentration at the Surfaces.These is no sediment exchange at the surfaces, and sediment concentration at the surfaces can be deduced from (36)ωsS+εs∂S∂zz=zfree=0(  4) Sediment Concentration at the Near-Bed Region.The near-bed region is the layer in the water but right near the river bed, the sediment concentration of the near-bed region Sb is a significant hypothesis for the 3D model, and it can be calculated by (37)S=Sb-Sb∗+ce-ω/εszAnd replacing the variables with the known equations: z=δb,S=Sb,z=zk-δb,and S=Sk, 40 can be rewritten as follows:(38)Sb=Sk+Sb∗1-e-ω/εszk-δb(  5) Sediment Carrying Capacity of the Near-Bed Region.As a more commonly used method at present, van Rijn is used to calculate sediment carrying capacity. Sediment carrying capacity of the near-bed region can be deduced from (39)Svb∗=0.015d50T1.5aD∗0.3(40)D∗=d50ρs-ρ/ρgν21/3(41)a=0.5ks(42)a≥0.01d50


## 4.3.1. Entrance Boundary

Suspended sediment concentration distribution in river entrance section used measured date, or vertical sediment concentration distribution formula. Bed load sediment concentration in entrance section is given to keep no deformation occurring in the section.


## 4.3.2. Export Boundary

The normal gradient of suspended sediment concentration in export section is zero.


## 4.3.3. Wall Surface and Bottom Surface Boundaries

(  1) Current Velocity at the Surfaces.The walls and bottom of the river are coarse surfaces; they would hinder the water flow away. In this 3D model, the “law of the wall” (Ban, 1995) is used to simulate the current velocity at the interfaces, which is the water and solid interface.(27)uu∗=1kln⁡zu∗ν+Bwhere u∗ is friction velocity; k=0.41; B=5.0.At the meanwhile, the shear stress at the water up surface is supposed to be zero, and the shear stress at the water and soil interface can be deduced from(28)τwρ=ν∂u∂z(  2) Stress at the Water-Soil Interface.The stress at the water-soil interface can be calculated by (29)τbx=ρCdubub2+vb2τby=ρCdvbub2+vb2In the moving bed model, current velocity can be calculated by (29) use logarithmic wall function:(30)ubu∗=1kln⁡zz0where friction velocity u∗=τb/ρ. Thus,Cd can be calculated by (31)~(34):(31)Cd=11/kln⁡zb/z02(32)z0=νEu∗(33)E=eκB-ΔB(34)ΔB=0Re<2.25B-8.5+1κln⁡Resin⁡0.428+ln⁡Re-0.8112.25≤Re<90B-8.5+1κln⁡ReRe≥90(35)Re=u∗ksνwhere z0 is the height when current velocity is equal to 0; E is rough parameters of bed surface; B=5.2; ks is rough height of river bottom or riverbank, ks=2.5d [30].(  3) Sediment Concentration at the Surfaces.These is no sediment exchange at the surfaces, and sediment concentration at the surfaces can be deduced from (36)ωsS+εs∂S∂zz=zfree=0(  4) Sediment Concentration at the Near-Bed Region.The near-bed region is the layer in the water but right near the river bed, the sediment concentration of the near-bed region Sb is a significant hypothesis for the 3D model, and it can be calculated by (37)S=Sb-Sb∗+ce-ω/εszAnd replacing the variables with the known equations: z=δb,S=Sb,z=zk-δb,and S=Sk, 40 can be rewritten as follows:(38)Sb=Sk+Sb∗1-e-ω/εszk-δb(  5) Sediment Carrying Capacity of the Near-Bed Region.As a more commonly used method at present, van Rijn is used to calculate sediment carrying capacity. Sediment carrying capacity of the near-bed region can be deduced from (39)Svb∗=0.015d50T1.5aD∗0.3(40)D∗=d50ρs-ρ/ρgν21/3(41)a=0.5ks(42)a≥0.01d50


## 4.4. Computational Meshes

Unstructured grid mesh method is used in the proposed 3D model, and three typical meandering reaches have been chosen to verify the mesh technique:(  1) Tiaoguan CR. Tiaoguan CR is about 11 km in length and 1.1~2km in width, the mesh size is 20~30 m in the horizontal direction, and along the water channel depth direction there are 10 layers grids, and there are a total of 5925 grids, shown as Figure 12(a).Figure 12
Mesh details of: (a) Tiaoguan CR, (b) Laijiapu CR, and (c) Qigongling CR.
(a) (b) (c)(  2) Laijiapu CR.Laijiapu CR is about 8 km in length and 0.7~1.3km in width, the mesh size is 17~30 m in the horizontal direction, and along the water channel depth direction there are 10 layers grids, and there are a total of 3604 grids, shown as Figure 12(b).(  3) Chibakou CR.Chibakou CR is about 9 km in length and 0.7~1.7km in width, the mesh size is 20~50 m in the horizontal direction, and along the water channel depth direction there are 10 layers grids, and there are a total of 4905 grids, shown as Figure 12(c).


## 4.5. Verification

### 4.5.1. Verification by Flume Test

Onishi classic model test [36] was used to validate the mathematical model. The given water level at the model entrance was 0.1333 m, the average current velocity was 0.536 m/s, and the average diameter of the particles was 0.25 mm, shown as Figure 13. Verification results show that the average water lever errors between calculated value and testing value by total sediment transport model (TSTM) are less than 0.06 m in all sections, while errors are less than 0.04 m by suspended sediment transport model (SSTM). Thus, TSTM is chosen to simulate the river, shown as Figure 14. Vertical flow velocity errors are mainly produced by fluctuation of water level. The errors are less than 0.2m/s, shown as Figure 15. The predicted sediment concentration in the top layer was almost the same as the testing data, while the predicted sediment concentration of in bottom layer was vibrating within 100 g/m3 at the testing data, shown as Figure 16.Figure 13
Onishi classic model test plane layout.Figure 14
Water depth value contrast. ○: the first test set; △: the second test set; —: calculated value by TSTM. – –: calculated value by SSTM.Figure 15
Vertical flow velocity contrast. —: calculated value. ○: test value.Figure 16
Sediment concentration value contrast. —: calculated value. □: test value.Briefly, the parameters calibration results indicate that the proposed 3D model displayed a reasonable agreement with the testing data, and could be used to do comprehensive verification for fixed bed model and moving bed model.


### 4.5.2. Verification by Measured Data

The fixed bed model is vivificated by measured data of February 2010. Figures17 and 18 show the velocity comparison of cross section and longitudinal section in chosen reaches. Velocity errors of cross section are within 0.2m/s. Velocity errors of longitudinal section are within 0.4m/s.Figure 17
Velocity contrast of cross section.Figure 18
Velocity contrast of longitudinal section: (a) Tiaoguan CR, (b) Laijiapu CR, and (c) Qigongling CR.
(a) (b) (c)The moving bed model used the geomorphic data of February 2012 as the original inputs and used the geomorphic data of October 2012 as the verification materials. Figure19 shows that the predicted scouring or depositing tendency of the reaches is the same as the measured data and the predicted variation of the sediment is slightly vibrating within 20% error near the measured data, which already meet the calculation accuracy of engineering projects.Figure 19
Scouring or depositing distribution diagram from February, 2012 to October, 2012. (a) Tiaoguan CR, calculated value. (b) Tiaoguan CR, measured data. (c) Laijiapu CR, calculated value. (d) Laijiapu CR, measured data. (e) Qigongling CR, calculated value. (f) Qigongling CR, measured data.
(a) (b) (c) (d) (e) (f)According to the verification of moving bed model, the proposed 3D model behaves reliably in predict the transportation properties of water, sediment, and even the deformation of the river bed. It can be used in further study on river evolution mechanism and can be used in broader prediction of river evolution disciplines.


## 4.5.1. Verification by Flume Test

Onishi classic model test [36] was used to validate the mathematical model. The given water level at the model entrance was 0.1333 m, the average current velocity was 0.536 m/s, and the average diameter of the particles was 0.25 mm, shown as Figure 13. Verification results show that the average water lever errors between calculated value and testing value by total sediment transport model (TSTM) are less than 0.06 m in all sections, while errors are less than 0.04 m by suspended sediment transport model (SSTM). Thus, TSTM is chosen to simulate the river, shown as Figure 14. Vertical flow velocity errors are mainly produced by fluctuation of water level. The errors are less than 0.2m/s, shown as Figure 15. The predicted sediment concentration in the top layer was almost the same as the testing data, while the predicted sediment concentration of in bottom layer was vibrating within 100 g/m3 at the testing data, shown as Figure 16.Figure 13
Onishi classic model test plane layout.Figure 14
Water depth value contrast. ○: the first test set; △: the second test set; —: calculated value by TSTM. – –: calculated value by SSTM.Figure 15
Vertical flow velocity contrast. —: calculated value. ○: test value.Figure 16
Sediment concentration value contrast. —: calculated value. □: test value.Briefly, the parameters calibration results indicate that the proposed 3D model displayed a reasonable agreement with the testing data, and could be used to do comprehensive verification for fixed bed model and moving bed model.


## 4.5.2. Verification by Measured Data

The fixed bed model is vivificated by measured data of February 2010. Figures17 and 18 show the velocity comparison of cross section and longitudinal section in chosen reaches. Velocity errors of cross section are within 0.2m/s. Velocity errors of longitudinal section are within 0.4m/s.Figure 17
Velocity contrast of cross section.Figure 18
Velocity contrast of longitudinal section: (a) Tiaoguan CR, (b) Laijiapu CR, and (c) Qigongling CR.
(a) (b) (c)The moving bed model used the geomorphic data of February 2012 as the original inputs and used the geomorphic data of October 2012 as the verification materials. Figure19 shows that the predicted scouring or depositing tendency of the reaches is the same as the measured data and the predicted variation of the sediment is slightly vibrating within 20% error near the measured data, which already meet the calculation accuracy of engineering projects.Figure 19
Scouring or depositing distribution diagram from February, 2012 to October, 2012. (a) Tiaoguan CR, calculated value. (b) Tiaoguan CR, measured data. (c) Laijiapu CR, calculated value. (d) Laijiapu CR, measured data. (e) Qigongling CR, calculated value. (f) Qigongling CR, measured data.
(a) (b) (c) (d) (e) (f)According to the verification of moving bed model, the proposed 3D model behaves reliably in predict the transportation properties of water, sediment, and even the deformation of the river bed. It can be used in further study on river evolution mechanism and can be used in broader prediction of river evolution disciplines.


## 5. Conclusions

(1) Jingjiang reach is a typical meandering reach containing many local cured reaches downstream the TGP. Impounding of TGR has changed the incoming water and sediment condition in downstream river, leading to an obvious change of evolution disciplines. According to the measured data, curved reaches in Jingjiang could be divided into Type A and Type B. Type A is widely distributed. Its concave bank is scoured while convex bank is deposited before impounding of TGR, but concave bank is deposited and convex bank is scoured after impounding of TGP. Type B is mainly distributed in sharp curved reaches. Its convex bank scouring and even cutting-bend before impounding of TGR, however the evolution tendency had aggravated after impounding of TGR.(2) Evolution mechanism of curved reaches is affected by flow and sediment transport in longitudinal, transverse and vertical directions. Moreover, Jingjiang reach bed consists of fine sand. It makes transport mechanism of suspended sediment and bed load more complicated. Thus, new 2D model and 3D mathematical models which adapting to the Jingjiang curved reachs were proposed to simulate evolution disciplines of curved reaches.(3) Measured data and Onishi test results were used to verify numerical models. The 2D model verification showed that water level errors between calculated values and measure values were within 0.11 m; branch diversion ratio errors were within 2%; flow velocity errors were within 0.2 m/s.(4) The 3D model verification showed that water lever errors were within 0.04 m; vertical flow velocity errors were within 0.2 m/s; sediment concentration was within 100 g/m3; variable quantity of sediment errors was within 15.2%. Comparing calculated value with measure data, velocity errors of cross section were within 0.2 m/s; velocity errors of longitudinal section were within 0.4 m/s; sediment variation errors were within 20%. The results indicated that the proposed 3D model displayed a reasonable agreement with the testing and measure data and could be used to research the evolution mechanism of curved reaches.



---

*Source: 5980609-2018-11-12.xml*