## Abstract

A two-dimensional water quality numerical model has been developed to analyse the influence of adsorption boundary on the variation of water pollutant degradation. The characteristic finite element method is adopted to solve the discrete pollutant concentration equations, so the time-monotonicity of concentration is ensured in the convection-dominated flow field with various adsorption conditions. The second-order Runge–Kutta method is adopted to track the upstream concentration source points along the streamline, so the precision of the calculated pollutant concentration is improved even under the streamline sharp-bending condition. The flow roughness and pollutant degradation coefficient are formulated by elements according to the actual situation, so the model can quantitatively analyse the effect of various ecological restoration schemes on purification capacity in water areas. The adaptability and robustness of the model were verified by a typical engineering case, where the calculated results coincide with the measured data.

## INTRODUCTION

At present, the Chinese government has increased support for water pollution control and water ecological restoration to promote the construction of ecological civilization, with many water environment renovation projects having been implemented. The main purpose of water ecological restoration projects are either to promote the regional water environmental capacity, or to improve the self-purification capacity of a water area (Jiang *et al.* 2014), if only considering the landscape effect, without the quantitative analysis of water quality and the later intensive management of pollutant discharge, the project is likely to be unable to achieve the desired goal (Qin *et al.* 2002). Therefore, the quantitative analysis of water quality is very important for the determination of the layout of the water ecological restoration project and the corresponding water body operation scheme (Xu *et al.* 2017).

In a water ecological restoration project, in addition to taking measures to increase the water area and improve the replacement rate of water (Wang *et al.* 2017), various adsorption boundaries are generally used to purify the water quality, which has an important influence on the later maintenance of the ecological restoration project. For example, in order to improve the existing water purification capacity, an artificial wetland is constructed in the surrounding lake (Hao *et al.* 2017; Huang *et al.* 2017), by circulating water adsorption, to achieve the partial purification of lake water; for the lifting of the water landscape and improving water quality, the landscape planning department usually adopts floating plants and aquatic plants in the waters during the growing season (Guan *et al.* 2014; Lyu *et al.* 2017; Zhang *et al.* 2017), which can adsorb contaminants in water and enhance its purification capacity, however, the plants inevitably decompose due to death and decay (Tang *et al.* 2013; Qi *et al.* 2017).

With the development of science and technology some new measures appeared, such as artificial aeration, adsorption film, nano-zeolite bed, physical and chemical reagents, etc. (Chen & Lu 2007; Liu *et al.* 2011; Qiang *et al.* 2014; Wang *et al.* 2016; Wu *et al.* 2017). The common feature of these measures is that the pollutants are partially adsorbed or degraded when they pass through the above interface, thus, improving the self-purification capacity of the water. In addition, some adsorption boundary in the change of water self-purification ability will also change the water flow conditions at the same time, such as aquatic plants making the local roughness increase and the turbulent diffusion coefficient decrease, and it is easy to form a retention area or a backwater region (Li *et al.* 2017; Zheng *et al.* 2017) which has a negative impact on the overall efficiency of water replacement. To solve the above problems, a two-dimensional flow and water quality numerical model has been developed to quantitatively analyse the influence of an adsorption boundary on the flow field and pollutant degradation law in the water area, which establishes a foundation for future engineering application.

## BASIC THEORY OF THE MATHEMATICAL MODEL

### Governing equations

*u*and

*v*are flow velocity, (m.s

^{−1});

*h*and are the water depth and bed elevation, (m);

*E*is the turbulent exchange coefficient, (m

^{2}.s

^{−1});

*C*is the Chezy coefficient, ,

*n*is roughness;

*c*is the concentration of pollutants in the water, (mg.L

^{−1}); is the diffusion coefficient for pollutants (m

^{2}.s

^{−1});

*q*is the source/sink discharge (s

^{−1}); is the concentration of pollutants in the source/sink discharge (mg.L

^{−1}). The pollutant degradation term is often described by the first-order kinetic model: , where

*k*is the degradation coefficient of pollutant in water, (d

^{−1}).

In the pollutant concentration equation, the effects of the adsorption boundary on the distribution of pollutant concentration are mainly described in the pollutant degradation term, where the degradation coefficient *k* value is dependent on the water temperature, velocity, adsorption boundary characteristics and the initial concentration of pollutants (Guo *et al.* 2008). For this reason, the finite element method is used to solve the discrete equations, so that the parameters of the flow and concentration can be specified separately for each numerical element, and by modifying the code, we can also use a specific pollutant degradation function.

Under the action of the flow field and the adsorption boundary, the pollutant concentration should be attenuated along the streamline and, for any point in the water, the pollutant concentration increases with time and then tends to the equilibrium value. Thus, a characteristic finite element method is adopted to guarantee the monotonicity of concentration under the strong convection flow condition and the adsorption boundary.

The flow boundary conditions are usually the flow velocity of the nodes at the inflow boundary and the water level of the nodes at the outflow boundary. The pollutant concentration boundary conditions include two situations: one is to set the pollutant concentration of the nodes on the inflow boundary, which can simulate the transport and degradation of upstream pollutants along the river flow. The other is to set the concentration of the nodes on the sidewall boundary, which can simulate the discharge process of pollutant sources along the river. The water flow can be either constant or unsteady, but in any case, the movement of pollutant is always an unsteady process.

### Characteristic finite element method

For the convection-dominated transport problem, the variation of pollutant concentration at any point in the flow field is mainly relative to upstream pollutant transportation besides local concentration diffusion. However, the traditional finite element method is equivalent to the central difference, which needs to increase the additional numerical viscosity to suppress the oscillations that easily cause local distortion. Therefore, a characteristic finite element method is adopted in this paper. It only needs to source the characteristic point along the streamline and to simplify the transport equation. The specific methods are as follows.

*K*at the time, and the position coordinates of the point can be expressed as: The transport equation of pollutant concentration is discrete as follows: By substituting Equation (5) into the above equations, the finite discrete element equation in each element is obtained:

In the equations above, the shape functions in the element are expressed as ; the concentrations of the three nodes are ; the characteristic concentrations of the nodes are ; and the source sink node concentrations are . Through numerical integration, the element coefficient matrix and the right end vector are obtained for each element and assembled into the global matrix to calculate the node concentration in the whole domain.

### The determination of characteristic points

*K*should be obtained through the integrated solution in Equation (6). The specific method is shown in Figure 1, starting from the node position, tracked upstream along the streamline and finding the characteristic point. Considering the error caused by streamline bending, the second-order Runge–Kutta method is used: where and are the position coordinates of the concentration characteristic point

*K*; and are the position coordinates of the node

*P*; and are flow velocities for the node; and are velocities of the halfway characteristic point

*M*. Using the upstream tracking along the streamline to find the characteristic points, the accuracy of the characteristic concentration is improved by one order.

After obtaining the coordinates of the characteristic points (also including its halfway characteristic points), we need to find the element in which the points are located, better known as the characteristic element. And so, the characteristic concentration value can be interpolated according to the shape function and the three node concentrations at the time . In the program design, the numbers of characteristic elements for all nodes are stored in the first step. At the next step, the program initially starts the characteristic element check from the pre-stored element. If the characteristic point is not in this element, the program performs a bi-directional search from the stored number of characteristic elements to accelerate the computing speed. For the constant problem, the search for the characteristic element number is only executed once.

## MODEL VERIFICATION AND ANALYSIS

The Luan-he river regulation project of Chengde city in the Hebei Province of China was taken as an example. The layout of the regulation project is shown in Figure 2. The main stream is about 4.7 km long, and averages 300 m wide. The tributary is about 1.6 km long and averages 200 m wide. The research scheme includes two cases. The first case is to analyse the variation of the flow field and pollutant distribution after the traditional river course has been carried out. The second case analyses the effects of the ecological restoration projects in specific regions on the purification capacity of the water area. In the calculation, using the river organic pollutant chemical oxygen demand (COD) as research index, according to the annual environment bulletin of water quality in Chengde city in 2011–2012, the water quality at the monitoring sections of the upstream Shang-ban-cheng bridge at inflow boundary I and the downstream Wu-long-ji bridge at the outflow boundary are all Grade IV–V, so the concentration of COD at inflow boundary I of the main stream is 40 mg.L^{−1}, which corresponds to Grade V water quality. The water quality of the tributary is clearly better than the main stream because of the completed ecological restoration, so the pollutant concentration at inflow boundary II of the tributary is 2.0 mg.L^{−1}. The flow of the main stream at inflow boundary I is 500 m^{3}.s^{−1}, the tributary flow at inflow boundary II is 100 m^{3}.s^{−1}, and the water level at the downstream outflow boundary is controlled at 275 m, so that both flows from the main channel and tributary will be mixed at the junction area and diluted and attenuated at the lower reaches of the river channel. The computational domain is divided into four regions where the different parameters are provided for water flow and water quality calculations according to the design schemes. A series of formulas for the turbulent exchange coefficient were evaluated by Najafzadeh & Tafarojnoruz (2016), wherein the Kashefipour & Falconer (2002) formula is adopted to calculate the turbulent exchange coefficient, and the diffusion coefficient of pollutant is 0.9 times the turbulent exchange coefficient. River roughness and pollutant degradation coefficient are selected in different parts of the river (Guo *et al.* 2008; Zhang 2012), with the detailed parameter settings shown in Table 1. The calculation conditions above remain constant.

Case . | River functional area . | Turbulent exchange coefficient (m^{2}.s^{−}^{1})
. | Roughness . | Pollutant degradation coefficient (d^{−1})
. | Case description . |
---|---|---|---|---|---|

1 | (1)Beach land | 0.9 | 0.035 | 0.18 | No ecological restoration measures, the degradation coefficient of pollutant remains unchanged |

(2)River band | 0.6 | 0.04 | 0.18 | ||

(3)Main channel | 1.2 | 0.022 | 0.18 | ||

(4)Beach land | 0.9 | 0.035 | 0.18 | ||

2 | (1)Beach land | 0.9 | 0.035 | 0.18 | The degradation coefficient of pollutant was increased due to the planting of aquatic macrophytes in region (4) |

(2)River band | 0.6 | 0.04 | 0.18 | ||

(3)Main channel | 1.2 | 0.022 | 0.18 | ||

(4)Ecological restoration | 0.6 | 0.08 | 0.35 |

Case . | River functional area . | Turbulent exchange coefficient (m^{2}.s^{−}^{1})
. | Roughness . | Pollutant degradation coefficient (d^{−1})
. | Case description . |
---|---|---|---|---|---|

1 | (1)Beach land | 0.9 | 0.035 | 0.18 | No ecological restoration measures, the degradation coefficient of pollutant remains unchanged |

(2)River band | 0.6 | 0.04 | 0.18 | ||

(3)Main channel | 1.2 | 0.022 | 0.18 | ||

(4)Beach land | 0.9 | 0.035 | 0.18 | ||

2 | (1)Beach land | 0.9 | 0.035 | 0.18 | The degradation coefficient of pollutant was increased due to the planting of aquatic macrophytes in region (4) |

(2)River band | 0.6 | 0.04 | 0.18 | ||

(3)Main channel | 1.2 | 0.022 | 0.18 | ||

(4)Ecological restoration | 0.6 | 0.08 | 0.35 |

Figures 3 and 4 show the river flow field and the distribution of pollutant concentration after the conventional river regulation scheme. In the downstream of the junction area a high-velocity zone appears with a velocity magnitude of approximately 0.14 m.s^{−1} and is mainly located in the middle of the river. Due to the clear water discharging in the tributary, there appear obvious concentration demarcations in the junction area. The main pollutants in the river are on the right side, and along the left bank the water quality is better. This is until the downstream outlet section where the peak value of pollutant concentration is up to 30–36 mg.L^{−1}, which still belongs to the Grade IV–V water.

Figures 5 and 6 are changes of river flow field and pollutant concentration distribution after the implementation of ecological restoration in a specific region. In the ecological restoration area on the right bank, the aquatic plants are densely planted to increase the local roughness, which causes the high-velocity area on the left bank side, that subsequently increases the pollutant mixture of mainstream and tributary. At the same time, the degradation coefficient increases and the pollutant concentration decreases in the ecological restoration area. The peak value of pollutant concentration decreases to 25–30 mg.L^{−1} at the downstream outlet, which belongs to Grade III–IV water.

Following comprehensive comparison, ecological restoration in the second case was performed by design in 2012. According to the environmental bulletin of Chengde city in 2015–2016, the water quality monitoring section of the downstream Wu-long-ji bridge was of Grades IV and III respectively, wherein the COD decreased to 20–30 mg.L^{−1}. This tendency shows consistency with the conclusion of the model prediction.

## CONCLUSIONS

Adsorption boundaries exist widely in nature and bring great differences in river pollutant degradation. Now, in many river regulation projects of China, some ecological or artificial restoration measures with adsorption or absorption characteristics have been used to enhance water purification capacity. The quantitative study of these adsorption boundary effects is very important for ecological environment preservation in the future. Therefore, a two-dimensional flow water quality mathematical model is developed to analyse the influence of an adsorption boundary on the concentration degradation of water pollutants. By using the characteristic finite element method, the model can ensure the time-monotonicity of the calculated concentration, even in the convection-dominated flow field with various adsorption conditions. The second-order Runge–Kutta method is used to track the location of the characteristic points along the streamline, and the precision of the calculated concentration is improved under the streamline bending conditions. By specifying the flow roughness and pollutant degradation coefficient in each element, the actual action of various sub-regional water quality improvement measures can be considered. The ecological restoration case, an actual river regulation project, is analysed and the calculated results show that the COD concentration at the downstream outlet decreases from the original 30–36 mg.L^{−1} to 20–30 mg.L^{−1} after implementing ecological restoration, which shows consistency with the measured data. In practice, the model can be used to quantitatively evaluate various ecological restoration measures in the water area, and provide support during the decision-making process.

## ACKNOWLEDGEMENTS

The work described in this paper was supported by a grant from the National Key Research and Development Program of China (No. 2016YFC0401708) and the National Natural Science Foundations of China (No. 51709278 and No. 51679261).