Open access peer-reviewed compact

Simulation Technology for Hydrodynamic and Water Quality in the Main Canal

By Yu Tian, Hezhen Zheng, Xiaohui Lei and Wei Dai

Reviewed: October 15th 2018Published: November 21st 2018

DOI: 10.5772/intechopen.82022

Downloaded: 167

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.

Keywords

  • hydrodynamic
  • water quality
  • simulation

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:

At+Qx=qlE1
tQA+xβQ22A2+ghx+gSfS0=0E2

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

Sf can be determined according to a discharge modulus:

Sf=QQK2E3

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 θand φare a time weight coefficient and a space weight coefficient, respectively, with the range of values: 0 ≤ θ≤ 1, 0 ≤ φ≤ 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:

fL=θfji+1+1θfjiE4
fR=θfj+1i+1+1θfj+1iE5
fU=φfj+1i+1+1φfji+1E6
fD=φfj+1i+1φfjiE7

Then, the functional value for point M is:

fM=φfR+1φfLE8

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

fM=φθfj+1i+1+1θfj+1i+1φθfji+1+1θfjiE9

Temporal discretization:

ftMfUfDΔt=φfj+1i+1+1φfji+1φfj+1i+1φfjiΔt=φfj+1i+1fj+1iΔt+1φfji+1fjiΔtE10

Spatial discretization:

fxMfRfLΔx=θfj+1i+1+1θfj+1iθfji+1+1θfjiΔx=θfj+1i+1fji+1Δx+1θfj+1ifjiΔxE11

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

Continuity equation:

φΔt(Aj+1i+1Aj+1i)+1φΔt(Aji+1Aji)+θΔx(Qj+1i+1Qji+1)+1θΔx(Qj+1iQji) θ[φqj+1i+1+(1φ)qji+1](1θ)[φqj+1i+(1φ)qji]=0E12

Momentum equation:

φΔt(Qj+1i+1Aj+1i+1Qj+1iAj+1i)+1φΔt(Qji+1Aji+1QjiAji)+θΔx[12(βj+1i+1Qj+1i+1Aj+1i+1)212(βji+1Qji+1Aji+1)2]+1θΔx[βj+1i2(Qj+1iAj+1i)2βji2(QjiAji)2]+θgΔx(hj+1i+1hji+1)+(1θ)gΔx(hj+1ihji)+ θg[Sf,j+1i+1+(1φR)Sf,ji+1]+(1θ)g[Sf,j+1i+(1φR)Sf,ji]=0E13

where φRis the space weight coefficient of hydraulic gradient, with a value different from that of φ.

When the water in a channel flows slowly, the characteristic root symbols of the control equation are opposite, i.e., λ1>0 or λ2< 0; therefore, Cr=λ1ΔtΔx>0 or Cr=λ2ΔtΔx< 0. Because Cr>0 and Cr< 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:

Cr=ughΔtΔx>φ0.50.5θ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., Δhand ΔQ. Firstly, area and flow are transformed into an increment form:

Aji+1=Aj+ΔAj=Aj+BjΔhjE15
Qji+1=Qj+ΔQjE16

where * stands for the value of the variable at the previous cycle step; ΔA, Δh, and ΔQare 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:

ajΔhj+bjΔQj+cjΔhj+1+djΔQj+1=pjE17

where

aj=1φBj/ΔtE18
bj=θ/ΔxE19
cj=φBj+1/ΔtE20
dj=θ/ΔtE21
pj=φΔtAj+1Aj+1i1φΔtAjAjiθΔxQj+1Qj1θΔxQj+1iQji+ θφqj+1i+1+1φqji+1+1θφqj+1i+1φqjiE22

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

Qji+12=Qj2+2QjΔQjE23
1Kji+12=1Kj22Kj3KhjΔhjE24
Sf,ji+1=Sf,j+2QjKj2ΔQj2Sf,jKjKhjΔhjE25
Qji+1Aji+1=QjAj+1AjΔQjQjBjAj2ΔhjE26
Qji+1Aji+12=QjAj2+2QjAj2ΔQj2Qj2BjAj3ΔhjE27

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

aj+1Δhj+bj+1ΔQj+cj+1Δhj+1+dj+1ΔQj+1=pj+1E28

where

aj+1=1φΔtQjBjAj2+θΔxQj2BjAj3θgΔx2θ1φRgSf,jKjKhjE29
bj+1=1φΔt1AjθΔxQjAj2+2θ1φRgQjKj2E30
cj+1=φΔtQj+1Bj+1Aj+12θΔxQj+12Bj+1Aj+13+θgΔx2θφRgSf,j+1Kj+1Khj+1E31
dj+1=φΔt1Aj+1+θΔxQj+1Aj+12+2θφRgQj+1Kj+12E32
pj+1=φΔtQj+1Aj+1Qj+1iAj+1i1φΔtQjAjQjiAjiθΔx12Qj+1Aj+1212QjAj21θΔx12Qj+1iAj+1i212QjiAji2θgΔxhj+1hj1θgΔxhj+1ihjiθgSf,j+11θgSf,jE33

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:

h1=htE34

where h1is the upstream water depth.

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

Δh1=h1i+1h1E35

where h1i+1is the upstream water depth at time (i + 1); h1is 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:

c1Δh1+d1ΔQ1=p1E36

where c1= 1.0, d1= 0, and p1= h1i+1h1.

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:

Q1=QtE37

where Q1is the upstream flow.

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

ΔQ1=Q1i+1Q1E38

where Q1i+1is the upstream flow at time (i + 1) and Q1is 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:

c1Δh1+d1ΔQ1=p1E39

where c1= 0, d1= 1.0, and p1=Q1i+1Q1.

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:

hn=htE40

where hnis the downstream water depth.

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

Δhn=hni+1hnE41

where hni+1is the downstream water depth at time (i + 1); hnis 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:

anΔhn+bnΔQn=pnE42

where an=1.0, bn= 0, and pn= hni+1hn.

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:

Qn=QtE43

where Qnis the upstream flow.

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

ΔQn=Qni+1QnE44

where Qni+1is the upstream flow at time (i + 1); Qnis 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:

anΔhn+bnΔQn=pnE45

where an= 0, bn= 1.0, and pn= Qni+1Qn.

2.4.1.5 Upstream water depth-flow boundary

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

Q1=fh1E46

where h1is the upstream water depth and Q1is 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:

ΔQ1+Q1=fh1+dfdhΔh1E47

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

c1Δh1+d1ΔQ1=p1E48

where c1=dfdhand d1= 1.0; and p1=fh1Q1.

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:

Qn=fhnE49

where hnis the downstream water depth and Qnis 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:

ΔQn+Qn=fhn+dfdhΔhnE50

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

anΔhn+bnΔQn=pnE51

where an=dfdhand bn= 1.0; and pn= fhnQn.

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:

Qj=Qj+1+QfE52
hj+Zj+12gQjAj2=hj+1+Zj+1+12gQj+1Aj+12E53

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; Qjand Qj+1are the upstream and downstream cross section flows, respectively; Qfis the divided water flow at the dividing gate; Zjand Zj+1are the elevation of the channel bottom of upstream and downstream cross sections at the dividing gate.

After linearization of the aforesaid formula, we obtain:

ΔQjΔQj+1=Qj+1Qj+QfE54
1Qj2BjgAj3Δhj+QjgAj2ΔQj1Qj+12Bj+1gAj+13Δhj+1Qj+1gAj+12ΔQj+1=hjZj12gQjAj2+hj+1+Zj+1+12gQj+1Aj+12E55

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

aj=0,bj=1.0,cj=0,dj=1.0,pj=Qj+1Qj+Qf.E56
aj+1=1Qj2BjgAj3,bj+1=QjgAj2,cj+1=1Qj+12Bj+1gAj+13,dj+1=Qj+1gAj+12E57
pj+1=hjZj12gQjAj2+hj+1+Zj+1+12gQj+1Aj+12.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:

Qj=Qj+1E59
hj+Zj+12gQjAj2=hj+1+Zj+1+12gQj+1Aj+12+ζ12gQjAj2Qj+1Aj+12+ξ12gQjAj2E60

Figure 3.

Schematic diagram of transition section.

where ζis a local loss coefficient and ξ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:

ΔQjΔQj+1=Qj+1QjE61
11ζQj2BjgAj3Δhj+1ζQjgAj2ΔQj1ζQj+1gAj+12+ξQj+1gAj+12ΔQj+111ζQj+12Bj+1gAj+13ξQj+1Qj+1Bj+1gAj+13Δhj+1=hjZj1ζ2gQjAj2+hj+1+Zj+1+1ζ2gQj+1Aj+12+ξQj+1Qj+12gAj+12E62
11+ζQj2BjgAj3Δhj+1+ζQjgAj2ΔQj1+ζQj+1gAj+12+ξQj+1gAj+12ΔQj+111+ζQj+12Bj+1gAj+13ξQj+1Qj+1Bj+1gAj+13Δhj+1=hjZj1+ζ2gQjAj2+hj+1+Zj+1+1+ζ2gQj+1Aj+12+ξQj+1Qj+12gAj+12E63

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

aj=0,bj=1.0,cj=0,dj=1.0,pj=Qj+1Qj.E64

The coefficient for the converge section:

aj+1=(1(1ζ)(Qj)2Bjg(Aj)3),bj+1=(1ζ)Qjg(Aj)2,cj+1=((1ζ)Qj+1g(Aj+1)2+ξQj+1g(Aj+1)2),dj+1=(1(1ζ)(Qj+1)2Bj+1g(Aj+1)3ξQj+1Qj+1Bj+1g(Aj+1)3)E65
pj+1=hjZj1ζ2gQjAj2+hj+1+Zj+1+1ζ2gQj+1Aj+12+ξQj+1Qj+12gAj+12.E66

The coefficient for the divergence section:

aj+1=(1(1+ζ)(Qj)2Bjg(Aj)3),bj+1=(1+ζ)Qjg(Aj)2,cj+1=((1+ζ)Qj+1g(Aj+1)2+ξQj+1g(Aj+1)2),dj+1=(1(1+ζ)(Qj+1)2Bj+1g(Aj+1)3ξQj+1Qj+1Bj+1g(Aj+1)3)E67
pj+1=hjZj1+ζ2gQjAj2+hj+1+Zj+1+1+ζ2gQj+1Aj+12+ξQj+1Qj+12gAj+12.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:

Qj=Qj+1E69
hj+Zj+12gQjAj2=hj+1+Zj+1+12gQj+1Aj+12+ζ112gQjAj2+ζ212gQj+1Aj+12+Qj+1Qj+1Kj+12LE70

Figure 4.

Schematic diagram of inverted siphon.

where ζ1and ζ2are 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:

ΔQjΔQj+1=Qj+1QjE71
11ζ1Qj2BjgAj3Δhj1+ζ2Qj+1gAj+12+2Qj+1Kj+12LΔQj+111+ζ2Qj+12Bj+1gAj+13Δhj+1+1ζ1QjgAj2ΔQj=hjZj1ζ12gQjAj2+hj+1+Zj+1+1+ζ22gQj+1Aj+12+Qj+1Qj+1Kj+12LE72

where K=AnR23.

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

aj=0,bj=1.0,cj=0,dj=1.0,pj=Qj+1Qj.E73
aj+1=(1(1ζ1)(Qj)2Bjg(Aj)3),bj+1=(1ζ1)Qjg(Aj)2,cj+1=(1(1+ζ2)(Qj+1)2Bj+1g(Aj+1)3),dj+1=(1(1+ζ2)(Qj+1)2Bj+1g(Aj+1)3)E74
pj+1=hjZj1ζ12gQjAj2+hj+1+Zj+1+1+ζ22gQj+1Aj+12+Qj+1Qj+1Kj+12LE75

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:

Qj=Qj+1E76
Qj=σcσsmbfHjHj+1eE77

Figure 5.

Schematic diagram of arc check gate.

where σcis the side-contract coefficient; σsis the submergence coefficient; mis the flow coefficient; bis the width of sluice; Hjand Hj+1are the downstream water levels; eis the opening of sluice.

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

ΔQjΔQj+1=Qj+1QjE78
σcσsmbfHjΔHjΔQj+σcσsmbfHj+1ΔHj+1=QjσcσsmbfE79

where f=σcσsmbfHjHj+1e.

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

aj=0,bj=1.0,cj=0,dj=1.0,pj=Qj+1Qj.E80
aj+1=σcσsmbfHj,bj+1=1,cj+1=σcσsmbfHj+1,dj+1=0E81
pj+1=QjσcσsmbfE82

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:

a2iΔhi+b2iΔQi+c2iΔhi+1+d2iΔQi+1=D2iE83
ai+1Δhi+b2i+1ΔQi+c2i+1Δhi+1+d2iΔQi+1=D2i+1E84

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:

AX=DE85

where

A=c1d1a2b2c2d2a3b3c3d3a4b4c4d4a5b5c5d5a2n2b2n2c2n2d2n2a2n1b2n1c2n1d2n1a2nb2n,X=Δh1ΔQ1Δh2ΔQ2Δh3ΔQn1ΔhnΔQn,D=D0D1D2D3D4D2n2D2n1D2n.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 t0, the following process can be used to determine the water depth and flow at time t0+Δt:

  1. assume the initial value of the water depth and flow of each cross section of the channel at time t0+Δt;

  2. calculate the value of all elements in matrices Aand D;

  3. solve equation AX=Dand obtain X; and

  4. judge that all elements in Xmeet xiε(i = 0, 1, 2 … 2n). If all elements meet the conditions, the obtained result is the solution at time t0+Δ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]:

ACt=xEACxQCx+ASE87

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; Sis 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-α; (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 (dCdt) 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-α

The reaction process of chlorophyll-αis indicated as follows:

dCAdt=μCAρACAσ1HCAE88

where μis algae growth rate; ρAis the algae respiratory rate; σ1is the algae settling velocity; and His 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:

dCN1dt=α1ρACAKN1CN1+σ3AE89

where α1is the ratio of ammonia nitrogen in algae biomass; σ3is the release rate of ammonia nitrogen from underwater wildlife which is not considered in the general main canal; KN1is the oxidization rate of ammonia nitrogen; and other symbols have the same meaning as before.

Nitrite nitrogen:

dCN2dt=KN1CN1KN2CN2E90

where KN2is the oxidization rate and other symbols have the same meaning as before.

Nitrate nitrogen:

dCN3dt=KN2CN2α1μCAE91

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:

dCpdt=α2ρACAα2μCA+σ2AE92

where α2is the ratio of phosphorus in algae biomass and σ2is 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:

dLdt=K1LK3LE93

where K1is the degradation rate of carbonized BOD and K3is 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:

dOdt=K2OSO+α3μα4ρACAK1Lα5KN1CN1α6KN2CN2K4AE94

where OSis the saturation concentration of DO; Ois the concentration of DO; α3is the rate of oxygen production by algae per unit via photosynthesis; α4is the rate of oxygen consumption by algae per unit via respiration action; α5is the rate of oxygen consumption in the oxidization of ammonia nitrogen per unit; α6is the rate of oxygen consumption in the oxidization of nitrite nitrogen per unit; K2is the reoxygenation coefficient; K4is 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:

dCRdt=K6CRE95

where K6is 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.

Sext=qC+SCE96

where qis the flow of per unit along the channel; Cis the concentration of quality variables of water flowing out of the channel; and SCis 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:

Ct+uCx=xECx+qCA+SCC1t+uC1x=xEC1x+qC1A+SC1K6C1Lt+uLx=xELx+qLA+SLK1LK3LCAt+uCAx=xECAx+qCAA+SCA+μCAρACAσ1HCACPt+uCPx=xECPx+qCPA+SCP+α2ρACAα2μCACN1t+uCN1x=xECN1x+qCN1A+SCN1+α1ρACAKN1CN1CN2t+uCN2x=xECN2x+qCN2A+SCN2+KN1CN1KN2CN2CN3t+uCN3x=xECN3x+qCN3A+SCN3+KN2CN2α1μCAOt+uOx=xEOx+qOA+SO+K2OSO+α3μα4ρACAK1Lα5KN1CN1α6KN2CN2E97

where C,C1,L,CA,CP,CN1,CN2,CN3,Oare the concentration of degradable substances, nondegradable substances, BOD, chlorophyll-α, dissoluble phosphorus, ammonia nitrogen, nitrite nitrogen, nitrate nitrogen, and DO, respectively; SC,SC1,SL,SCA,SCP,SCN1,SCN2,SCN3,SOare 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:

Vj=14Aj1/2+AjΔxj1+14Aj+1/2+AjΔxjE98

Figure 7.

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

where Vjis the volume of the equilibrium area at cross section J; Aj12,Aj+12are the areas of the cross sections of the inlet and outlet in the equilibrium area, respectively (Aj12=Aj1+Aj2, Aj+12=Aj+Aj+12). Then, the volume of the equilibrium area is:

Vj=18Aj1/2+3AjΔxj1+18Aj+1/2+3AjΔxjE99

Then, the mass variation of water quality variables in the equilibrium area at Δttime interval:

Δm=Vji+1Cji+1VjiCjiE100

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 Fxof water quality variables at some point in the x direction is:

Fx=uCE101

where uis the time average velocity of some point in the x direction and Cis the time average concentration of water quality variables at some point.

Then, the mass of water quality variables that passes some cross section at Δttime interval due to advection is:

m=AFxΔt=Au¯C¯Δt=Q¯C¯ΔtE102

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 Δttime interval are, respectively, as follows:

m11=Qj12Cj12ΔtE103
m12=Qj+12Cj+12ΔtE104

where Qj12=Qj1+Qj2, Qj+12=Qj+Qj+12; Cj12=θCj1+1θCj, Cj+12=θCj+1θCj+1; θis the upwind factor (0θ1). If Qj12>0, θ12;

If Qj120, θ12; if Qj12>0, θ=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:

Mm=EmCxE105

where Mmis the molecule diffusive flux at some point in the x direction and Emis the molecular diffusion coefficient and its range of values generally is 109108m2/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 Δttime interval are, respectively, as follows:

m21=Em,j12Aj12CjCj1Δxj1ΔtE106
m22=Em,j+12Aj+12Cj+1CjΔxjΔtE107

where Em,j12and Em,j+12are the molecular diffusion coefficients at xj12and xj+12, respectively; Em,j12=Em,j1+Em,j2; and Em,j+12=Em,j+Em,j+12.

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:

Mt=EtxCxE108

where Mtis the turbulent diffusion flux at some point in the x direction; Etxis the turbulent diffusion coefficient; for a turbulent flow field whose Reynolds number is about Re=104, the turbulent diffusion coefficient is approximately 3.36 × 10−4 m2/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 Δttime interval are, respectively, as follows:

m31=Etx,j12Aj12CjCj1Δxj1ΔtE109
m32=Etx,j+12Aj+12Cj+1CjΔxjΔtE110

where Etx,j12and Etx,j+12are the turbulent diffusion coefficients at xj12and xj+12, respectively; Etx,j12=Etx,j1+Etx,j2; and Etx,j+12=Etx,j+Etx,j+12.

3.3.1.4 Dispersion

The expression of dispersion is as follows:

Md=EdCxE111

where Mdis the discrete transport flux at some point in the x direction and Edis the dispersion coefficient with a range of values of 10∼103 m2/s.

Then, consider that the masses of pollutants that flow into and out of an equilibrium area by means of dispersion at Δttime interval are, respectively, as follows:

m41=Ed,j12Aj12CjCj1Δxj1ΔtE112
m42=Ed,j+12Aj+12Cj+1CjΔxjΔtE113

where Ed,j12and Ed,j+12are the dispersion coefficients at xj12and xj+12, respectively; Ed,j12=Ed,j1+Ed,j2; and Ed,j+12=Ed,j+Ed,j+12.

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:

Δm1=m11m12+m41m42=ΔtCj1θQj1/2+Ed,j1/2Aj1/2Δxj1+Cj1θQj1/2θQj+1/2Ed,j1/2Aj1/2Δxj1Ed,j+1/2Aj+1/2Δxj+Cj+11θQj+1/2+Ed,j+1/2Aj+1/2ΔxjE114

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 α0and α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 Δttime interval is:

Δm2=ΔtVjα0α1CjE115

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 Cq, and the mass of water quality variables input at a sudden point source per unit time is mm, then, consider the concentration change of water quality variables in an equilibrium area owning to lateral inflow or the sudden point source at Δttime interval is:

Δm3=ΔtCqqΔxj12+Δxj2+mmΔtE116

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 Δttime 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.

Δm=Δm1+Δm2+Δm3E117

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

VjCjVjiCji=ΔtCj1θQj1/2+Ed,j1/2Aj1/2Δxj1+Cj1θQj1/2θQj+1/2Ed,j1/2Aj1/2Δxj1Ed,j+1/2Aj+1/2Δxj+Cj+11θQj+1/2+Ed,j+1/2Aj+1/2Δxj+ΔtVjα0α1Cj+ΔtCqqΔxj12+Δxj2+mmΔtE118

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

After arrangement, we obtain:

Cj1θQj1/2+Ed,j1/2Aj1/2Δxj1+Cj1θQj1/2θQj+1/2Ed,j1/2Aj1/2Δxj1Ed,j+1/2Aj+1/2Δxjα1VjVjΔt+Cj+11θQj+1/2+Ed,j+1/2Aj+1/2Δxj=VjiΔtCji+α0VjCqqΔxj12+Δxj2mmE119

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

AjCj1+BjCj+DjCj+1=EjE120

where

Aj=θQj1/2+Ed,j1/2Aj1/2Δxj1E121
Bj=1θQj1/2θQj+1/2Ed,j1/2Aj1/2Δxj1Ed,j+1/2Aj+1/2Δxjα1VjVjΔtE122
Dj=1θQj+1/2+Ed,j+1/2Aj+1/2ΔxjE123
Ej=VjiΔtCji+α0VjCqqΔxj12+Δxj2mmE124

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=C1tE125

From the formula above, we obtain the concentration C1,0at a canal head cross section x1at any time; after transforming it into the form of the above formula (120), we obtain:

B1C1+D1C2=E1E126

where B1=1.0, D1=0, and E1=C1,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 x1per unit time is known, i.e.,

F1=Q1CE127

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:

C1θQ1+1/2Ed,1+1/2A1+1/2Δx1α1V1V1Δt+C21θQ1+1/2+Ed,1+1/2A1+1/2Δx1=Q1CV1iΔtC1i+α0V1Δx12CqqmmE128

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

B1C1+D1C2=E1E129

where B1=θQ1+1/2Ed,1+1/2A1+1/2Δx1α1V1V1Δt, D1=1θQ1+1/2+Ed,1+1/2A1+1/2Δx1, and E1=Q1CV1iΔtC1i+α0V1Δx12Cqqmm; 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=CmtE130

From the formula above, we obtain the concentration Cm,0at a canal end cross section xmat any time. After transforming it into the form of the above formula (120), we obtain:

AmCm1+BmCm=EmE131

where Am=0, Bm=1.0, and Em=Cm,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 xmper unit time is:

Fm=QmCmE132

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:

Cm1θQm1/2+Ed,m1/2Am1/2Δxm1+Cm1θQm1/2QmEd,m1/2Am1/2Δxm1α1VmVmΔt=VmiΔtCmi+α0VmCqqΔxm12mmE133

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

AmCm1+BmCm=EmE134

where Am=θQm1/2+Ed,m1/2Am1/2Δxm1, Bm=1θQm1/2QmEd,m1/2Am1/2Δxm1α1VmVmΔt, and Em=VmiΔtCmi+α0VmCqqΔxm12mm; 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 xjand the downstream cross-section is xj+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:

Cj1θQj1/2+Ed,j1/2Aj1/2Δxj1+Cj1θQj1/2QjEd,j1/2Aj1/2Δxj1α1VjVjΔt=VjiΔtCji+α0VjCqqΔxj12mmE135

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

AjCj1+BjCj=EjE136

where Aj=θQj1/2+Ed,j1/2Aj1/2Δxj1, Bj=1θQj1/2QjEd,j1/2Aj1/2Δxj1α1VjVjΔt, and Ej=VjiΔtCji+α0VjCqqΔxj12mm; 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 Cj; then by imitating the above formula (128), we obtain:

CjQj+1+Cj+1θQj+1+1/2Vj+1ΔtEd,j+1+1/2Aj+1+1/2Δxj+1α1Vj+1+Cj+21θQj+1+1/2+Ed,j+1+1/2Aj+1+1/2Δxj+1=Vj+1iΔtCj+1i+α0Vj+1Δxj+12CqqmmE137

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

Aj+1Cj+Bj+1Cj+Dj+1Cj+2=Ej+1E138

where Aj+1=Qj+1, Bj+1=θQj+1+1/2Vj+1ΔtEd,j+1+1/2Aj+1+1/2Δxj+1α1Vj+1, Dj+1=1θQj+1+1/2+Ed,j+1+1/2Aj+1+1/2Δxj+1, and Ej+1=Vj+1iΔtCj+1i+ α0Vj+1Δxj+12Cqqmm; 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:

Ed=γuA2H3E139

where γis an empirical coefficient generally having a range of values of 0.5–0.6; uis the frictional velocity; u=gHJ; Jis the hydraulic gradient; Ais the cross-section area; and His the cross section water depth.

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

PollutantReaction coefficientValueUnitRange of values
Nondegradable substancesu08day−1
u10.8gm−3 day−1
Algaeμ0.075day−11.0–3.0
ρA0.1day−10.005–0.5
σ11mday−10.153–1.83
Nitrogenɑ10.085gg−10.08–0.09
KN10.03day−10.01–0.5
KN21.5day−10.5–2.0
Phosphorusɑ20.0135gg−10.012–0.015
CBODK10.15day−10.1–0.2
K30.18day−1
OxygenOS9.08gg−1
ɑ31.6gg−11.4–1.8
ɑ42gg−11.6–2.3
ɑ53.5gg−13.0–4.0
ɑ61.07gg−11.0–1.14
K25day−10–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 m3/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 hestimated according to the open channel uniform flow formula Q=AR2/3in1= 11.20 m;

    2. the area of flow section A1069.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=gRi= 0.11 m/s;

    6. estimated channel diffusion coefficient Dx=Dy=0.15hu= 0.18 m2/s; and

    7. estimated longitudinal diffusion (including dispersion) coefficient D=6.01hu= 7.4 m2/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 m3/s within 15 min, and opened to 120 m3/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 structurePile no. (km)Length of channel pool (km)Elevation of sluice bottom (m)Design flow (m3/s)Design water level (m)
36Check gate in Anyang River717.04514.32185.6023592.67
Release sluice in Zhanghe River730.62385.19120/
Check gate in Zhanghe River731.36682.5023591.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.

ItemContentQuantity
Pollutant categoryNondegradable substance1
Mode of occurrenceAn instantaneous point source1
The channel pool where an accident occursThe 36th channel pool2
Pollutant mass1, 5, and 10t3
Place of occurrence10, 30, 50, 70, and 90% of channel pool length5
Channel pool flow30, 50, and 70% of design flow3
Simulation time/step24 h/10 min1

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 (T0), concentration peak (Cmax), and the occurrence time of concentration peak (TCmax), 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 24; 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 accidentPollutant mass (t)30% of design flow50% of design flow70% 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 L11000.3843560600.3679380400.3507260
5701.9212560401.8397380301.7535260
10603.8425560403.6793380203.507260
0.3 L1600.4107450400.3939310200.3771210
5402.0535450201.9693310101.8856210
10304.1071450203.9385310103.7711210
0.5 L1200.4965280100.475200100.4577130
5102.4823280102.3752200102.2887130
10104.9647280104.7504200104.5773130
0.7 L1100.5387220100.5146150100.4975100
5102.6935220102.5729150102.4873100
10105.3871220105.1457150104.9745100
0.9 L1100.8862110100.819980100.796350
5104.4308110104.099680103.981650
10108.8615110108.199380107.963350

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.

© 2018 The Author(s). Licensee IntechOpen. This compact is distributed under the terms of the Creative Commons Attribution-NonCommercial 4.0 License, which permits use, distribution and reproduction for non-commercial purposes, provided the original is properly cited.

How to cite and reference

Link to this compact Copy to clipboard

Cite this compact Copy to clipboard

Yu Tian, Hezhen Zheng, Xiaohui Lei and Wei Dai (November 21st 2018). Simulation Technology for Hydrodynamic and Water Quality in the Main Canal, Emergency Operation Technologies for Sudden Water Pollution Accidents in the Middle Route of South-to-North Water Diversion Project, Xiaohui Lei, IntechOpen, DOI: 10.5772/intechopen.82022. Available from:

compact statistics

167total compact downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Emergency Operation Technologies for Sudden Water Pollution Accidents in the Middle Route of South-to-North Water Diversion Project

Edited by Xiaohui Lei

Next compact

Traceability Technology for Sudden Water Pollution Accidents in Rivers

By Weihong Liao, Zhiguo Gan, Jiabiao Wang and Hezhen Zheng

Related Book

First chapter

Spatial Decision Support System (SDSS) for Stormwater Management and Water Quality Assessment

By Nally Kaunda-Bukenya, Wubishet Tadesse, Yujian Fu, Teferi Tsegaye and Mezemir Wagaw

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us