Simulation Technology for Hydrodynamic and Water Quality in the Main Canal

Abstract

The hydrodynamic and water quality simulation technology can be used for predicting the pollutant diffusion process after a sudden water pollution accident, and for analyzing the effect of emergency operation measures. The MRP features a long route, a variety of buildings, etc.; therefore, a set of hydrodynamic and water quality models that are applicable to the main canal of the MRP was independently developed based on 1-D open canal hydrodynamic and water quality theory and with various types of buildings as inner boundaries. Through calibration and verification, these models can be applied to the simulation of hydraulic and water quality response process under any operation conditions in the main canal of the MRP.

1. Introduction

The hydraulic model is used to obtain the temporal and spatial distribution information of water level and flow by solving a numerical method, and to describe the change relations among water body quality factors with time and space under the influence of various factors [1]. After years of research and development, the hydrodynamic water quality simulation technology has become very mature and can be divided into 0-D, 1-D, 2-D, and 3-D models from the perspective of spatial dimension. The software commonly used for hydrodynamic water quality simulation are MIKE [2], EFDC [1], WASP [3], QUAL [4], DELFT3D [5], etc.

The MRP is characterized by a long route, a variety of buildings (such as inverted siphon, aqueduct, gate, etc.) and variable water distribution at dividing gates, which lead to complicated boundary conditions and the frequent operation of gate/pump groups. These characteristics result in the fact that conventional software for hydrodynamic water quality simulation cannot be well applied to the simulation of the hydraulic and water quality response process under the operation of the gate/pump groups for the main canal of the MRP [6, 7].

The water depth and water surface width of the main canal are extremely small compared with the length, and the average slope of the main canal is 1/25000 (mild slope); in case of any sudden water pollution accident, the pollutant can mix with water uniformly at a section after moving a short distance. Therefore, the pollution can be often simplified to a 1-D water quality problem for treatment, i.e., the concentration of the pollutant in the section is uniform and changes only with water flow direction. Therefore, we independently developed a set of hydrodynamic and water quality models that are applicable to the main canal of the MRP based on 1-D open canal hydrodynamic and water quality theory.

2. Hydraulic model

2.1 Control equations

When the Saint-Venant equations [8] are used to describe a 1-D unsteady canal flow model, the following assumptions should be made:

1. 1-D flow in a channel and uniform distribution of cross section of flow;

2. undulating water surface changes gradually with a very small acceleration in the vertical direction, and the pressure on cross section of flow is consistent with the distribution law of hydrostatic pressure;

3. the frictional head loss is considered only with the local head loss neglected, and the constant can be used for calculation; and

4. cosθ ≈ 1 for a small channel bottom slope.

The de Saint-Venant system of equations consists of a continuity equation and a momentum equation, and with water level Z and flow Q as variables. Its specific forms are as follows:

$$\frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=q_l$$E1

$$\frac{\partial }{\partial t}\left(\frac{Q}{A}\right)+\frac{\partial }{\partial x}\left(\beta \frac{Q^2}{2A^2}\right)+g\frac{\partial h}{\partial x}+g\left({\mathrm{S}}_f-{\mathrm{S}}_0\right)=0$$E2

where A = cross-section area; Q = cross-section flow; S₀ = channel bottom slope; S_f = hydraulic gradient; x = spatial coordinate; t = time coordinate; q_l = lateral inflow of channel of unit length; h = water depth; β = correction factor for uniform flow rate on cross section.

S_f can be determined according to a discharge modulus:

$$S_f=\frac{Q\left|Q\right|}{K^2}$$E3

where K = discharge modulus.

2.2 Equation discretization

The finite difference scheme in de Saint-Venant system of equations, the scheme applied mostly is the four-point implicit scheme raised by Preissmann, also called a space-time eccentricity scheme (Figure 1) [9]. The de Saint-Venant system of equations is dispersed by using this method in this book.

Figure 1.

Schematic diagram of Preissmann four-point implicit scheme.

See the following figure:

where $\theta$ and $\phi$ are a time weight coefficient and a space weight coefficient, respectively, with the range of values: 0 ≤ $\theta$ ≤ 1, 0 ≤ $\phi$ ≤ 1; j and j + 1 are the nodes in the space direction; i and i + 1 are the nodes in the time direction. Then, the functional values at points L, R, U, and D can be determined:

$$f_L=\theta f_j^{i+1}+\left(1-\theta \right)f_j^i$$E4

$$f_R=\theta f_{j+1}^{i+1}+\left(1-\theta \right)f_{j+1}^i$$E5

$$f_U=\phi f_{j+1}^{i+1}+\left(1-\phi \right)f_j^{i+1}$$E6

$$f_D=\phi f_{j+1}^i+\left(1-\phi \right)f_j^i$$E7

Then, the functional value for point M is:

$$f_M=\phi f_R+\left(1-\phi \right)f_L$$E8

Substituting Formulas (4) and (5) into (8), and then, we obtain:

$$f_M=\phi \left(\theta f_{j+1}^{i+1}+\left(1-\theta \right)f_{j+1}^i\right)+\left(1-\phi \right)\left(\theta f_j^{i+1}+\left(1-\theta \right)f_j^i\right)$$E9

Temporal discretization:

$${\left(\frac{\partial f}{\partial t}\right)}_M\approx \frac{f_U-f_D}{\Delta }=\frac{\left(\phi f_{j+1}^{i+1}+\left(1-\phi \right)f_j^{i+1}\right)-\left(\phi f_{j+1}^i+\left(1-\phi \right)f_j^i\right)}{\Delta }=\phi \frac{f_{j+1}^{i+1}-f_{j+1}^i}{\Delta }+\left(1-\phi \right)\frac{f_j^{i+1}-f_j^i}{\Delta }$$E10

Spatial discretization:

$${\left(\frac{\partial f}{\partial x}\right)}_M\approx \frac{f_R-f_L}{\Delta }=\frac{\left(\theta f_{j+1}^{i+1}+\left(1-\theta \right)f_{j+1}^i\right)-\left(\theta f_j^{i+1}+\left(1-\theta \right)f_j^i\right)}{\Delta }=\theta \frac{f_{j+1}^{i+1}-f_j^{i+1}}{\Delta }+\left(1-\theta \right)\frac{f_{j+1}^i-f_j^i}{\Delta }$$E11

Then, according to Formulas (9), (10) and (11), the de Saint-Venant system of equations can be dispersed into the following form:

Continuity equation:

$$\begin{array}{c}\frac{\phi }{\Delta t}\left(A_{j+1}^{i+1}-A_{j+1}^i\right)+\frac{1-\phi }{\Delta t}\left(A_j^{i+1}-A_j^i\right)+\frac{\theta }{\Delta x}\left(Q_{j+1}^{i+1}-Q_j^{i+1}\right)+\frac{1-\theta }{\Delta x}\left(\text{}{Q}_{j+1}^i-Q_j^i\right)\\ - \theta \left[\phi q_{j+1}^{i+1}+\left(1-\phi \right)q_j^{i+1}\right]-\left(1-\theta \right)\left[\phi q_{j+1}^i+\left(1-\phi \right)q_j^i\right]=0\end{array}$$E12

Momentum equation:

$$\begin{array}{l}\frac{\phi }{\Delta t}\left(\frac{Q_{j+1}^{i+1}}{A_{j+1}^{i+1}}-\frac{Q_{j+1}^i}{A_{j+1}^i}\right)+\frac{1-\phi }{\Delta t}\left(\frac{Q_j^{i+1}}{A_j^{i+1}}-\frac{Q_j^i}{A_j^i}\right)+\frac{\theta }{\Delta x}\left[\frac{1}{2}{\left({\beta }_{j+1}^{i+1}\frac{Q_{j+1}^{i+1}}{A_{j+1}^{i+1}}\right)}^2-\frac{1}{2}{\left({\beta }_j^{i+1}\frac{Q_j^{i+1}}{A_j^{i+1}}\right)}^2\right]\hfill \\ +\frac{1-\theta }{\Delta x}\left[\frac{{\beta }_{j+1}^i}{2}{\left(\frac{Q_{j+1}^i}{A_{j+1}^i}\right)}^2-\frac{{\beta }_j^i}{2}{\left(\frac{Q_j^i}{A_j^i}\right)}^2\right]+\frac{\theta g}{\Delta x}\left(h_{j+1}^{i+1}-h_j^{i+1}\right)+\frac{\left(1-\theta \right)g}{\Delta x}\left(h_{j+1}^i-h_j^i\right)\hfill \\ \begin{array}{l}+ \theta g\left[S_{f,j+1}^{i+1}+\left(1-{\phi }_R\right)S_{f,j}^{i+1}\right]+\left(1-\theta \right)g\left[S_{f,j+1}^i+\left(1-{\phi }_R\right)S_{f,j}^i\right]=0\hfill \end{array}\hfill \end{array}$$E13

where ${\phi }_R$ is the space weight coefficient of hydraulic gradient, with a value different from that of $\phi$.

When the water in a channel flows slowly, the characteristic root symbols of the control equation are opposite, i.e., ${\lambda }_1>$ 0 or ${\lambda }_2$ < 0; therefore, $C_r={\lambda }_1\frac{\Delta \mathit{t}}{\Delta }>$ 0 or $C_r={\lambda }_2\frac{\Delta \mathit{t}}{\Delta }$ < 0. Because $C_r>$ 0 and $C_r$ < 0 exist at the same time, two discrete equations are not likely to go into an unconditional stability state at the same time. The stability condition for discrete equations is as follows:

$$C_r=\left(u-\sqrt{\mathit{gh}}\right)\frac{\Delta \mathit{t}}{\Delta }>\frac{\phi -0.5}{0.5-\theta }$$E14

2.3 Discrete equation linearization

The aforesaid (12) and (13) are the continuity equation and momentum equation that have been dispersed. The discrete equations after discretization are nonlinear, so they still need linearization and are determined. In this book, discrete equations are determined by solving water level and flow increment, i.e., $\Delta \mathit{h}$ and $\Delta \mathit{Q}$. Firstly, area and flow are transformed into an increment form:

$$A_j^{i+1}=A_j^{\ast }+\Delta A_j=A_j^{\ast }+B_j^{\ast }\Delta h_j$$E15

$$Q_j^{i+1}=Q_j^{\ast }+\Delta Q_j$$E16

where * stands for the value of the variable at the previous cycle step; $\Delta \mathit{A}$, $\Delta \mathit{h}$, and $\Delta \mathit{Q}$ are the flow area, channel water depth, and flow increment, respectively; and B is the water surface width.

After substituting the aforesaid Formulas (15) and (16) into the dispersed continuity equation (12), we obtain:

$$a_j\Delta h_j+b_j\Delta Q_j+c_j\Delta h_{j+1}+d_j\Delta Q_{j+1}=p_j$$E17

where

$$a_j=\left(1-\phi \right)B_j^{\ast }/\Delta \mathit{t}$$E18

$$b_j=-\theta /\Delta \mathit{x}$$E19

$$c_j=\phi B_{j+1}^{\ast }/\Delta \mathit{t}$$E20

$$d_j=\theta /\Delta \mathit{t}$$E21

$$\begin{array}{ll}{p}_j=& -\frac{\phi }{\Delta }\left(A_{j+1}^{\ast }-A_{j+1}^i\right)-\frac{1-\phi }{\Delta }\left(A_j^{\ast }-A_j^i\right)-\frac{\theta }{\Delta }\left(Q_{j+1}^{\ast }-Q_j^{\ast }\right)-\frac{1-\theta }{\Delta }\left(Q_{j+1}^i-Q_j^i\right)\\ & + \theta \left[\phi q_{j+1}^{i+1}+\left(1-\phi \right)q_j^{i+1}\right]+\left(1-\theta \right)\left[\phi q_{j+1}^i+\left(1-\phi \right)q_j^i\right]\end{array}$$E22

During linearization of the momentum equation, the following formulas are needed:

$${\left(Q_j^{i+1}\right)}^2={\left(Q_j^{\ast }\right)}^2+2Q_j^{\ast }\Delta Q_j$$E23

$$\frac{1}{{\left(K_j^{i+1}\right)}^2}=\frac{1}{{\left(K_j^{\ast }\right)}^2}-\frac{2}{{\left(K_j^{\ast }\right)}^3}{\left(\frac{\partial K}{\partial h}\right)}_j^{\ast }\Delta h_j$$E24

$$S_{f,j}^{i+1}=S_{f,j}^{\ast }+\frac{2\left|Q_j^{\ast }\right|}{{\left(K_j^{\ast }\right)}^2}\Delta Q_j-\frac{2S_{f,j}^{\ast }}{K_j^{\ast }}{\left(\frac{\partial K}{\partial h}\right)}_j^{\ast }\Delta h_j$$E25

$$\frac{Q_j^{i+1}}{A_j^{i+1}}=\frac{Q_j^{\ast }}{A_j^{\ast }}+\frac{1}{A_j^{\ast }}\Delta Q_j-\frac{Q_j^{\ast }{B}_j^{\ast }}{{\left(A_j^{\ast }\right)}^2}\Delta h_j$$E26

$${\left(\frac{Q_j^{i+1}}{A_j^{i+1}}\right)}^2={\left(\frac{Q_j^{\ast }}{A_j^{\ast }}\right)}^2+\frac{2Q_j^{\ast }}{{\left(A_j^{\ast }\right)}^2}\Delta Q_j-\frac{2{\left(Q_j^{\ast }\right)}^2{B}_j^{\ast }}{{\left(A_j^{\ast }\right)}^3}\Delta h_j$$E27

After substituting Formulas (23)–(27) into the dispersed momentum equation (13), we obtain a linearized momentum equation:

$$a_{j+1}\Delta h_j+b_{j+1}\Delta Q_j+c_{j+1}\Delta h_{j+1}+d_{j+1}\Delta Q_{j+1}=p_{j+1}$$E28

where

$$a_{j+1}=-\frac{1-\phi }{\Delta }\frac{Q_j^{\ast }{B}_j^{\ast }}{{\left(A_j^{\ast }\right)}^2}+\frac{\theta }{\Delta }\frac{{\left(Q_j^{\ast }\right)}^2{B}_j^{\ast }}{{\left(A_j^{\ast }\right)}^3}-\frac{\mathit{\theta g}}{\Delta }-2\theta \left(1-{\phi }_R\right)g\frac{S_{f,j}^{\ast }}{K_j^{\ast }}{\left(\frac{\partial K}{\partial h}\right)}_j^{\ast }$$E29

$$b_{j+1}=\frac{1-\phi }{\Delta }\frac{1}{A_j^{\ast }}-\frac{\theta }{\Delta }\frac{Q_j^{\ast }}{{\left(A_j^{\ast }\right)}^2}+2\theta \left(1-{\phi }_R\right)g\frac{\left|Q_j^{\ast }\right|}{{\left(K_j^{\ast }\right)}^2}$$E30

$$c_{j+1}=-\frac{\phi }{\Delta }\frac{Q_{j+1}^{\ast }{B}_{j+1}^{\ast }}{{\left(A_{j+1}^{\ast }\right)}^2}-\frac{\theta }{\Delta }\frac{{\left(Q_{j+1}^{\ast }\right)}^2{B}_{j+1}^{\ast }}{{\left(A_{j+1}^{\ast }\right)}^3}+\frac{\mathit{\theta g}}{\Delta }-2{\mathit{\theta \phi }}_{R}g\frac{S_{f,j+1}^{\ast }}{K_{j+1}^{\ast }}{\left(\frac{\partial K}{\partial h}\right)}_{j+1}^{\ast }$$E31

$$d_{j+1}=\frac{\phi }{\Delta }\frac{1}{A_{j+1}^{\ast }}+\frac{\theta }{\Delta }\frac{Q_{j+1}^{\ast }}{{\left(A_{j+1}^{\ast }\right)}^2}+2{\mathit{\theta \phi }}_{R}g\frac{\left|Q_{j+1}^{\ast }\right|}{{\left(K_{j+1}^{\ast }\right)}^2}$$E32

$$\begin{array}{ll}{p}_{j+1}& =-\frac{\phi }{\Delta }\left(\frac{Q_{j+1}^{\ast }}{A_{j+1}^{\ast }}-\frac{Q_{j+1}^i}{A_{j+1}^i}\right)-\frac{1-\phi }{\Delta }\left(\frac{Q_j^{\ast }}{A_j^{\ast }}-\frac{Q_j^i}{A_j^i}\right)-\frac{\theta }{\Delta }\left[\frac{1}{2}{\left(\frac{Q_{j+1}^{\ast }}{A_{j+1}^{\ast }}\right)}^2-\frac{1}{2}{\left(\frac{Q_j^{\ast }}{A_j^{\ast }}\right)}^2\right]\\ & -\frac{1-\theta }{\Delta }\left[\frac{1}{2}{\left(\frac{Q_{j+1}^i}{A_{j+1}^i}\right)}^2-\frac{1}{2}{\left(\frac{Q_j^i}{A_j^i}\right)}^2\right]-\frac{\mathit{\theta g}}{\Delta }\left(h_{j+1}^{\ast }-h_j^{\ast }\right)-\frac{\left(1-\theta \right)g}{\Delta }\left(h_{j+1}^i-h_j^i\right)\\ & -\theta {\mathit{gS}}_{f,j+1}^{\ast }-\left(1-\theta \right){\mathit{gS}}_{f,j}^{\ast }\end{array}$$E33

The aforesaid Formulas (17) and (28) are the continuity equation and the momentum equation in de Saint-Venant system of equations, respectively, after linearization via an increment method.

2.4 Boundary conditions

With the main canal of the MRP as a whole system, if all the hydraulic elements are to be determined via solving the de Saint-Venant system of equations, the various types of boundary conditions need to be defined. Normally, boundary conditions are divided into outer and inner ones. With regard to the actual situation of the general main canal of the MRP, the upstream water depth boundary and downstream flow boundary are the outer boundary conditions for value simulation, while the water level discharge process of structures such as the dividing gates, transition sections, inverted siphons, check gates, and so on in the channel system are the inner boundary conditions for value simulation.

2.4.1 Outer boundary conditions

The characteristics of upstream and downstream water connected with the general main canal are considered as the outer boundary conditions of the system. They can be roughly divided into three types: water depth boundary, flow boundary, and water depth-flow boundary. Because the system has a water connection with the upstream and downstream, the outer boundary conditions of the system should be considered, respectively.

2.4.1.1 Upstream water depth boundary

The upstream water depth generally changes with time, so it can be considered as the function of time t, namely:

$$h_1=h\left(t\right)$$E34

where $h_1$ is the upstream water depth.

Transform the aforesaid Formula (34) to an increment form, and then, we obtain:

$$\Delta h_1=h_1^{i+1}-h_1^{\ast }$$E35

where $h_1^{i+1}$ is the upstream water depth at time (i + 1); $h_1^{\ast }$ is the value of the variable at the previous cycle step.

Transform the aforesaid Formula (35) into a form identical to linearized de Saint-Venant system of equations, and then, we obtain:

$$c_1\Delta h_1+d_1\Delta Q_1=p_1$$E36

where $c_1$ = 1.0, $d_1$ = 0, and $p_1$ = $h_1^{i+1}-h_1^{\ast }$.

2.4.1.2 Upstream flow boundary

The upstream flow generally changes with time, so it can be considered as the function of time t, namely:

$$Q_1=Q\left(t\right)$$E37

where $Q_1$ is the upstream flow.

Transform the aforesaid Formula (34) to an increment form, and then, we obtain:

$$\Delta Q_1=Q_1^{i+1}-Q_1^{\ast }$$E38

where $Q_1^{i+1}$ is the upstream flow at time (i + 1) and $Q_1^{\ast }$ is the value of the variable at the previous cycle step.

Transform the aforesaid Formula (38) into a form identical to linearized de Saint-Venant system of equations, and then, we obtain:

$$c_1\Delta h_1+d_1\Delta Q_1=p_1$$E39

where $c_1$ = 0, $d_1$ = 1.0, and $p_1$ =$Q_1^{i+1}-Q_1^{\ast }$.

2.4.1.3 Downstream water depth boundary

Like the upstream water depth boundary, the downstream water depth generally changes with time too, so it can be considered as the function of time t, namely:

$$h_n=h\left(t\right)$$E40

where $h_n$ is the downstream water depth.

Transform the aforesaid Formula (40) to an increment form, and then, we obtain:

$$\Delta h_n=h_n^{i+1}-h_n^{\ast }$$E41

where $h_n^{i+1}$ is the downstream water depth at time (i + 1); $h_n^{\ast }$ is the value of the variable at the previous cycle step.

Transform the aforesaid Formula (41) into a form identical to linearized de Saint-Venant system of equations, and then, we obtain:

$$a_n\Delta h_n+b_n\Delta Q_n=p_n$$E42

where $a_n$ =1.0, $b_n$ = 0, and $p_n$ = $h_n^{i+1}-h_n^{\ast }$.

2.4.1.4 Downstream flow boundary

Like the upstream flow boundary, the downstream flow generally changes with time too, so it can be considered as the function of time t, namely:

$$Q_n=Q\left(t\right)$$E43

where $Q_n$ is the upstream flow.

Transform the aforesaid Formula (43) to an increment form, and then, we obtain:

$$\Delta Q_n=Q_n^{i+1}-Q_n^{\ast }$$E44

where $Q_n^{i+1}$ is the upstream flow at time (i + 1); $Q_n^{\ast }$ is the value of the variable at the previous cycle step.

Transform the aforesaid Formula (44) into a form identical to linearized de Saint-Venant system of equations, and then, we obtain:

$$a_n\Delta h_n+b_n\Delta Q_n=p_n$$E45

where $a_n$ = 0, $b_n$ = 1.0, and $p_n$ = $Q_n^{i+1}-Q_n^{\ast }$.

2.4.1.5 Upstream water depth-flow boundary

The upstream water depth and flow are considered conforming to the following functional relation:

$$Q_1=f\left(h_1\right)$$E46

where $h_1$ is the upstream water depth and $Q_1$ is the upstream flow.

Transform the left of the aforesaid formula into an increment form, and make Taylor’s series expansion of the right, and then, we obtain:

$$\Delta Q_1+Q_1^{\ast }=f\left(h_1^{\ast }\right)+\frac{\mathit{df}}{\mathit{dh}}\Delta h_1$$E47

Transform the aforesaid Formula (47) into a form identical to linearized de Saint-Venant system of equations, and then, we obtain:

$$c_1\Delta h_1+d_1\Delta Q_1=p_1$$E48

where $c_1=-\frac{\mathit{df}}{\mathit{dh}}$ and $d_1$ = 1.0; and $p_1$ =$f\left(h_1^{\ast }\right)-Q_1^{\ast }$.

2.4.1.6 Downstream water depth-flow boundary

Like the upstream treatment, the downstream water depth and flow are also considered conforming to the following functional relation:

$$Q_n=f\left(h_n\right)$$E49

where $h_n$ is the downstream water depth and $Q_n$ is the downstream flow.

Transform the left of the aforesaid formula into an increment form and make Taylor’s series expansion of the right, and then, we obtain:

$$\Delta Q_n+Q_n^{\ast }=f\left(h_n^{\ast }\right)+\frac{\mathit{df}}{\mathit{dh}}\Delta h_n$$E50

Transform the aforesaid Formula (50) into a form identical to linearized de Saint-Venant system of equations, and then, we obtain:

$$a_n\Delta h_n+b_n\Delta Q_n=p_n$$E51

where $a_n=-\frac{\mathit{df}}{\mathit{dh}}$ and $b_n$ = 1.0; and $p_n$ = $f\left(h_n^{\ast }\right)-Q_n^{\ast }$.

2.4.2 Inner boundary conditions

Because the MRP involves a variety of buildings along the main route, the whole general main canal can be considered to be a channel pool system which connects different hydraulic structures in series, with all check gates used to divide different channel pools. Therefore, when considering inner boundary conditions, we mainly focus on four structures inside the channel system: dividing gate (release sluice), transition section, inverted siphon, and check gate.

2.4.2.1 Boundary conditions of dividing gate

The situation of water flow at the dividing gate (Figure 2) is described by using the water balance equation and the energy conservation equation:

$$Q_j=Q_{j+1}+Q_f$$E52

$$h_j+Z_j+\frac{1}{2g}{\left(\frac{Q_j}{A_j}\right)}^2=h_{j+1}+Z_{j+1}+\frac{1}{2g}{\left(\frac{Q_{j+1}}{A_{j+1}}\right)}^2$$E53

Figure 2.

Schematic diagram of dividing gate.

where j and j + 1 are the upstream and downstream cross section numbers of the dividing gate, respectively; $Q_j$ and $Q_{j+1}$ are the upstream and downstream cross section flows, respectively; $Q_f$ is the divided water flow at the dividing gate; $Z_j$ and $Z_{j+1}$ are the elevation of the channel bottom of upstream and downstream cross sections at the dividing gate.

After linearization of the aforesaid formula, we obtain:

$$\mathrm{\Delta }{Q}_j-\mathrm{\Delta }{Q}_{j+1}=Q_{j+1}^{\ast }-Q_j^{\ast }+Q_f$$E54

$$\begin{array}{l}\left(1-\frac{{\left(Q_j^{\ast }\right)}^2{B}_j^{\ast }}{g{\left(A_j^{\ast }\right)}^3}\right)\mathrm{\Delta }{h}_j+\frac{Q_j^{\ast }}{g{\left(A_j^{\ast }\right)}^2}\mathrm{\Delta }{Q}_j-\left(1-\frac{{\left(Q_{j+1}^{\ast }\right)}^2{B}_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}\right)\mathrm{\Delta }{h}_{j+1}-\frac{Q_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^2}\mathrm{\Delta }{Q}_{j+1}\\ =-h_j^{\ast }-Z_j-\frac{1}{2g}{\left(\frac{Q_j^{\ast }}{A_j^{\ast }}\right)}^2+h_{j+1}^{\ast }+Z_{j+1}+\frac{1}{2g}{\left(\frac{Q_{j+1}^{\ast }}{A_{j+1}^{\ast }}\right)}^2\end{array}$$E55

Transform them into the form of Formulas (17) and (28), and all respective factors are as follows:

$$a_j=0,b_j=1.0,c_j=0,d_j=-1.0,p_j=Q_{j+1}^{\ast }-Q_j^{\ast }+Q_f.$$E56

$$a_{j+1}=\left(1-\frac{{\left(Q_j^{\ast }\right)}^2{B}_j^{\ast }}{g{\left(A_j^{\ast }\right)}^3}\right),b_{j+1}=\frac{Q_j^{\ast }}{g{\left(A_j^{\ast }\right)}^2},c_{j+1}=-\left(1-\frac{{\left(Q_{j+1}^{\ast }\right)}^2{B}_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}\right),d_{j+1}=-\frac{Q_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^2}$$E57

$$p_{j+1}=-h_j^{\ast }-Z_j-\frac{1}{2g}{\left(\frac{Q_j^{\ast }}{A_j^{\ast }}\right)}^2+h_{j+1}^{\ast }+Z_{j+1}+\frac{1}{2g}{\left(\frac{Q_{j+1}^{\ast }}{A_{j+1}^{\ast }}\right)}^2.$$E58

2.4.2.2 Boundary conditions of transition section

Like the treatment of the dividing gate, the water balance equation and energy conservation equation used to describe hydraulic characteristics at the transition section (Figure 3) are shown below:

$$Q_j=Q_{j+1}$$E59

$$\begin{array}{l}{h}_j+Z_j+\frac{1}{2g}{\left(\frac{Q_j}{A_j}\right)}^2=h_{j+1}+Z_{j+1}+\frac{1}{2g}{\left(\frac{Q_{j+1}}{A_{j+1}}\right)}^2\\ +\zeta \frac{1}{2g}\mid \left({\left(\frac{Q_j}{A_j}\right)}^2-{\left(\frac{Q_{j+1}}{A_{j+1}}\right)}^2\right)\mid +\xi \frac{1}{2g}{\left(\frac{Q_j}{A_j}\right)}^2\end{array}$$E60

Figure 3.

Schematic diagram of transition section.

where $\zeta$ is a local loss coefficient and $\xi$ is an other loss coefficient. However, when the hydraulic characteristics of the transition section are considered, the transition section should be divided into two cases: converge section and divergence section, which should be treated differently.

After linearization of the aforesaid formula, we obtain:

$$\mathrm{\Delta }{Q}_j-\mathrm{\Delta }{Q}_{j+1}=Q_{j+1}^{\ast }-Q_j^{\ast }$$E61

$$\begin{array}{l}\left(1-\left(1-\zeta \right)\frac{{\left(Q_j^{\ast }\right)}^2{B}_j^{\ast }}{g{\left(A_j^{\ast }\right)}^3}\right)\mathrm{\Delta }{h}_j+\frac{\left(1-\zeta \right)Q_j^{\ast }}{g{\left(A_j^{\ast }\right)}^2}\mathrm{\Delta }{Q}_j-\left(\frac{\left(1-\zeta \right)Q_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^2}+\frac{\xi \mid Q_{j+1}^{\ast }\mid }{g{\left(A_{j+1}^{\ast }\right)}^2}\right)\mathrm{\Delta }{Q}_{j+1}\\ -\left(1-\left(1-\zeta \right)\frac{{\left(Q_{j+1}^{\ast }\right)}^2{B}_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}-\frac{\xi Q_{j+1}^{\ast }\mid Q_{j+1}^{\ast }\mid B_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}\right)\mathrm{\Delta }{h}_{j+1}\\ =-h_j^{\ast }-Z_j-\frac{1-\zeta }{2g}{\left(\frac{Q_j^{\ast }}{A_j^{\ast }}\right)}^2+h_{j+1}^{\ast }+Z_{j+1}+\frac{1-\zeta }{2g}{\left(\frac{Q_{j+1}^{\ast }}{A_{j+1}^{\ast }}\right)}^2+\frac{\xi Q_{j+1}^{\ast }\mid Q_{j+1}^{\ast }\mid }{2g{\left(A_{j+1}^{\ast }\right)}^2}\end{array}$$E62

$$\begin{array}{l}\left(1-\left(1+\zeta \right)\frac{{\left(Q_j^{\ast }\right)}^2{B}_j^{\ast }}{g{\left(A_j^{\ast }\right)}^3}\right)\mathrm{\Delta }{h}_j+\frac{\left(1+\zeta \right)Q_j^{\ast }}{g{\left(A_j^{\ast }\right)}^2}\mathrm{\Delta }{Q}_j-\left(\frac{\left(1+\zeta \right)Q_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^2}+\frac{\xi \mid Q_{j+1}^{\ast }\mid }{g{\left(A_{j+1}^{\ast }\right)}^2}\right)\mathrm{\Delta }{Q}_{j+1}\\ -\left(1-\left(1+\zeta \right)\frac{{\left(Q_{j+1}^{\ast }\right)}^2{B}_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}-\frac{\xi Q_{j+1}^{\ast }\mid Q_{j+1}^{\ast }\mid B_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}\right)\mathrm{\Delta }{h}_{j+1}\\ =-h_j^{\ast }-Z_j-\frac{1+\zeta }{2g}{\left(\frac{Q_j^{\ast }}{A_j^{\ast }}\right)}^2+h_{j+1}^{\ast }+Z_{j+1}+\frac{1+\zeta }{2g}{\left(\frac{Q_{j+1}^{\ast }}{A_{j+1}^{\ast }}\right)}^2+\frac{\xi Q_{j+1}^{\ast }\mid Q_{j+1}^{\ast }\mid }{2g{\left(A_{j+1}^{\ast }\right)}^2}\end{array}$$E63

Transform them into the form of Formulas (17) and (28), and all respective factors are as follows:

$$a_j=0,b_j=1.0,c_j=0,d_j=-1.0,p_j=Q_{j+1}^{\ast }-Q_j^{\ast }.$$E64

The coefficient for the converge section:

$$\begin{array}{l}{a}_{j+1}=\left(1-\left(1-\zeta \right)\frac{{\left(Q_j^{\ast }\right)}^2{B}_j^{\ast }}{g{\left(A_j^{\ast }\right)}^3}\right),b_{j+1}=\frac{\left(1-\zeta \right)Q_j^{\ast }}{g{\left(A_j^{\ast }\right)}^2},c_{j+1}=-\left(\frac{\left(1-\zeta \right)Q_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^2}+\frac{\xi \mid Q_{j+1}^{\ast }\mid }{g{\left(A_{j+1}^{\ast }\right)}^2}\right),\hfill \\ d_{j+1}=-\left(1-\left(1-\zeta \right)\frac{{\left(Q_{j+1}^{\ast }\right)}^2{B}_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}-\frac{\xi Q_{j+1}^{\ast }\mid Q_{j+1}^{\ast }\mid B_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}\right)\hfill \end{array}$$E65

$$p_{j+1}=-h_j^{\ast }-Z_j-\frac{1-\zeta }{2g}{\left(\frac{Q_j^{\ast }}{A_j^{\ast }}\right)}^2+h_{j+1}^{\ast }+Z_{j+1}+\frac{1-\zeta }{2g}{\left(\frac{Q_{j+1}^{\ast }}{A_{j+1}^{\ast }}\right)}^2+\frac{\xi Q_{j+1}^{\ast }\mid Q_{j+1}^{\ast }\mid }{2g{\left(A_{j+1}^{\ast }\right)}^2}.$$E66

The coefficient for the divergence section:

$$\begin{array}{l}{a}_{j+1}=\left(1-\left(1+\zeta \right)\frac{{\left(Q_j^{\ast }\right)}^2{B}_j^{\ast }}{g{\left(A_j^{\ast }\right)}^3}\right),b_{j+1}=\frac{\left(1+\zeta \right)Q_j^{\ast }}{g{\left(A_j^{\ast }\right)}^2},c_{j+1}=-\left(\frac{\left(1+\zeta \right)Q_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^2}+\frac{\xi \mid Q_{j+1}^{\ast }\mid }{g{\left(A_{j+1}^{\ast }\right)}^2}\right),\hfill \\ d_{j+1}=-\left(1-\left(1+\zeta \right)\frac{{\left(Q_{j+1}^{\ast }\right)}^2{B}_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}-\frac{\xi Q_{j+1}^{\ast }\mid Q_{j+1}^{\ast }\mid B_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}\right)\hfill \end{array}$$E67

$$p_{j+1}=-h_j^{\ast }-Z_j-\frac{1+\zeta }{2g}{\left(\frac{Q_j^{\ast }}{A_j^{\ast }}\right)}^2+h_{j+1}^{\ast }+Z_{j+1}+\frac{1+\zeta }{2g}{\left(\frac{Q_{j+1}^{\ast }}{A_{j+1}^{\ast }}\right)}^2+\frac{\xi Q_{j+1}^{\ast }\mid Q_{j+1}^{\ast }\mid }{2g{\left(A_{j+1}^{\ast }\right)}^2}.$$E68

2.4.2.3 Boundary conditions of inverted siphon

When describing the hydraulic characteristics at the inverted siphon (Figure 4), we consider that the water head difference at the inlet and outlet consists of the local head loss and middle frictional head loss; then, the water balance equation and energy conservation equation used to describe the hydraulic characteristics are indicated below:

$$Q_j=Q_{j+1}$$E69

$$\begin{array}{l}{h}_j+Z_j+\frac{1}{2g}{\left(\frac{Q_j}{A_j}\right)}^2=h_{j+1}+Z_{j+1}+\frac{1}{2g}{\left(\frac{Q_{j+1}}{A_{j+1}}\right)}^2\\ +{\zeta }_1\frac{1}{2g}{\left(\frac{Q_j}{A_j}\right)}^2+{\zeta }_2\frac{1}{2g}{\left(\frac{Q_{j+1}}{A_{j+1}}\right)}^2+\frac{Q_{j+1}\mid Q_{j+1}\mid }{{\left(K_{j+1}\right)}^2}L\end{array}$$E70

Figure 4.

Schematic diagram of inverted siphon.

where ${\zeta }_1$ and ${\zeta }_2$ are the loss coefficient at the inlet and outlet; and L is the middle length of the inverted siphon.

After linearization of the aforesaid formula, we obtain:

$$\mathrm{\Delta }{Q}_j-\mathrm{\Delta }{Q}_{j+1}=Q_{j+1}^{\ast }-Q_j^{\ast }$$E71

$$\begin{array}{l}\left(1-\left(1-{\zeta }_1\right)\frac{{\left(Q_j^{\ast }\right)}^2{B}_j^{\ast }}{g{\left(A_j^{\ast }\right)}^3}\right)\mathrm{\Delta }{h}_j-\left(\frac{\left(1+{\zeta }_2\right)Q_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^2}+\frac{2\mid Q_{j+1}^{\ast }\mid }{{\left(K_{j+1}^{\ast }\right)}^2}L\right)\mathrm{\Delta }{Q}_{j+1}\\ -\left(1-\left(1+{\zeta }_2\right)\frac{{\left(Q_{j+1}^{\ast }\right)}^2{B}_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}\right)\mathrm{\Delta }{h}_{j+1}+\frac{\left(1-{\zeta }_1\right)Q_j^{\ast }}{g{\left(A_j^{\ast }\right)}^2}\mathrm{\Delta }{Q}_j\\ =-h_j^{\ast }-Z_j-\frac{1-{\zeta }_1}{2g}{\left(\frac{Q_j^{\ast }}{A_j^{\ast }}\right)}^2+h_{j+1}^{\ast }+Z_{j+1}+\frac{1+{\zeta }_2}{2g}{\left(\frac{Q_{j+1}^{\ast }}{A_{j+1}^{\ast }}\right)}^2+\frac{Q_{j+1}^{\ast }\mid Q_{j+1}^{\ast }\mid }{{\left(K_{j+1}^{\ast }\right)}^2}L\end{array}$$E72

where $K=\frac{A}{n}{R}^{\frac{2}{3}}$.

Transform them into the form of Formulas (17) and (28), and all respective factors are as follows:

$$a_j=0,b_j=1.0,c_j=0,d_j=-1.0,p_j=Q_{j+1}^{\ast }-Q_j^{\ast }.$$E73

$$\begin{array}{l}{a}_{j+1}=\left(1-\left(1-{\zeta }_1\right)\frac{{\left(Q_j^{\ast }\right)}^2{B}_j^{\ast }}{g{\left(A_j^{\ast }\right)}^3}\right),b_{j+1}=\frac{\left(1-{\zeta }_1\right)Q_j^{\ast }}{g{\left(A_j^{\ast }\right)}^2},c_{j+1}=-\left(1-\left(1+{\zeta }_2\right)\frac{{\left(Q_{j+1}^{\ast }\right)}^2{B}_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}\right),\hfill \\ d_{j+1}=-\left(1-\left(1+{\zeta }_2\right)\frac{{\left(Q_{j+1}^{\ast }\right)}^2{B}_{j+1}^{\ast }}{g{\left(A_{j+1}^{\ast }\right)}^3}\right)\hfill \end{array}$$E74

$$p_{j+1}=-h_j^{\ast }-Z_j-\frac{1-{\zeta }_1}{2g}{\left(\frac{Q_j^{\ast }}{A_j^{\ast }}\right)}^2+h_{j+1}^{\ast }+Z_{j+1}+\frac{1+{\zeta }_2}{2g}{\left(\frac{Q_{j+1}^{\ast }}{A_{j+1}^{\ast }}\right)}^2+\frac{Q_{j+1}^{\ast }\mid Q_{j+1}^{\ast }\mid }{{\left(K_{j+1}^{\ast }\right)}^2}L$$E75

2.4.2.4 Boundary conditions of check gate

The equations used to describe the hydraulic characteristics of water flowing through the check gate (Figure 5) are shown below:

$$Q_j=Q_{j+1}$$E76

$$Q_j={\sigma }_c{\sigma }_s\mathit{mbf}\left(H_j,H_{j+1},e\right)$$E77

Figure 5.

Schematic diagram of arc check gate.

where ${\sigma }_c$ is the side-contract coefficient; ${\sigma }_s$ is the submergence coefficient; $m$ is the flow coefficient; $b$ is the width of sluice; $H_j$ and $H_{j+1}$ are the downstream water levels; $e$ is the opening of sluice.

Transform Formulas (76) and (77) to an increment form, and then, we obtain:

$$\mathrm{\Delta }{Q}_j-\mathrm{\Delta }{Q}_{j+1}=Q_{j+1}^{\ast }-Q_j^{\ast }$$E78

$${\sigma }_c{\sigma }_s\mathit{mb}\frac{\partial f^{\ast }}{\partial H_j^{\ast }}\mathrm{\Delta }{H}_j-\mathrm{\Delta }{Q}_j+{\sigma }_c{\sigma }_s\mathit{mb}\frac{\partial f^{\ast }}{\partial H_{j+1}^{\ast }}\mathrm{\Delta }{H}_{j+1}=Q_j^{\ast }-{\sigma }_c{\sigma }_s{\mathit{mbf}}^{\ast }$$E79

where $f^{\ast }={\sigma }_c{\sigma }_s\mathit{mbf}\left(H_j^{\ast },H_{j+1}^{\ast },e\right)$.

Transform them into the form of Formulas (17) and (28), and all respective factors are as follows:

$$a_j=0,b_j=1.0,c_j=0,d_j=-1.0,p_j=Q_{j+1}^{\ast }-Q_j^{\ast }.$$E80

$$a_{j+1}={\sigma }_c{\sigma }_s\mathit{mb}\frac{\partial f^{\ast }}{\partial H_j^{\ast }},b_{j+1}=-1,c_{j+1}={\sigma }_c{\sigma }_s\mathit{mb}\frac{\partial f^{\ast }}{\partial H_{j+1}^{\ast }},d_{j+1}=0$$E81

$$p_{j+1}=Q_j^{\ast }-{\sigma }_c{\sigma }_s{\mathit{mbf}}^{\ast }$$E82

2.5 Equation solution

Assuming that a channel has n cross sections; each of which has two variables—flow and water level. That is to say, it has totally 2n variables. The channel is divided by n cross sections into n−1 channel section. According to the de Saint-Venant system of equations, two equations can be set up for each channel section, and totally 2(n−1) equations can be obtained; and after combining them with upstream and downstream boundary conditions, totally 2n equations can be obtained, forming a closed system of equations. The coefficient of each equation in the system of equations can be calculated by an initial value and a previous iterative value; therefore, the system of equations is a constant coefficient linear system of equations which can be transformed into a matrix form.

Formulas (17) and (28) are transformed into the following forms:

$$a_{2i}\Delta h_i+b_{2i}\Delta Q_i+c_{2i}\Delta h_{i+1}+d_{2i}\Delta Q_{i+1}=D_{2i}$$E83

$$a_{i+1}\Delta h_i+b_{2i+1}\Delta Q_i+c_{2i+1}\Delta h_{i+1}+d_{2i}\Delta Q_{i+1}=D_{2i+1}$$E84

where all coefficients correspond to Formulas (17) and (28); i = 0, 1, 2, 3…n−1, n, the number of each channel section; i = 0 and i = n stand for the upstream and downstream boundary conditions, respectively. Then, the amended system of equations can be transformed into the following matrix form:

$$\mathit{AX}=D$$E85

where

$$A=\left\{\begin{array}{cccccccc}{c}_1& d_1& & & & & & \\ a_2& b_2& c_2& d_2& & & & \\ a_3& b_3& c_3& d_3& & & & \\ & & a_4& b_4& c_4& d_4& & \\ & & a_5& b_5& c_5& d_5& & \\ & & ⋮& ⋮& ⋮& ⋮& & \\ & & & & a_{2n-2}& b_{2n-2}& c_{2n-2}& d_{2n-2}\\ & & & & a_{2n-1}& b_{2n-1}& c_{2n-1}& d_{2n-1}\\ & & & & & & a_{2n}& b_{2n}\end{array}\right\},X=\left[\begin{array}{c}\mathrm{\Delta }{h}_1\\ \mathrm{\Delta }{Q}_1\\ \mathrm{\Delta }{h}_2\\ \mathrm{\Delta }{Q}_2\\ \mathrm{\Delta }{h}_3\\ ⋮\\ \mathrm{\Delta }{Q}_{n-1}\\ \mathrm{\Delta }{h}_n\\ \mathrm{\Delta }{Q}_n\end{array}\right],D=\left[\begin{array}{c}{D}_0\\ D_1\\ D_2\\ D_3\\ D_4\\ ⋮\\ D_{2n-2}\\ D_{2n-1}\\ D_{2n}\end{array}\right].$$E86

After such transformation, the system of equations can be solved.

When calculating the water depth and flow of each cross section of a channel with a known time $t_0$, the following process can be used to determine the water depth and flow at time $t_0+\Delta \mathit{t}$:

1. assume the initial value of the water depth and flow of each cross section of the channel at time $t_0+\Delta \mathit{t}$;

2. calculate the value of all elements in matrices $A$ and $D$;

3. solve equation $\mathit{AX}=D$ and obtain $X$; and

4. judge that all elements in $X$ meet $\left|x_i\right|\le \epsilon$ (i = 0, 1, 2 … 2n). If all elements meet the conditions, the obtained result is the solution at time $t_0+\Delta \mathit{t}$; if not, use the obtained value as a new iterative value and repeat steps (2)–(4).

3. Water quality model

3.1 Control equations

When the 1-D water quality of the main canal of the MRP is described by reference to a QUAL-II comprehensive water quality model, the following equilibrium-dispersion mass equation [10]:

$$\frac{\partial \mathit{AC}}{\partial t}=\frac{\partial }{\partial x}\left(\mathit{EA}\frac{\partial C}{\partial x}\right)-\frac{\partial \mathit{QC}}{\partial x}+\mathit{AS}$$E87

where C is the average concentration of pollutant cross section; A is the cross section area; Q is the average flow of cross section; E is the longitudinal dispersion coefficient; $S$ is the source and drain items.

3.2 Correlation between all water quality variables

When the water quality model for the main route project of the MRP is established according to the actual situation and needs of the MRP and by reference to the QUAL-II comprehensive water quality model [10], the following nine types of water quality variables are considered: (1) dissolved oxygen (DO); (2) ammonia nitrogen; (3) nitrite nitrogen; (4) nitrate nitrogen; (5) biochemical oxygen demand (BOD); (6) chlorophyll-$\alpha$; (7) dissoluble phosphorus; (8) a kind of supposed degradable substance; and (9) a kind of supposed nondegradable substance. One or more water quality variables may be chosen from the model for simulation according to needs. The relation between all variables is shown in Figure 6 given below.

Figure 6.

Diagram of the relation between all water quality variables.

The arrows on the figure are used to show the relation between seven types of water quality variables except for supposed degradable and nondegradable substances. The interaction relation between all water quality variables centers around the content of dissolved oxygen (DO) in water body; a series of chemical reactions occur with physical processes including sediment and the like at the same time. The meaning of all processes is as follows: (1) reoxygenation by water body; (2) consumption of oxygen by sediment; (3) carbonized BOD (CBOD); (4) production of oxygen by organisms including algae and the like via photosynthesis; (5) consumption of oxygen during the process of ammonia-nitrogen oxidation; (6) consumption of oxygen during the process of nitrous acid nitrogen oxidation; (7) precipitation of CBOD; (8) absorption of nitrate nitrogen by organisms including algae and the like; (9) production of dissoluble phosphorus by organisms including algae and the like via respiration action; (10) absorption of dissoluble phosphorus by organisms including algae and the like; (11) death and sediment of organisms including algae and the like; (12) production of ammonia nitrogen by organisms including algae and the like via respiration action; (13) release of ammonia nitrogen by sediment; (14) conversion of ammonia nitrogen via oxidation into nitrite nitrogen; (15) conversion of nitrite nitrogen via oxidation into nitrate nitrogen; and (16) release of dissoluble phosphorus by sediment.

The general main canal is an artificial open channel, the MRP sets high requirements for water quality, and there is no sediment in the canal. Therefore, the release of water quality variables (13 and 16) by sediment and the process of oxygen consumption by sediment (2) in the figure above are not considered when the water quality model is established for the general main canal.

For the aforesaid nine types of water quality variables, they have the same form of transport equations and the difference only in the item of growth or attenuation ($\frac{\mathit{dC}}{\mathit{dt}}$) of chemical reaction; then, for each different water quality variable, its item of growth or attenuation of chemical reaction is shown below:

3.2.1 Chlorophyll-$\alpha$

The reaction process of chlorophyll-$\alpha$ is indicated as follows:

$$\frac{{\mathit{dC}}_A}{\mathit{dt}}=\mu C_A-{\rho }_A{C}_A-\frac{{\sigma }_1}{H}{C}_A$$E88

where $\mu$ is algae growth rate; ${\rho }_A$ is the algae respiratory rate; ${\sigma }_1$ is the algae settling velocity; and $H$ is the average water depth.

3.2.2 Nitrogen cycle

Considering that there are normally three different forms of nitrogen in water body: ammonia nitrogen, nitrite nitrogen, and nitrate nitrogen. Ammonia nitrogen can be oxidized to nitrite nitrogen and transformed to nitrate nitrogen after oxidization. The reaction term of these three substances is shown below, respectively:

Ammonia nitrogen:

$$\frac{{\mathit{dC}}_{N1}}{\mathit{dt}}={\alpha }_1{\rho }_A{C}_A-K_{N1}{C}_{N1}+\frac{{\sigma }_3}{A}$$E89

where ${\alpha }_1$ is the ratio of ammonia nitrogen in algae biomass; ${\sigma }_3$ is the release rate of ammonia nitrogen from underwater wildlife which is not considered in the general main canal; $K_{N1}$ is the oxidization rate of ammonia nitrogen; and other symbols have the same meaning as before.

Nitrite nitrogen:

$$\frac{{\mathit{dC}}_{N2}}{\mathit{dt}}=K_{N1}{C}_{N1}-K_{N2}{C}_{N2}$$E90

where $K_{N2}$ is the oxidization rate and other symbols have the same meaning as before.

Nitrate nitrogen:

$$\frac{{\mathit{dC}}_{N3}}{\mathit{dt}}=K_{N2}{C}_{N2}-{\alpha }_1\mu C_A$$E91

The symbols have the same meaning as before.

3.2.3 Phosphorus cycle

Unlike the nitrogen cycle, the phosphorus cycle is not very complex. Therefore, only the correlation between dissoluble phosphorus and algae as well as phosphorus released by sediment is considered. The item of sediment is not considered in the main route project of the MRP. The reaction item is shown below:

$$\frac{{\mathit{dC}}_p}{\mathit{dt}}={\alpha }_2{\rho }_A{C}_A-{\alpha }_2\mu C_A+\frac{{\sigma }_2}{A}$$E92

where ${\alpha }_2$ is the ratio of phosphorus in algae biomass and ${\sigma }_2$ is the release rate of phosphorus from sediment which is not considered in the main route project.

3.2.4 Carbonized BOD (CBOD)

For the change rate of carbonized BOD, refer to the first-order reaction kinetics, and the reaction item is shown below:

$$\frac{\mathit{dL}}{\mathit{dt}}=-K_{1}L-K_{3}L$$E93

where $K_1$ is the degradation rate of carbonized BOD and $K_3$ is the rate of carbonized BOD consumption due to precipitation. During the change process of carbonized BOD, dissolved oxygen (DO) is consumed only for degradation than for sediment.

3.2.5 Dissolved oxygen (DO)

The reaction item of DO is shown below:

$$\frac{\mathit{dO}}{\mathit{dt}}=K_2\left(O_S-O\right)+\left({\alpha }_3\mu -{\alpha }_4{\rho }_A\right)C_A-K_{1}L-{\alpha }_5{K}_{N1}{C}_{N1}-{\alpha }_6{K}_{N2}{C}_{N2}-\frac{K_4}{A}$$E94

where $O_S$ is the saturation concentration of DO; $O$ is the concentration of DO; ${\alpha }_3$ is the rate of oxygen production by algae per unit via photosynthesis; ${\alpha }_4$ is the rate of oxygen consumption by algae per unit via respiration action; ${\alpha }_5$ is the rate of oxygen consumption in the oxidization of ammonia nitrogen per unit; ${\alpha }_6$ is the rate of oxygen consumption in the oxidization of nitrite nitrogen per unit; $K_2$ is the reoxygenation coefficient; $K_4$ is the coefficient of oxygen consumption by sediment which is not considered in the main route project; and other symbols have the same meaning as before.

3.2.6 Degradable and nondegradable substances

The reaction item of a degradable substance randomly chosen is shown below:

$$\frac{{\mathit{dC}}_R}{\mathit{dt}}=-K_6{C}_R$$E95

where $K_6$ is the degradation rate of the degradable substance. When it is 0, the reaction equation corresponding to a nondegradable substance randomly chosen is obtained.

3.2.7 External source items of water quality

For the external source items in the model, the following two aspects are mainly considered: (1) the flow of water out of a channel from a dividing gate (release sluice) and (2) input of water quality variable point sources.

$$S_{\mathit{ext}}=\mathit{qC}+S_C$$E96

where $q$ is the flow of per unit along the channel; $C$ is the concentration of quality variables of water flowing out of the channel; and $S_C$ is the concentration of water quality variables at a point source.

3.2.8 Basic equation for the model

Substitute the reaction item of the aforesaid different water quality variables into the basic control equation for water quality, and then, we obtain the basic equation for the model:

$$\left\{\begin{array}{l}\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}=\frac{\partial }{\partial x}\left(E\frac{\partial C}{\partial x}\right)+\frac{\mathit{qC}}{A}+S_C\\ \frac{\partial C_1}{\partial t}+u\frac{\partial C_1}{\partial x}=\frac{\partial }{\partial x}\left(E\frac{\partial C_1}{\partial x}\right)+\frac{{\mathit{qC}}_1}{A}+S_{C_1}-K_6{C}_1\\ \frac{\partial L}{\partial t}+u\frac{\partial L}{\partial x}=\frac{\partial }{\partial x}\left(E\frac{\partial L}{\partial x}\right)+\frac{\mathit{qL}}{A}+S_L-K_{1}L-K_{3}L\\ \frac{\partial C_A}{\partial t}+u\frac{\partial C_A}{\partial x}=\frac{\partial }{\partial x}\left(E\frac{\partial C_A}{\partial x}\right)+\frac{{\mathit{qC}}_A}{A}+S_{C_A}+\mu C_A-{\rho }_A{C}_A-\frac{{\sigma }_1}{H}{C}_A\\ \frac{\partial C_P}{\partial t}+u\frac{\partial C_P}{\partial x}=\frac{\partial }{\partial x}\left(E\frac{\partial C_P}{\partial x}\right)+\frac{{\mathit{qC}}_P}{A}+S_{C_P}+{\alpha }_2{\rho }_A{C}_A-{\alpha }_2\mu C_A\\ \frac{\partial C_{N1}}{\partial t}+u\frac{\partial C_{N1}}{\partial x}=\frac{\partial }{\partial x}\left(E\frac{\partial C_{N1}}{\partial x}\right)+\frac{{\mathit{qC}}_{N1}}{A}+S_{C_{N1}}+{\alpha }_1{\rho }_A{C}_A-K_{N1}{C}_{N1}\\ \frac{\partial C_{N2}}{\partial t}+u\frac{\partial C_{N2}}{\partial x}=\frac{\partial }{\partial x}\left(E\frac{\partial C_{N2}}{\partial x}\right)+\frac{{\mathit{qC}}_{N2}}{A}+S_{C_{N2}}+K_{N1}{C}_{N1}-K_{N2}{C}_{N2}\\ \frac{\partial C_{N3}}{\partial t}+u\frac{\partial C_{N3}}{\partial x}=\frac{\partial }{\partial x}\left(E\frac{\partial C_{N3}}{\partial x}\right)+\frac{{\mathit{qC}}_{N3}}{A}+S_{C_{N3}}+K_{N2}{C}_{N2}-{\alpha }_1\mu C_A\\ \frac{\partial O}{\partial t}+u\frac{\partial O}{\partial x}=\frac{\partial }{\partial x}\left(E\frac{\partial O}{\partial x}\right)+\frac{\mathit{qO}}{A}+S_O+K_2\left(O_S-O\right)+\left({\alpha }_3\mu -{\alpha }_4{\rho }_A\right)C_A-K_{1}L\\ -{\alpha }_5{K}_{N1}{C}_{N1}-{\alpha }_6{K}_{N2}{C}_{N2}\end{array}\right.$$E97

where $C,C_1,L,C_A,C_P,C_{N1},C_{N2},C_{N3},O$ are the concentration of degradable substances, nondegradable substances, BOD, chlorophyll-$\alpha$, dissoluble phosphorus, ammonia nitrogen, nitrite nitrogen, nitrate nitrogen, and DO, respectively; $S_C,S_{C_1},S_L,S_{C_A},S_{C_P},S_{C_{N1}},S_{C_{N2}},S_{C_{N3}},S_O$ are the concentration of the aforesaid nine types of water quality variables at a point source, respectively; E is the longitudinal dispersion coefficient; and other parameters refer to all subsections above.

3.3 Discrete equation derivation

Without considering the boundary of a channel or river, a water quality equation can be derived by using the Law of Conservation of Mass in an equilibrium area, and then, a discrete equation form can be obtained. In the equilibrium area, substances are fully mixed and the concentration at all points is uniform, as shown below:

The shaded area in Figure 7 is an equilibrium area. From the figure, the volume of the equilibrium area is:

$$V_j=\frac{1}{4}\left(A_{j-1/2}+A_j\right)\mathrm{\Delta }{x}_{j-1}+\frac{1}{4}\left(A_{j+1/2}+A_j\right)\mathrm{\Delta }{x}_j$$E98

Figure 7.

Schematic diagram of equilibrium area of random nonboundary cross section of a channel.

where $V_j$ is the volume of the equilibrium area at cross section J; $A_{j-\frac{1}{2}},A_{j+\frac{1}{2}}$ are the areas of the cross sections of the inlet and outlet in the equilibrium area, respectively ($A_{j-\frac{1}{2}}=\frac{A_{j-1}+A_j}{2}$, $A_{j+\frac{1}{2}}=\frac{A_j+A_{j+1}}{2}$). Then, the volume of the equilibrium area is:

$$V_j=\frac{1}{8}\left(A_{j-1/2}+3A_j\right)\mathrm{\Delta }{x}_{j-1}+\frac{1}{8}\left(A_{j+1/2}+3A_j\right)\mathrm{\Delta }{x}_j$$E99

Then, the mass variation of water quality variables in the equilibrium area at $\Delta \mathit{t}$ time interval:

$$\mathrm{\Delta }m=V_j^{i+1}{C}_j^{i+1}-V_j^i{C}_j^i$$E100

where i and i + 1 stand for layer i and layer i + 1, respectively.

The change of the concentration of water quality variables in a water body is a comprehensive physical, chemical, and biological processes that are very complicated. The physical process includes advection, molecular diffusion, turbulent diffusion, dispersion, absorption and desorption, sedimentation and resuspension, etc.; and the chemical and biological processes include oxidation, respiration action, photosynthesis, etc. Besides, the concentration of water quality variables in the water body might be affected by external sources such as lateral inflow, sudden point source, etc.

3.3.1 The effect of advection and dispersion on the water quality variables in an equilibrium area

3.3.1.1 Advection

The advection flux $F_x$ of water quality variables at some point in the x direction is:

$$F_x=\mathit{uC}$$E101

where $u$ is the time average velocity of some point in the x direction and $C$ is the time average concentration of water quality variables at some point.

Then, the mass of water quality variables that passes some cross section at $\Delta \mathit{t}$ time interval due to advection is:

$$m={\mathit{AF}}_x\Delta \mathit{t}=A\overline{u}\overline{C}\Delta \mathit{t}=\overline{Q}\overline{C}\Delta \mathit{t}$$E102

where the variables with “—” are the average value of the cross section with the superscript omitted below.

Then, consider that the masses of water quality variables that flow into and out of an equilibrium area by means of advection at $\Delta \mathit{t}$ time interval are, respectively, as follows:

$$m_{11}=Q_{j-\frac{1}{2}}{C}_{j-\frac{1}{2}}\Delta \mathit{t}$$E103

$$m_{12}=Q_{j+\frac{1}{2}}{C}_{j+\frac{1}{2}}\Delta \mathit{t}$$E104

where $Q_{j-\frac{1}{2}}=\frac{Q_{j-1}+Q_j}{2}$, $Q_{j+\frac{1}{2}}=\frac{Q_j+Q_{j+1}}{2}$; $C_{j-\frac{1}{2}}=\theta C_{j-1}+\left(1-\theta \right)C_j$, $C_{j+\frac{1}{2}}=\theta C_j+\left(1-\theta \right)C_{j+1}$; $\theta$ is the upwind factor ($0\le \theta \le 1$). If $Q_{j-\frac{1}{2}}>0$, $\theta \ge \frac{1}{2}$;

If $Q_{j-\frac{1}{2}}\le 0$, $\theta \le \frac{1}{2}$; if $Q_{j-\frac{1}{2}}>0$, $\theta =1$; then this case is called a complete upwind scheme, and the concentration of water quality variables flowing into the equilibrium area under such case is the concentration at an upstream inflow node.

3.3.1.2 Molecular diffusion

Molecular diffusion conforms to Fick’s first law:

$$M_m=-E_m\frac{\partial C}{\partial x}$$E105

where $M_m$ is the molecule diffusive flux at some point in the x direction and $E_m$ is the molecular diffusion coefficient and its range of values generally is ${10}^{-9}\sim {10}^{-8}{m}^2/s$.

Then, consider that the masses of water quality variables that flow into and out of an equilibrium area by means of molecular diffusion at $\Delta \mathit{t}$ time interval are, respectively, as follows:

$$m_{21}=-E_{m,j-\frac{1}{2}}{A}_{j-\frac{1}{2}}\frac{C_j-C_{j-1}}{\Delta x_{j-1}}\Delta \mathit{t}$$E106

$$m_{22}=-E_{m,j+\frac{1}{2}}{A}_{j+\frac{1}{2}}\frac{C_{j+1}-C_j}{\Delta x_j}\Delta \mathit{t}$$E107

where $E_{m,j-\frac{1}{2}}$ and $E_{m,j+\frac{1}{2}}$ are the molecular diffusion coefficients at $x_{j-\frac{1}{2}}$ and $x_{j+\frac{1}{2}}$, respectively; $E_{m,j-\frac{1}{2}}=\frac{E_{m,j-1}+E_{m,j}}{2}$; and $E_{m,j+\frac{1}{2}}=\frac{E_{m,j}+E_{m,j+1}}{2}$.

3.3.1.3 Turbulent diffusion

Like the molecular diffusion, turbulent diffusion also can be expressed by using the form of Fick’s first law:

$$M_t=-E_{\mathit{tx}}\frac{\partial C}{\partial x}$$E108

where $M_t$ is the turbulent diffusion flux at some point in the x direction; $E_{\mathit{tx}}$ is the turbulent diffusion coefficient; for a turbulent flow field whose Reynolds number is about $Re={10}^4$, the turbulent diffusion coefficient is approximately 3.36 × 10⁻⁴ m²/s.

Then, consider that the masses of water quality variables that flow into and out of an equilibrium area by means of turbulent diffusion at $\Delta \mathit{t}$ time interval are, respectively, as follows:

$$m_{31}=-E_{\mathrm{tx},j-\frac{1}{2}}{A}_{j-\frac{1}{2}}\frac{C_j-C_{j-1}}{\Delta x_{j-1}}\Delta \mathit{t}$$E109

$$m_{32}=-E_{\mathrm{tx},j+\frac{1}{2}}{A}_{j+\frac{1}{2}}\frac{C_{j+1}-C_j}{\Delta x_j}\Delta \mathit{t}$$E110

where $E_{\mathrm{tx},j-\frac{1}{2}}$ and $E_{\mathrm{tx},j+\frac{1}{2}}$ are the turbulent diffusion coefficients at $x_{j-\frac{1}{2}}$ and $x_{j+\frac{1}{2}}$, respectively; $E_{\mathrm{tx},j-\frac{1}{2}}=\frac{E_{\mathrm{tx},j-1}+E_{\mathrm{tx},j}}{2}$; and $E_{\mathrm{tx},j+\frac{1}{2}}=\frac{E_{\mathrm{tx},j}+E_{\mathrm{tx},j+1}}{2}$.

3.3.1.4 Dispersion

The expression of dispersion is as follows:

$$M_d=-E_d\frac{\partial C}{\partial x}$$E111

where $M_d$ is the discrete transport flux at some point in the x direction and $E_d$ is the dispersion coefficient with a range of values of 10∼10³ m²/s.

Then, consider that the masses of pollutants that flow into and out of an equilibrium area by means of dispersion at $\Delta \mathit{t}$ time interval are, respectively, as follows:

$$m_{41}=-E_{\mathrm{d},j-\frac{1}{2}}{A}_{j-\frac{1}{2}}\frac{C_j-C_{j-1}}{\Delta x_{j-1}}\Delta \mathit{t}$$E112

$$m_{42}=-E_{\mathrm{d},j+\frac{1}{2}}{A}_{j+\frac{1}{2}}\frac{C_{j+1}-C_j}{\Delta x_j}\Delta \mathit{t}$$E113

where $E_{\mathrm{d},j-\frac{1}{2}}$ and $E_{\mathrm{d},j+\frac{1}{2}}$ are the dispersion coefficients at $x_{j-\frac{1}{2}}$ and $x_{j+\frac{1}{2}}$, respectively; $E_{\mathrm{d},j-\frac{1}{2}}=\frac{E_{\mathrm{d},j-1}+E_{\mathrm{d},j}}{2}$; and $E_{\mathrm{d},j+\frac{1}{2}}=\frac{E_{\mathrm{d},j}+E_{\mathrm{d},j+1}}{2}$.

3.3.1.5 The effect of absorption and desorption and of sedimentation and resuspension on the water quality variables in an equilibrium area

Because the water body in the Danjiangkou Reservoir within the water source region is clear, contains less sand, and has high requirements for water quality, the content of solid matters in the channel such as sediment and the like is extremely low. Therefore, the action of absorption and desorption and of sedimentation and resuspension arising therefrom is omitted during water quality simulation.

3.3.1.6 The effect of advection and dispersion on the water quality variables in an equilibrium area

From the above, we learn that the diffusion coefficient is far greater than the molecular diffusion coefficient and the turbulent diffusion coefficient. Therefore, when a 1-D water quality problem of a channel is solved, advection and dispersion are considered only, with molecular diffusion and turbulent diffusion often omitted. Then, the mass change of water quality variables in an equilibrium area caused by advection and dispersion is:

$$\begin{array}{ll}\mathrm{\Delta }{m}_1& =m_{11}-m_{12}+m_{41}-m_{42}\\ & =\Delta \mathit{t}\left\{C_{j-1}\left(\theta Q_{j-1/2}+\frac{E_{d,j-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}\right)\right.\\ & +C_j\left[\left(1-\theta \right)Q_{j-1/2}-\theta Q_{j+1/2}-\frac{E_{d,j-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}-\frac{E_{d,j+1/2}{A}_{j+1/2}}{\Delta x_j}\right]\\ & +\left.C_{j+1}\left[-\left(1-\theta \right)Q_{j+1/2}+\frac{E_{d,j+1/2}{A}_{j+1/2}}{\Delta x_j}\right]\right\}\end{array}$$E114

3.3.2 The effect of a source-sink term on the water quality variables in an equilibrium area

3.3.2.1 The effect of water quality variable degradation and the correlation between water quality variables on the water quality variables in an equilibrium area

According to section 3.2 correlation between all water quality variables conforms to the first-order reaction kinetics equation. Therefore, when the water quality variable degradation and the correlation between water quality variables are solved, zero- and first-order rates are used to characterize, and the coefficients corresponding to them are ${\alpha }_0$ and ${\alpha }_1$, respectively; at the same time, “+” (plus sign) means reduction of water quality variables, and then, consider the concentration increment of water quality variables in an equilibrium area because of water quality variable degradation and the correlation between water quality variables at $\Delta \mathit{t}$ time interval is:

$$\mathrm{\Delta }{m}_2=\Delta {\mathit{tV}}_j\left(-{\alpha }_0-{\alpha }_1{C}_j\right)$$E115

3.3.2.2 The effect of lateral inflow and a sudden point source on the water quality variables in an equilibrium area

When lateral inflow containing corresponding water quality variables exists along the route of an open channel, or when a sudden point source in the open channel is input, this will bring effect to water quality variables in the open channel. Assuming that the lateral inflow per unit is $q$, the water quality concentration corresponding thereto is $C_q$, and the mass of water quality variables input at a sudden point source per unit time is $\mathit{mm}$, then, consider the concentration change of water quality variables in an equilibrium area owning to lateral inflow or the sudden point source at $\Delta \mathit{t}$ time interval is:

$$\mathrm{\Delta }{m}_3=\Delta {\mathit{tC}}_{q}q\left(\frac{\Delta x_{j-1}}{2}+\frac{\Delta x_j}{2}\right)+\mathit{\text{mm}}\Delta t$$E116

3.3.3 Final form of discrete equation

Considering the conservation of mass in an equilibrium area, the mass change of water quality variables in an equilibrium area at $\Delta \mathit{t}$ time interval is equal to the sum of the mass of water quality variables flowing into and out of the equilibrium area at such time interval.

$$\mathrm{\Delta }m=\mathrm{\Delta }{m}_1+\mathrm{\Delta }{m}_2+\mathrm{\Delta }{m}_3$$E117

Substitute the formulas above into this formula, and then, we obtain:

$$\begin{array}{ll}{V}_j{C}_j-V_j^i{C}_j^i& =\Delta \mathit{t}\left\{C_{j-1}\left(\theta Q_{j-1/2}+\frac{E_{d,j-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}\right)\right.\\ & +C_j\left[\left(1-\theta \right)Q_{j-1/2}-\theta Q_{j+1/2}-\frac{E_{d,j-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}-\frac{E_{d,j+1/2}{A}_{j+1/2}}{\Delta x_j}\right]\\ & +\left.C_{j+1}\left[-\left(1-\theta \right)Q_{j+1/2}+\frac{E_{d,j+1/2}{A}_{j+1/2}}{\Delta x_j}\right]\right\}+\Delta {\mathit{tV}}_j\left(-{\alpha }_0-{\alpha }_1{C}_j\right)\\ & +\Delta {\mathit{tC}}_{q}q\left(\frac{\Delta x_{j-1}}{2}+\frac{\Delta x_j}{2}\right)+\mathit{\text{mm}}\Delta \mathit{\text{t}}\end{array}$$E118

The items without a superscript are the variable values on layer i + 1.

After arrangement, we obtain:

$$\begin{array}{l}{C}_{j-1}\left(\theta Q_{j-1/2}+\frac{E_{d,j-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}\right)+\\ C_j\left[\left(1-\theta \right)Q_{j-1/2}-\theta Q_{j+1/2}-\frac{E_{d,j-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}-\frac{E_{d,j+1/2}{A}_{j+1/2}}{\Delta x_j}-{\alpha }_1{V}_j-\frac{V_j}{\Delta }\right]\\ +C_{j+1}\left[-\left(1-\theta \right)Q_{j+1/2}+\frac{E_{d,j+1/2}{A}_{j+1/2}}{\Delta x_j}\right]=\\ -\frac{V_j^i}{\Delta }{C}_j^i+{\alpha }_0{V}_j-C_{q}q\left(\frac{\Delta x_{j-1}}{2}+\frac{\Delta x_j}{2}\right)-\mathit{mm}\end{array}$$E119

By reference to the form of Formulas (17) and (28), the formula above can be transformed to:

$$A_j{C}_{j-1}+B_j{C}_j+D_j{C}_{j+1}=E_j$$E120

where

$$A_j=\theta Q_{j-1/2}+\frac{E_{d,j-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}$$E121

$$B_j=\left(1-\theta \right)Q_{j-1/2}-\theta Q_{j+1/2}-\frac{E_{d,j-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}-\frac{E_{d,j+1/2}{A}_{j+1/2}}{\Delta x_j}-{\alpha }_1{V}_j-\frac{V_j}{\Delta }$$E122

$$D_j=-\left(1-\theta \right)Q_{j+1/2}+\frac{E_{d,j+1/2}{A}_{j+1/2}}{\Delta x_j}$$E123

$$E_j=-\frac{V_j^i}{\Delta }{C}_j^i+{\alpha }_0{V}_j-C_{q}q\left(\frac{\Delta x_{j-1}}{2}+\frac{\Delta x_j}{2}\right)-\mathit{mm}$$E124

3.4 Boundary conditions

Like considering boundary conditions during hydrodynamic calculation, boundary conditions also can be divided into external and internal boundary conditions according to the actual situation of the general main canal when water quality is simulated. The external boundary conditions are mainly about the upstream and downstream channel system; while the internal boundary conditions mainly involve the effect of hydraulic structures in the channel system on water quality variables during the operation of such structures. Only the dividing gate is discussed here according to the actual situation of the general main canal.

3.4.1 Generalization of channels of the middle route

As described above, there are multiple types of structures along the route of the MRP, including 64 check gates, 88 dividing gates, release sluices, inverted siphons, aqueducts, culverts, etc., and those structures are combined in series, with a close hydraulic relation among upstream and downstream structures. Therefore, generalization should be performed in combination of the features of all structures along the route.

Among many structures (excluding lateral outflow such as dividing gate, release sluice, and the like) along the route, other structures do not cause the water body in the channel to flow out of the channel, and at the same time, there is almost no possibility of lateral inflow under the MRP. Therefore, other structures (excluding lateral outflow such as dividing gate, release sluice, and the like) can be generalized to channels during the generalization. Then, by reference to the classification of nodes in the QUAL—II model, the MRP totally includes the following types of nodes: (1) source node; (2) normal channel node; (3) the node containing a point source; (4) upstream lateral outflow node; (5) downstream lateral outflow node; and (6) the node at a channel end. The information of generalization is shown in Figure 8 below:

Figure 8.

(a) Channel generalization diagram and (b) node type distribution diagram.

3.4.2 Outer boundary conditions

3.4.2.1 Upstream boundary conditions

When upstream boundary conditions are considered, this mainly targets canal head cross sections and their equilibrium area, as indicated by the gray area in Figure 9 above.

Figure 9.

Schematic diagram for upstream boundary conditions.

3.4.2.1.1 Upstream boundary conditions with a known upstream concentration

The upstream concentration is considered to change with time, namely:

$$C=C_1\left(t\right)$$E125

From the formula above, we obtain the concentration $C_{1,0}$ at a canal head cross section $x_1$ at any time; after transforming it into the form of the above formula (120), we obtain:

$$B_1{C}_1+D_1{C}_2=E_1$$E126

where $B_1=1.0$, $D_1=0$, and $E_1=C_{1,0}$.

3.4.2.1.2 Upstream boundary conditions with a known upstream water quality variable flux

A known flux means that the mass of water quality variables passing through cross section $x_1$ per unit time is known, i.e.,

$$F_1=Q_{1}C$$E127

Then, according to the conservation of mass in an equilibrium area and by reference to the process of derivation of the discrete equation in the last section, we know:

$$\begin{array}{l}{C}_1\left[-\theta Q_{1+1/2}-\frac{E_{d,1+1/2}{A}_{1+1/2}}{\Delta x_1}-{\alpha }_1{V}_1-\frac{V_1}{\Delta }\right]+C_2\left[-\left(1-\theta \right)Q_{1+1/2}+\frac{E_{d,1+1/2}{A}_{1+1/2}}{\Delta x_1}\right]\\ =-Q_{1}C-\frac{V_1^i}{\Delta }{C}_1^i+{\alpha }_0{V}_1-\frac{\Delta x_1}{2}{C}_{q}q-\mathit{mm}\end{array}$$E128

After transforming it into the form of the above formula (120), we obtain:

$$B_1{C}_1+D_1{C}_2=E_1$$E129

where $B_1=-\theta Q_{1+1/2}-\frac{E_{d,1+1/2}{A}_{1+1/2}}{\Delta x_1}-{\alpha }_1{V}_1-\frac{V_1}{\Delta }$, $D_1=-\left(1-\theta \right)Q_{1+1/2}+\frac{E_{d,1+1/2}{A}_{1+1/2}}{\Delta x_1}$, and $E_1=-Q_{1}C-\frac{V_1^i}{\Delta }{C}_1^i+{\alpha }_0{V}_1-\frac{\Delta x_1}{2}{C}_{q}q-\mathit{mm}$; all the other parameters have the same meaning as before.

3.4.2.2 Downstream boundary conditions

Like the upstream boundary conditions, when downstream boundary conditions are considered, this mainly targets canal end cross sections and their equilibrium area, as indicated by the gray area in Figure 10 below.

Figure 10.

Schematic diagram for downstream boundary conditions.

3.4.2.2.1 Downstream boundary conditions with a known downstream concentration

The downstream concentration is also considered to change with time, namely:

$$C=C_m\left(t\right)$$E130

From the formula above, we obtain the concentration $C_{m,0}$ at a canal end cross section $x_m$ at any time. After transforming it into the form of the above formula (120), we obtain:

$$A_m{C}_{m-1}+B_m{C}_m=E_m$$E131

where $A_m=0$, $B_m=1.0$, and $E_m=C_{m,0}$.

3.4.2.2.2 Downstream boundary conditions with a known downstream water quality variable flux

For a canal, the concentration of water quality variables in the water body passing through the canal end cross section is the concentration of external water quality variables. Therefore, the mass of water quality variables passing through cross section $x_m$ per unit time is:

$$F_m=Q_m{C}_m$$E132

Then, according to the conservation of mass in an equilibrium area and by reference to the process of derivation of the discrete equation in the last section, we know:

$$\begin{array}{l}{C}_{m-1}\left(\theta Q_{m-1/2}+\frac{E_{d,\mathrm{m}-1/2}{A}_{m-1/2}}{\Delta x_{m-1}}\right)+\\ C_m\left[\left(1-\theta \right)Q_{m-1/2}-Q_m-\frac{E_{d,\mathrm{m}-1/2}{A}_{m-1/2}}{\Delta x_{m-1}}-{\alpha }_1{V}_m-\frac{V_m}{\Delta }\right]\\ =-\frac{V_m^i}{\Delta }{C}_m^i+{\alpha }_0{V}_m-C_{q}q\frac{\Delta x_{m-1}}{2}-\mathit{mm}\end{array}$$E133

After transforming it into the form of the above formula (120), we obtain:

$$A_m{C}_{m-1}+B_m{C}_m=E_m$$E134

where $A_m=\theta Q_{m-1/2}+\frac{E_{d,\mathrm{m}-1/2}{A}_{m-1/2}}{\Delta x_{m-1}}$, $B_m=\left(1-\theta \right)Q_{m-1/2}-Q_m-\frac{E_{d,\mathrm{m}-1/2}{A}_{m-1/2}}{\Delta x_{m-1}}-{\alpha }_1{V}_m-\frac{V_m}{\Delta }$, and $E_m=-\frac{V_m^i}{\Delta }{C}_m^i+{\alpha }_0{V}_m-C_{q}q\frac{\Delta x_{m-1}}{2}-\mathit{mm}$; all the other parameters have the same meaning as before.

3.4.3 Inner boundary conditions

As mentioned above, there are various kinds of structures along the main route project under the MRP. Most of them, however, have less impact on the concentration of water quality variables in the water body. Therefore, when water quality is simulated, the change law of water quality variables at the dividing gate is analyzed only.

The dividing gate along the main route is mainly used to lead the water body out of the channel, dividing it to each water-receiving area. If the water body in the channel gets polluted and the group of polluted water passes through the dividing gate, a certain mass of pollutants will surely pass through the dividing gate and go out of the channel, resulting in the change in the total amount of pollutants in the water body. When treating the pollutants at the dividing gate, we herein think that the dividing gate divides the channel into two sections as shown in Figure 11 below, where the upstream cross section of the dividing gate is $x_j$ and the downstream cross-section is $x_{j+1}$.

Figure 11.

Schematic diagram for boundary conditions of dividing gate.

3.4.3.1 Upstream of dividing gate

When the upstream channel section of the dividing gate is treated, the treatment method is similar to external boundary conditions of the known downstream water quality variable flux within the external boundary of the channel; then by imitating the above formula (133), we obtain:

$$\begin{array}{l}{C}_{j-1}\left(\theta Q_{j-1/2}+\frac{E_{d,\mathrm{j}-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}\right)+\\ C_j\left[\left(1-\theta \right)Q_{j-1/2}-Q_j-\frac{E_{d,\mathrm{j}-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}-{\alpha }_1{V}_j-\frac{V_j}{\Delta }\right]\\ =-\frac{V_j^i}{\Delta }{C}_j^i+{\alpha }_0{V}_j-C_{q}q\frac{\Delta x_{j-1}}{2}-\mathit{mm}\end{array}$$E135

After transforming it into the form of the above formula (120), we obtain:

$$A_j{C}_{j-1}+B_j{C}_j=E_j$$E136

where $A_j=\theta Q_{j-1/2}+\frac{E_{d,\mathrm{j}-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}$, $B_j=\left(1-\theta \right)Q_{j-1/2}-Q_j-\frac{E_{d,\mathrm{j}-1/2}{A}_{j-1/2}}{\Delta x_{j-1}}-{\alpha }_1{V}_j-\frac{V_j}{\Delta }$, and $E_j=-\frac{V_j^i}{\Delta }{C}_j^i+{\alpha }_0{V}_j-C_{q}q\frac{\Delta x_{j-1}}{2}-\mathit{mm}$; all the other parameters have the same meaning as before.

3.4.3.2 Downstream of dividing gate

When the downstream channel section of the dividing gate is treated, this may be considered similar to external boundary conditions of the known upstream water quality variable flux within the external boundary of the channel, and the concentration of water quality variables is $C_j$; then by imitating the above formula (128), we obtain:

$$\begin{array}{l}{C}_j{Q}_{j+1}+C_{j+1}\left(-\theta Q_{j+1+1/2}-\frac{V_{j+1}}{\Delta }-\frac{E_{d,\mathrm{j}+1+1/2}{A}_{j+1+1/2}}{\Delta x_{j+1}}-{\alpha }_1{V}_{j+1}\right)+\\ C_{j+2}\left[-\left(1-\theta \right)Q_{j+1+1/2}+\frac{E_{d,\mathrm{j}+1+1/2}{A}_{j+1+1/2}}{\Delta x_{j+1}}\right]=\\ -\frac{V_{j+1}^i}{\Delta }{C}_{j+1}^i+{\alpha }_0{V}_{j+1}-\frac{\Delta x_{j+1}}{2}{C}_{q}q-\mathit{mm}\end{array}$$E137

After transforming it into the form of the above formula (120), we obtain:

$$A_{j+1}{C}_j+B_{j+1}{C}_j+D_{j+1}{C}_{j+2}=E_{j+1}$$E138

where $A_{j+1}=Q_{j+1}$, $B_{j+1}=-\theta Q_{j+1+1/2}-\frac{V_{j+1}}{\Delta }-\frac{E_{d,\mathrm{j}+1+1/2}{A}_{j+1+1/2}}{\Delta x_{j+1}}-{\alpha }_1{V}_{j+1}$, $D_{j+1}=-\left(1-\theta \right)Q_{j+1+1/2}+\frac{E_{d,\mathrm{j}+1+1/2}{A}_{j+1+1/2}}{\Delta x_{j+1}}$, and $E_{j+1}=-\frac{V_{j+1}^i}{\Delta }{C}_{j+1}^i+ {\alpha }_0{V}_{j+1}-\frac{\Delta x_{j+1}}{2}{C}_{q}q-\mathit{mm}$; all the other parameters have the same meaning as before.

3.5 Determination of model coefficient

Because of a complicated transport and diffusion law of water quality variables in a water body, the form of the equation combining coefficient and concentration is often used to describe such law. This also makes determining the coefficient of the equation becomes a very important step to describe the law.

There are many methods to determine a coefficient. The common ones are as follows—Method 1: set up a mathematical model and assume a group of coefficients; compare the result of numerical simulation and the real test result and make adjustment to the coefficients according to the result of such comparison; repeat this procedure a number of times to determine the values of all the coefficients; Method 2 (a site test method): put a tracer into a simulative water body, trace and monitor the tracer and then calculate the coefficients by using test information; Method 3: estimate the coefficients by using experimental formulas.

Owing to the specificity of the MRP, it is difficult to put a tracer and conduct a site test. However, water flew into the canal not long before and water quality monitoring data are incomplete, so determining parameters in the book is mainly based on research of the same time and experimental formulas.

The determination of a longitudinal dispersion coefficient is mainly based on experimental formulas in the book. The experimental formula proposed by Henry Liu is used:

$$E_d=\gamma \frac{u_{\ast }{A}^2}{H^3}$$E139

where $\gamma$ is an empirical coefficient generally having a range of values of 0.5–0.6; $u_{\ast }$ is the frictional velocity; $u_{\ast }=\sqrt{\mathit{gHJ}}$; $J$ is the hydraulic gradient; $A$ is the cross-section area; and $H$ is the cross section water depth.

The values of other coefficients are indicated in Table 1 below:

Pollutant	Reaction coefficient	Value	Unit	Range of values
Nondegradable substances	u0	8	day^−1	—
	u1	0.8	gm^−3 day^−1	—
Algae	μ	0.075	day^−1	1.0–3.0
	${\rho }_A$	0.1	day^−1	0.005–0.5
	σ1	1	mday^−1	0.153–1.83
Nitrogen	ɑ1	0.085	gg^−1	0.08–0.09
	KN1	0.03	day^−1	0.01–0.5
	KN2	1.5	day^−1	0.5–2.0
Phosphorus	ɑ2	0.0135	gg^−1	0.012–0.015
CBOD	K1	0.15	day^−1	0.1–0.2
	K3	0.18	day^−1	—
Oxygen	OS	9.08	gg^−1	—
	ɑ3	1.6	gg^−1	1.4–1.8
	ɑ4	2	gg^−1	1.6–2.3
	ɑ5	3.5	gg^−1	3.0–4.0
	ɑ6	1.07	gg^−1	1.0–1.14
	K2	5	day^−1	0–100

Table 1.

Values of parameters corresponding to all water quality variables of the model.

3.6 Model solving

Assuming that the channel totally has m cross sections to be calculated; excluding the cross sections at the head and end of the channel, there are m−2 cross sections; then, for each type of water quality variable, m−2 equations whose form is similar to that of Formula (120) for the whole channel are set up; and two equations for cross sections at the head and end of the channel by using the method to treat boundary conditions are established; and finally, m equations are obtained. By simultaneously solving such m equations, the concentration of corresponding water quality variables can be obtained.

4. Model calibration and verification

4.1 Hydraulic model

A design value of 0.015 is taken as the roughness coefficient of the channel. The local head loss coefficient is adjusted gradually in a reasonable range of values to minimize the error of the simulation and actually measured values of water level. The operation data actually measured at AM 8:00, December 25, 2015 were selected and input the model for calculation; the first 56 check gates were selected, and the actually measured and simulation values of water level of downstream check gate were compared, and the average error is 3.42 cm (Figure 12). The results indicate that the model was rational and appropriate for simulating the water surface profile and hydrodynamic change process along channels.

Figure 12.

Comparison between simulation value and actually measured value of water level of downstream check gate.

4.2 Water quality model

Because actually measured data cannot be obtained, this model and MIKE11 model [2] are compared with the same calculation example to verify the precision of this model.

Assume that there is a channel with a prismatic trapezium cross section, which is 10 km long and has a bottom width of 67.5 m, a slope of 2.5, a bottom slope of about 0.00015, and a roughness coefficient of 0.027; and that there is a reservoir with a constant flow rate of 2000 m³/s at a place far from the channel upstream. Suppose that 1 ton of soluble but refractory pollutants were leaked suddenly on the right bank of at an upstream cross section of a channel at some time, then try to estimate the process of pollutant concentration change at the downstream outflow cross section after pollution as well as the pollutant concentration distribution in the channel after 0.5, 1, 1.5, and 2 h.

1. Firstly, calculate the following hydraulic elements of the channel:

   1. the water depth $h$ estimated according to the open channel uniform flow formula $Q={\mathit{AR}}^{2/3}\sqrt{i}{n}^{-1}$ = 11.20 m;

   2. the area of flow section $A$ 1069.60 m;

   3. average velocity of cross section $u$ = 1.87 m/s;

   4. hydraulic radius of cross section $R$ = 8.34 m;

   5. frictional velocity of cross section $u^{\ast }=\sqrt{\mathit{gRi}}$ = 0.11 m/s;

   6. estimated channel diffusion coefficient $D_x=D_y=0.15{\mathit{hu}}^{\ast }$ = 0.18 m²/s; and

   7. estimated longitudinal diffusion (including dispersion) coefficient $D=6.01{\mathit{hu}}^{\ast }$ = 7.4 m²/s.

After obtaining the aforesaid hydraulic elements, estimate the pollutant concentration change process at the places (outflow section) 5 and 10 km downstream of the channel (See Figure 13), by using the water quality model and MIKE11 AD module mentioned in 2.3.

Figure 13.

The process of calculating channel pollutant concentration by using a water quality calculation method. (a) The concentration process at the place 5 km downstream of the channel and (b) the concentration process at the place 10 km downstream of the channel.

From the above figure, we see that no matter whether a numerical model under general conditions or MIKE11 is used for simulation, the result of simulation is near to the concentration peak or the time when such peak occurs, well-reflecting transport law of pollutants with water flow in the channel. Moreover, the numerical method applied under this chapter adopts the assumption of equilibrium areas. Under such assumption, if a group of sudden polluted water does not diffuse fully, but gathers in some space, the numerical method will lead to a low simulated concentration; on the other hand, for pollution due to a sudden point source, the treatment technology used in this chapter is on the basis of instantaneous and uniform mixing of concentration and omitted recession process of pollutants in the first channel segment. Therefore, the simulation by the numerical method results in a low-concentration peak.

5. Application of typical project cases

After calibration verification, the aforesaid hydrodynamic water quality model can be applied to contingent operation of the MRP; the application of the hydrodynamic model and water quality model is further verified by taking channel drainage under common operating conditions and sudden water pollution accidents for examples.

5.1 Drainage simulation

In case of a sudden event under the MRP, the release sluice is the main structure for drainage. Therefore, the 36th channel pool is selected for application to analyze the hydraulic response process when the release sluice opens. The basic parameters of the channel pool are as follows:

Assuming that the 36th channel pool undergoes a sudden event, the Anyanghe and Zhanghe check gates are closed within 15 min. Assuming that the canal discharge is 30% of the design discharge, and that the upstream water level of the check gates is the design water level. To analyze the effect arising from different opening speeds of drainage gate, we assume three cases (no opening, opened to 60 m³/s within 15 min, and opened to 120 m³/s within 15 min). The time of simulation is 24 h with a step of 10 min.

After simulation of the drainage process of the 36th channel pool by using the hydraulic model, the change processes of the downstream water level of Anyanghe check gate and the upstream water level of Zhanghe check gate are shown in Figures 14 and 15, respectively. From these figures, we can see that (1) the downstream water level of Anyanghe check gate goes down first and then up, and the water level fluctuation range becomes small gradually and finally steady without opening the release sluice; while the release sluice is opened, the water level fluctuates slightly, but has a trend to become small on the whole until the channel pool is drained; and that (2) the upstream water level of Zhanghe check gate goes up first and then down, and the water level fluctuation range becomes small gradually and finally steady without opening the release sluice; while the release sluice is opened, the water level fluctuates slightly, but has a trend to become small on the whole until the channel pool is drained. Therefore, opening the release sluice makes it easy to reduce the likelihood of high water level in the channel due to fast closing of two check gates in the channel pool.

Figure 14.

Change process of downstream water level of Anyanghe check gate during drainage.

Figure 15.

Change process of upstream water level of Zhanghe check gate during drainage.

5.2 Simulation of sudden water pollution accidents

Since a sudden water pollution event happens generally at some place, assume that water pollution is instantaneously caused by a point source (Tables 2 and 3) and that pollutants are nondegradable. In order to analyze the effect caused by different pollutant amount levels, assume that there are three types of pollutant masses (1, 5, and 10 t). In order to analyze the effect of different channel operating conditions, assume that there are three types of channel pool flows (30, 50, and 70% of design flow) but that the upstream water level of each check gate is design flow. In order to analyze the effect of sudden water pollution at different places, assume that there are five accident places (10, 30, 50, 70, and 90% of channel pool length). Then, we obtain 45 cases, which have a simulation time of 24 h with a step of 10 min.

Channel pool no.	Control structure	Pile no. (km)	Length of channel pool (km)	Elevation of sluice bottom (m)	Design flow (m^3/s)	Design water level (m)
36	Check gate in Anyang River	717.045	14.321	85.60	235	92.67
	Release sluice in Zhanghe River	730.623		85.19	120	/
	Check gate in Zhanghe River	731.366		82.50	235	91.87

Table 2.

Table of basic information on the 36th channel pool.

‘/’ means that there is not Design Water Level of Release sluice in Zhanghe River.

Item	Content	Quantity
Pollutant category	Nondegradable substance	1
Mode of occurrence	An instantaneous point source	1
The channel pool where an accident occurs	The 36th channel pool	2
Pollutant mass	1, 5, and 10t	3
Place of occurrence	10, 30, 50, 70, and 90% of channel pool length	5
Channel pool flow	30, 50, and 70% of design flow	3
Simulation time/step	24 h/10 min	1

Table 3.

Sudden water pollution accident cases.

Considering the water supply features of the MRP, the most concern, after a sudden water pollution accident in the channel, is the water quality state at the dividing gate and the upstream water quality state of the lower check gate. This is because, these have a direct effect on the safety of supplied water quality at such dividing gate and in the downstream areas. Three characteristic parameters, pollutant arrival time (T₀), concentration peak (C_max), and the occurrence time of concentration peak (T_Cmax), are used to analyze the process of pollutant transport and diffusion. The pollutant arrival time is the time at which the concentration of pollutants at a water quality control point (dividing gate or upstream of lower check gate) exceeds 0.001 mg/L, because toxic effect will occur when the concentration of part of substances (for example, Hg and Cd) in a natural water body is above 0.001 mg/L [11].

In combination of the simulation results of the aforesaid cases, we can find that the concentration of the pollutants at each water quality control point goes up first and then down. These are shown by taking the concentration change process of the pollutants upstream of lower Zhanghe River check gate in case that the 36th channel pool undergoes a sudden water pollution accident at the place 10% of channel pool length and with a flow of 30% of design flow for example (Figure 16). The statistics of three characteristic parameters during the concentration change process of pollutants at each water quality control point under each case is made, as shown in Tables 2–4; and contrastive analysis of these three parameters is performed.

Figure 16.

Concentration change process of the pollutants upstream of Zhanghe River check gate in case that the 36th channel pool undergoes a sudden water pollution accident at the place 10% of channel pool length and with a flow of 30% of design flow.

Place of accident	Pollutant mass (t)	30% of design flow			50% of design flow			70% of design flow
		T0 (min)	Cmax (mg/L)	TCmax (min)	T0 (min)	Cmax (mg/L)	TCmax (min)	T0 (min)	Cmax (mg/L)	TCmax (min)
0.1 L	1	100	0.3843	560	60	0.3679	380	40	0.3507	260
	5	70	1.9212	560	40	1.8397	380	30	1.7535	260
	10	60	3.8425	560	40	3.6793	380	20	3.507	260
0.3 L	1	60	0.4107	450	40	0.3939	310	20	0.3771	210
	5	40	2.0535	450	20	1.9693	310	10	1.8856	210
	10	30	4.1071	450	20	3.9385	310	10	3.7711	210
0.5 L	1	20	0.4965	280	10	0.475	200	10	0.4577	130
	5	10	2.4823	280	10	2.3752	200	10	2.2887	130
	10	10	4.9647	280	10	4.7504	200	10	4.5773	130
0.7 L	1	10	0.5387	220	10	0.5146	150	10	0.4975	100
	5	10	2.6935	220	10	2.5729	150	10	2.4873	100
	10	10	5.3871	220	10	5.1457	150	10	4.9745	100
0.9 L	1	10	0.8862	110	10	0.8199	80	10	0.7963	50
	5	10	4.4308	110	10	4.0996	80	10	3.9816	50
	10	10	8.8615	110	10	8.1993	80	10	7.9633	50

Table 4.

Characteristic parameters of concentration change process of the pollutants upstream of Zhanghe River check gate in all sudden water pollution accident cases that the 36th channel pool undergoes.

5.2.1 Pollutant arrival time

When a sudden water pollution accident happens at the same place with the same pollutant mass, a greater flow in the channel pool means that the pollutants arrive at the Zhanghe River check gate at an earlier time. Because the calculation step is 10 min, pollutants under some cases will diffuse fast to the check gate; all the statistical results are 10 min. When a sudden water pollution accident happens at the same place at the same channel pool flow, the greater the pollutant mass is, the earlier the pollutants arrive the check gate. When pollutant mass and channel pool flow are the same, the farther the place of a sudden water pollution accident is from Anyang River check gate, the earlier the pollutants arrive the Zhanghe River check gate.

5.2.2 Concentration peak

When a sudden water pollution accident happens at the same place at the same pollutant mass, the greater the channel pool flow is, the smaller the concentration peak at Zhanghe River check gate is. When a sudden water pollution accident happens at the same place at the same channel pool flow, the greater the pollutant mass is, the greater the concentration peak at Zhanghe River check gate is; and the peak has a proportional relation with the pollutant mass. When pollutant mass and channel pool flow are the same, the farther the place of a sudden water pollution accident is from Anyang River check gate, the greater the concentration peak at Zhanghe River check gate is.

5.2.3 Occurrence time of concentration peak

When a sudden water pollution accident happens at the same place at the same pollutant mass, the greater the channel pool flow is, the earlier the concentration peak at Zhanghe River check gate occurs. When a sudden water pollution accident happens at the same place with the same channel pool flow but with different pollutant mass, the concentration peak at Zhanghe River check gate occurs at the same time. When pollutant mass and channel pool flow are the same, the farther the place of a sudden water pollution accident is from Anyang River check gate, the earlier the concentration peak at Zhanghe River check gate occurs.

6. Summary

For needs predicted via simulation of sudden water pollution accidents in the main canal, this chapter develops the 1-D hydrodynamic water quality coupling simulation model, which can realize the comprehensive simulation of the hydraulic response process under normal operating conditions, emergency operation operating conditions, etc., of the main canal and of nine types of water quality variables. Then, the chapter proves the precision of the hydrodynamic model by using the data actually measured and the precision of the water quality model by using a MIKE11 model. Furthermore, the chapter explains the applicability of the 1-D hydrodynamic water quality coupling simulation model in the operation of sudden water pollution accidents via simulation of channel pool drainage and sudden water pollution accident cases.

However, the water quality variables are not enough, and the 1-D model is less accurate than 2-D or 3-D model. In the future, a hydrodynamic water quality coupling simulation 2-D or 3-D model for the main canal should be developed, and the model can consider other water quality variables such as oil and algae, which can analyze the hydrodynamic and water quality process more precisely.