## Abstract

This study explores the effects of water quality simulation results by embedding axial dispersion into the classical advective-reactive model in a water distribution system. The Eulerian-Lagrangian operator-splitting method is employed to solve the model with axial dispersion. Satisfactory results were obtained after the numerical solution was tested against the analytical and other numerical solutions. The water quality simulation results without the reaction item show that when water velocity is low (or Peclet numbers are small), dispersion is dominant and significantly affects the numerical simulation results. The contaminant concentration at downstream node gradually increased with time along the upstream pipelines from the source, which is particularly reflected in the terminal. The simulation results show that the biomass concentration may demonstrate synthetic effects of axial dispersion and reaction, i.e. mutual promotion, given the multicomponent (substrate, residual chlorine, and microbial biomass) reaction-transport processes. It is particularly reflected in the low flow velocity.

## INTRODUCTION

Contamination is one of the main threats to drinking water security because of invasion, microbial growth, and pipeline corrosion, which lead to water quality deterioration and potential risk to human health (Machell & Boxall 2014; Davis & Janke 2015). Contaminant detection technology with water quality simulations is a technical challenge; this has been proposed to detect and manage contaminant events.

*C*is the concentration of contaminant (mg/L),

*u*is the average velocity flow of the pipe (m/s),

*x*is the axial distance in the pipe (m),

*t*is the time (s) and the left and right sides of the equation represent the transport equations and the reaction term, that is,

*r*[

*C*(

*x,t*)], respectively.

*Pe*), which are defined as

*ux*/

*D*, represent the ratio of advection and dispersion, where

_{s}*D*is the dispersion coefficient (m

_{s}^{2}/s). In the water distribution system (WDS),

*Pe*is much greater than 1, because

*D*(e.g. chlorine, bacteria and substrate) is much less than 1, thus, the advective effect is larger than the dispersive effect. However, when the time of accumulation is long, and the distance is far from the original source, the dispersion effect is very apparent, especially at low flow speed. Considering the actual situation of the WDS, water commonly flows at low speed, potentially because of the pipe network layout, daily and seasonal consumer water patterns and maintenance. Although the contribution of axial dispersion is frequently less important than that of the advection, the axial dispersion effect cannot be omitted (Axworthy & Karney 1996) from the general water quality simulation in WDS. The additional axial dispersion by the logical model in this research is proposed as follows (Axworthy & Karney 1996): Some researchers have considered dispersion with regard to water supply systems (Axworthy & Karney 1996; Zhang

_{s}*et al.*2004; Yang

*et al.*2008) in a simple pipe network and under steady state. However, the simulation of water quality with dispersion in an actual entire WDS under an unsteady state is attractive. The unsteady state (Kim & Kim 2017) is achieved when an actual hydraulic condition is under transient flow velocities and pressure associated with system operation, pressure control activities and customer behaviour.

The advective–dispersive–reactive equation can be solved as either an analytical or a numerical solution. The numerical solution is more effective and common because the analytical solution exists only for special and simplified cases. The Eulerian-Lagrangian (E-L) method is one of the most accurate and efficient solutions for the equation (Wang & Lacroix 1997; Basha & Malaeb 2007; Hu *et al.* 2014). The widely favourable approach in the literature is to use split operator technology that separates the advective, dispersive, and reactive components of the problem, such as the Eulerian-Lagrangian operator-splitting method (ELOS) (Xin & Wong 1989; Scardovelli & Zaleski 2003).

This study develops a more robust mathematical model that reflects water quality. The model is an embedded axial dispersion combined with a hydraulic network under unsteady-state flow conditions. To solve the advection–dispersion–reaction equations, the ELOS method is used to divide the grids of the time and space domains in a numerical solution.

## METHODOLOGY

A water supply network, as a node-link system, is composed of links (pipes, pumps, and valves), nodes (junctions and tanks), and sources. According to the network, a hydraulic model is applied to simulate the dynamic flows in pipelines during a set of hydraulic time steps. The flow dynamics at a pipe junction control solute mixing and concentrations in downstream pipes (Yu *et al.* 2016). The effect can lead to different outcomes of water quality modelling in a distribution network. The present study hypothesizes that all the water with contaminants from differently-sized pipes is mixed instantaneously and completely at the node, as well as in the tank and storage, because of performable and convenient comparison.

A water quality model determines how the contaminant concentration varies with time throughout the network under specified hydraulic conditions and contamination input patterns. As shown in Appendix II (available with the online version of this paper), the water quality model hypothesized the one-dimensional equation combined with dispersion as Equation (2) during a single or consecutive unsteady-state network flow hydraulics and ideal plug flow. The water quality model combines a description of the following processes: bacteria growth, mortality, attachment and detachment, inactivation by chlorine, substrate consumption, chlorine decay and radial mass transport.

### Numerical solution for water quality model

Before solving the numerical solution for water quality model, a schematic is utilized to demonstrate fluid-flow processes as shown in Figure 1, and the summary of fluid-flow processes is as follows.

- 1.
- 2.
Fixed profiles (

*P*_{1}*, P*_{2}…*P*_{n}_{−1}) and removable profiles (*S*_{P}_{1}*, S*_{P}_{2}…*S*_{Pn}_{−1}) are arranged based on space steps. Fixed profiles are cut in a fixed pipeline, and removable profiles are used to divide the bulk flow. - 3.
The initialized concentration of contaminants are distributed by profiles, as follows: the concentration distribution in bulk flow of the pipeline is

*C*(_{bj}*C*_{b}_{1}*, C*_{b}_{2}…*C*), and the concentration distribution in the wall zone of the pipeline is_{bn}*C*(_{wj}*C*_{w}_{1}*, C*_{w}_{2}…*C*)._{wn} - 4.
In Δ

*t*time, the amount of Δ*x*water enters the pipe (*C*), the water flows forward and the sharp profiles move to a position, as shown in Figure 1, state b. Synchronously, the amount of Δ_{0}*x*water outflows from the pipe, and its concentration is*C*._{bn} - 5.
Hereto, all the events are completely updated (dispersion, advection and reaction) within Δ

*t*time. The current state is the starting state of the next moment, followed by repetition of the above steps.

*t*(

*t*

_{n},

*t*

_{n}_{+1}). Hence, the guideline of the ELOS method is similar but different from that of Zhao

*et al.*(2011), the formula is as follows: where

*C**,

*C***,

*C****, and

*C***** are the transient concentrations, and

*n*is the time step.

#### Advective step

*C°*is at the Lagrangian mesh point, and

*u*× Δ

*t*is the distance between the Eulerian and Lagrangian meshes. The concentration of the characteristic curve is the same at every point of the dotted line, that is,

*C**(

*t*

_{n}_{+1/2},

*x*) =

_{j}*C°*. The value of

*C*° is obtained via the interpolation method, and then

*C**(

*t*

_{n}_{+1/2},

*x*) is obtained.

_{j}#### Dispersive step

*C***(

*t*) is obtained by

_{n},x_{j}*C**(

*t*

_{n}_{+1/2}

*,x*). However,

_{j}*C***(

*t*

_{n}_{+1}

*,x*),

_{j}*C***(

*t*

_{n}_{+1}

*,x*

_{j}_{+1}), and

*C***(

*t*

_{n}_{+1}

*,x*

_{j}_{−1}) are all unknown. The dispersion step is solved using the Thomas algorithm, which avoids the stability restrictions associated with explicit methods.

#### Mass transfer and reaction in radial direction

## APPLICATION

### Test a single pipeline

A single pipeline (*L* = 1,000 m) is used to test the water quality simulation with ELOS method in (1) water velocity *u* = 0.1 m/s and (2) water velocity *u* = 0.01 m/s. The concentration of the source contaminant is as follows: Type a *C*(0,*t*) = 1.0 mg/L and Type b *C*(0,*t*) = exp(–0.001 × *t*). Other parameters include Δ*x* = 5 m and Δ*t* = 100 s, and the dispersion coefficient is *D _{s}* = 0.1 m

^{2}/s.

To confirm the feasibility of ELOS method, the activated stabilized oxygen (ASO) (Zhang *et al.* 2004), the E-L method (Basha & Malaeb 2007), the finite difference method (FDM) used by Crank-Nicholson schemes and the analytical solution derived in Appendix I (available with the online version of this paper) are employed to simulate the water quality model, which is simplified to omit *r*[*C*(*x*,*t*)], and compared in a single pipeline. At the same time, EPAnet-MSX (Shang *et al.* 2011; Seyoum *et al.* 2013), as a numerical solution, is used to deal with the advective–reactive equation. When *r*[*C*(*x*,*t*)] = 0, EPAnet-MSX degrades into EPAnet.

- (i)
- (ii)
Waviness or total variation of error

*W*represents the oscillation degree of the error of the numerical solution compared with analytical solutions,

*C*

_{analytical}_{,i}is the analytical solution of the pipe section

*i*,

*C*

_{numerical}_{,i}is the numerical solution of the pipe section

*i*, and

*C*

_{analytical}_{,i}and

*C*are the normalization.

_{numerical,i}### Test a network

A test network is located at Guangdong Province, China, where the township area is approximately 130 km^{2}, and the population is roughly 100,000 persons. The schematic of the test network represents a network consisting of 102 pipes, 85 nodes, five pump elements, one tank, and two supply sources. The network hydraulic parameters (i.e. pipes, nodes, tanks, pumps and demand) are from an actual investigation (the data file in https://www.researchgate.net/profile/Guoqiang_Chen4**)**. Water patterns stand for the time pattern used to vary external source strength over time (Table 1). The network system operates with actual hydraulic conditions and hypothetically assumes a contaminant-intrusion source. The numerical simulation results of this model are compared with those of the models without axial dispersion.

Time . | Pattern . | |||||
---|---|---|---|---|---|---|

1–12 h | 1 | 1.2 | 1.4 | 1.6 | 1.4 | 1.2 |

12–24 h | 1 | 0.8 | 0.6 | 0.4 | 0.6 | 0.8 |

Time . | Pattern . | |||||
---|---|---|---|---|---|---|

1–12 h | 1 | 1.2 | 1.4 | 1.6 | 1.4 | 1.2 |

12–24 h | 1 | 0.8 | 0.6 | 0.4 | 0.6 | 0.8 |

The comparative results of the contaminant concentrations, which are simulated for representative nodes 3, 9, 25, and 35 are shown in Figure 3 for time simulation. The information from Figure 3 is as follows. Node 3 is a typical node affected by the water quality of source 1; node 25 represents the terminal equipment far from water source 1; nodes 9 and 35 denote two concentrated water areas, zones 1 and 2, respectively.

The following two cases: *i*, *r*[*C*(*x,t*)] = 0; *ii*, *r*[*C*(*x,t*)] are applied to study the multicomponent (substrate, biomass, and disinfectant) reaction-transport processes (Munavalli & Mohan Kumar 2004). This water quality model is shown in Appendix II. The dispersive coefficients of the model are set up by Lu *et al.* (1995) as follows: dispersive coefficient of substrate *D _{sb}*, 3.50 × 10

^{−9}m

^{2}/s; dispersive coefficient of biomass,

*D*3.33 × 10

_{b}^{−9}m

^{2}/s; dispersive coefficient of chlorine

*D*, 3.47 × 10

_{cl}^{−9}m

^{2}/s. The other parameters are set up by Munavalli & Mohan Kumar (2004). For comparison, the EPAnet-MSX simulates a water quality model without dispersion. Its parameters, such as hydraulic time step, water quality time step, solver, relative concentration tolerance and absolute concentration tolerance, are set as 2 h, 100 s, RK5, 0.01, and 0.001, respectively.

## RESULTS AND DISCUSSION

### Comparison with other solutions in a single pipeline

Figure 4 shows the analytical and numerical solutions with different marked points. The analytical solution is exhibited by a black solid line; the numerical solution without dispersion (EPAnet) is exhibited by a purple dotted line. The solid line constitutes a sharp profile, which is different from the dotted line in the middle of the pipeline where contaminants caused by the dispersion effect reached (Figure 4). When *u* is larger, the sharp profiles have increased gradient until the sharp profiles overlap with the dotted line. Three inferences are as follows:

- 1.
The gradient of the sharp profile is relevant to

*P*rather than_{e}*t*; if*P*is larger, then the gradient is larger, but the dispersion effect is smaller, and vice-versa._{e} - 2.
As a constant

*P*, the constant concentration profile is integrally moving along with_{e}*t*in the pipeline. - 3.
The two zones enclosed with the solid line and dotted lines are the same in size. According to the theory of mass conservation, the contaminant that spreads forward is transferred from the backward contaminant.

The numerical solutions with dispersion revolve around the analytical solution of the up and down shocks. Compared with the analytical solution, the ELOS method and other numerical solutions are competent to solve the advective–dispersive equation. The performance is judged by the standard error and waviness (Table 2). FDM lacks stability and performs the worst, whereas ELOS and E-L methods perform the best. However, E-L method is limited when the reaction term is complex. Therefore, the ELOS method is applied in the water quality model.

Type . | Velocity (m/s) . | ELOS . | E-L . | ASO . | FDM . | ||||
---|---|---|---|---|---|---|---|---|---|

S . | W . | S . | W . | S . | W . | S . | W . | ||

a | 0.1 | 3.0 × 10^{−4} | 4.4 × 10^{−4} | 3.1 × 10^{−4} | 8.2 × 10^{−4} | 3.5 × 10^{−3} | 1.9 × 10^{−2} | 3.7 × 10^{−3} | 2.0 × 10^{−2} |

0.01 | 2.3 × 10^{−4} | 2.3 × 10^{−3} | 4.1 × 10^{−4} | 4.0 × 10^{−3} | 4.8 × 10^{−4} | 1.0 × 10^{−2} | 4.9 × 10^{−4} | 1.2 × 10^{−2} | |

b | 0.1 | 2.9 × 10^{−4} | 4.3 × 10^{−3} | 2.9 × 10^{−4} | 8.3 × 10^{−4} | 3.4 × 10^{−3} | 1.9 × 10^{−2} | 3.8 × 10^{−3} | 2.0 × 10^{−2} |

0.01 | 2.4 × 10^{−4} | 2.2 × 10^{−3} | 4.0 × 10^{−4} | 4.0 × 10^{−3} | 4.9 × 10^{−4} | 1.1 × 10^{−2} | 5.0 × 10^{−4} | 1.2 × 10^{−2} |

Type . | Velocity (m/s) . | ELOS . | E-L . | ASO . | FDM . | ||||
---|---|---|---|---|---|---|---|---|---|

S . | W . | S . | W . | S . | W . | S . | W . | ||

a | 0.1 | 3.0 × 10^{−4} | 4.4 × 10^{−4} | 3.1 × 10^{−4} | 8.2 × 10^{−4} | 3.5 × 10^{−3} | 1.9 × 10^{−2} | 3.7 × 10^{−3} | 2.0 × 10^{−2} |

0.01 | 2.3 × 10^{−4} | 2.3 × 10^{−3} | 4.1 × 10^{−4} | 4.0 × 10^{−3} | 4.8 × 10^{−4} | 1.0 × 10^{−2} | 4.9 × 10^{−4} | 1.2 × 10^{−2} | |

b | 0.1 | 2.9 × 10^{−4} | 4.3 × 10^{−3} | 2.9 × 10^{−4} | 8.3 × 10^{−4} | 3.4 × 10^{−3} | 1.9 × 10^{−2} | 3.8 × 10^{−3} | 2.0 × 10^{−2} |

0.01 | 2.4 × 10^{−4} | 2.2 × 10^{−3} | 4.0 × 10^{−4} | 4.0 × 10^{−3} | 4.9 × 10^{−4} | 1.1 × 10^{−2} | 5.0 × 10^{−4} | 1.2 × 10^{−2} |

### Test a network

#### Case i

In case *i*, *r*[*C*(*x,t*)] = 0 is considered, it represents no reaction when a hypothetical contaminant invades the network. The contaminant is transmitted in the network by advection and dispersion.

Figure 5 shows that the contaminant concentration is changed with time at nodes 3, 9, 25, and 35 at 48 h. Pollutants spread throughout the network with the flow. Pollutant concentration changes with time as well as up and down shocks; the maximum and the minimum concentrations are 10 and 0 mg/L, respectively. The simulation results of EPAnet (without dispersion) and ELOS method (with dispersion) are basically consistent at nodes 3 and 35 but inconsistent at nodes 9 and 25. Node 25 is far from source 1; the water with contaminant without dispersion begins from source 1 and reaches node 25 after approximately 44 h, as shown by the continuous line in Figure 5 (node 25). However, after approximately 40 h, the pollutant concentration with dispersion begins to change, and its concentration is lower than that of the EPAnet-MSX. In the analysis, the dispersion effect occurs before the advection effect; the upstream mass is dispersed to the downstream based on material balance.

The analysis of flow direction at nodes 3, 9, 25 and 35 shows that the upstream pipelines link 52, 35, 16, and 3, whose flow direction is constant (Figure 6). The contaminant concentration of these nodes is mainly affected by the upstream pipelines rather than the other linked pipelines. The flow velocity of the upstream pipelines changes with time (Figure 7). The velocity of link 3 linked to node 3 is the largest (>0.55 m/s); the velocity of link 16 linked to node 25 is the smallest (0.01 m/s). The velocity of link 52 linked to node 9 ranges from 0.1 to 0.35 m/s; the velocity of link 70 linked to node 35 ranges from 0.2 to 0.6 m/s. The other evident phenomenon is that the velocities of all links are homologous to the water pattern. When the flow is in low velocity (during 0–5 and 24–30 h), the dispersion effect may be dominant. These results and the analysis in Figure 5 show that when the velocity is larger than 0.2 m/s, the dispersion effect is neglected. When the velocity is smaller than 0.1 m/s, the dispersion effect becomes prominent. The dispersion effect is not prominent during 0–5 h at node 9 because the contaminant is absent, but it becomes prominent during 24–30 h. The contaminant disperses forward along with the flow. The contaminant concentration at downstream node is accumulates with time along the upstream pipelines from the source. Thus, the contaminant concentration of the downstream nodes is gradually affected. Terminal node 25, which is far from the source, shows the largest difference between ELOS method and EPAnet-MSX.

#### Case ii

Case *ii* is applied to study the bacterial growth model (Munavalli & Mohan Kumar 2004), which was developed by multicomponent (substrate, biomass, and disinfectant) reaction-transport processes. The model utilizes the simplified expressions for the basic processes (in bulk flow and at pipe wall), e.g. bacterial growth and mortality, attachment to and detachment from the surface, substrate utilization and disinfectant action. Regular initial conditions are assumed, and the AOC, free-chlorine and free-biomass concentration in the water from source 1 are set to 1.065, 1.5 and 0.0006 mg/L, respectively. The conditions at source 2 are set to 0.168, 0.49, and 0.0001 mg/L.

Figure 8 shows that the concentrations of AOC and Chlorine using EPAnet-MSX and ELOS method are basically consistent, which agrees with the results in case *i*. However, the biomass concentration features a large difference. The simulation results of the ELOS method are larger than that of the EPAnet-MSX. Biomass is sensitive to AOC and chlorine concentrations. As shown in Figure 8(a), the concentrations of AOC and chlorine change weakly. The changes are slightly evident (the difference of AOC < 0.3 mg/L, chlorine <0.1 mg/L) only when the time is 35–45 h at node 9. However, the biomass concentration evidently changed (the concentration of ELOS is approximately twice that of EPAnet-MSX) because of the synthetic effects of axial dispersion, radial mass transfer and reaction. Figure 8(b) shows that node 25, which is the terminal point away from the sources, has similar characteristics to that in Figure 8(a). Node 3 is extremely close to source 1. Thus, no sufficient accumulation of the synthetic effects is achieved. The substrate and residual chlorine are the overall attenuation trend in the pipeline. Microbial metabolism consumes organic matter (Şeker *et al.* 1995). Residual chlorine decay (Ozdemir & Ger 1998) contains two aspects, as follows. Residual chlorine, which is a disinfectant, combines with microbial enzyme for microbial inactivation and generates the disinfection by-products of REDOX reaction, which also shows an attenuation trend. The attenuation is a function of the retention time. The microbial quantity changes with the concentration of substrate and chlorine. The dispersion effect on the multicomponent accelerates the radial mass transfer and reaction variation in unit time. When the flow velocity is low, the relative retention time is long; the overall variation of multicomponent is significant. In summary, it synthetically affects the axial dispersion and radial mass transfer. The reaction effect is larger than that of the model without dispersion.

## CONCLUSIONS

In the classical water quality model, axial dispersion is generally ignored during simulation in WDS. In this study, the classical model is embedded in axial dispersion and applied to numerically simulate the hypothetical contaminant-intrusion from the water source under hydraulic unsteady conditions. The ELOS model is used to solve the advection–dispersion–reaction equations and compared with the analytical solutions and other numerical solutions.

In single pipeline, the sharp profile caused by the dispersion effect can be observed. Its gradient is relevant to *P _{e}*, integrally moving along with

*t*in the pipeline. The contaminant that spreads forward is transferred from the backward contaminant. In an entire WDS, when

*u*> 0.2 m/s, the dispersion effect can almost be neglected. However, when

*u*< 0.1 m/s, the dispersion effect begins to become prominent. The contaminant concentration at downstream node is gradually accumulated with time along the upstream pipelines from the source, which is particularly reflected in the terminal. The simulation results compared with EPAnet-MSX show that biomass concentration may exhibit synthetic effects of axial dispersion and radial mass transfer and reaction given the multicomponent reaction transport processes. The dispersion effect features a large difference between the bulk flow and wall zone to mutually increase the radial mass transfer and reaction. Therefore, biomass concentration significantly changes, and this phenomenon is particularly reflected in the low flow velocity.

With further simulation of water quality in WDS, more realistic conditions with incomplete mixing type at the nodes and corrosion in cast iron pipes must be considered. The incomplete mixing type can provide more accurate simulations and may lead to more uneven concentration distributions, which is beneficial for dispersion. As complex biochemical reaction, the corrosion in pipes can change the multicomponent concentrations and is common in the actual WDS.

## ACKNOWLEDGEMENTS

This work is supported by the project in the National Science and Technology Pillar Program during the Twelfth Five-year Plan Period (2012BAJ25B06-003), the project of Chongqing Social Science Planning (2016BS081), and the Science and Technology Research Program of Chongqing Municipal Education Commission (KJ1706175). The valuable comments and suggestions from the editor and two anonymous reviewers are very much appreciated.

## REFERENCES

*.*