Impact of diffusion and dispersion of contaminants in water distribution networks modelling and monitoring

Water distribution systems connect consumers to water resources by using hydraulic components, such as pipes, valves and tanks, that are usually vulnerable elements; these elements may represent a pathway for contamination intrusion (Blokker et al. 2018).

Adequate water quality in distribution networks is a fundamental requirement that must be guaranteed to safeguard public health. The aim of maintaining water quality standards is important for the entire integrated water system (Freni & Sambito 2017). Water quality monitoring is an indispensable pre-requisite, and it can be achieved using a variety of methods. Several studies have been carried out in the literature on the optimal positioning of sensors in water distribution networks. Such studies utilise optimization methodologies often derived from other scientific sectors such as mathematics, statics or informatics (Salerno & Rabbeni 2018). Boccelli et al. (1998) formulated an optimisation model for the dynamic planning of disinfectant injections in order to minimise the total dose required to meet the water supply safety constraints. The proposed approach automatically generates the linear programming formulation of the optimal programming model, which is then solved using the simplex algorithm. The results of these applications suggest that the disinfection recall effect can reduce the amount of disinfectant needed to meet concentration constraints, providing an alternative to conventional disinfection applied only at the origin. In Tolson et al. (2004), an approach was studied that connects a genetic algorithm (GA), as an optimisation tool, with the first order reliability method (FORM), used for the estimation of the reliability of the network capacity. The latter is defined as the likelihood of meeting minimum permissible pressure constraints throughout the network under uncertain nodal demands and uncertain conditions of tube roughness. The combined use of a GA and the FORM led to an effective approach for the optimisation of reliability based on water distribution networks. Ozdemir & Ucaner (2005) studied the optimal positioning of a chlorination treatment station through the use of genetic algorithms and EPANET software (Rossman 2000). The objective of the study was to meet the minimum and maximum demand for residual chlorine at every point of the network, minimising the consumption of chlorine as much as possible (Ozdemir & Ucaner 2005). In Gibbs et al. (2006), a linear regression model, the Multi Perceptron layer (MPL), and two artificial neural networks used with general regression neural network (GRNN) were utilised to determine the optimal chlorine dose to avoid applying chlorine concentrations that were too low, which could result in bacterial regrowth, or chlorine concentrations that were too high, which could result in the formation of by-products. In Preis & Ostfeld (2008), the problem of optimal sensor positioning is addressed in order to improve the safety of a water distribution system against the intentional intrusion of a contaminant. In that study, a multi-objective model was developed for the optimal positioning of a sensor in a water distribution system using a non-dominant GA II and was validated using two increasingly complex water distribution systems.

In the literature, a model was developed for the analysis of the flow model called Comparison of Flow Pattern Distributions (CFPD). This method allows the user to compare flow patterns of a supply area and distinguish consistent from inconsistent changes in them. A CFPD analysis is made of all possible combinations of time blocks of a preselected length of the comparison frame within the complete dataset (van Thienen et al. 2013).

All the above-mentioned studies rely on models to present a robust simplified estimation of contaminant distribution and concentrations. The use of simplified models may make an unreliable estimation of contaminant propagation, thus taking to the deployment of inefficient monitoring networks (Piazza et al. 2017).

The great majority of literature research is based on hydraulic simulation tools, such as EPANET, which adopt a simplified approach to water quality based on advective transport and some simplified reaction kinetics.

Even if such simplifications are adequate for many practical applications, dispersive/diffusive transport processes become relevant when flow velocity is low and Reynolds numbers are under 50,000, such as frequently seen in urban water distribution networks during the night (Axworthy & Karney 1996). Moreover, optimisation studies were often based on numerical and modelling analyses without any comparison with experimental data.

Many authors have dealt with the problem of advection-dispersion transport. Axworthy & Karney (1996) analysed the solute transport problem using both the advective model and the advection-dispersion equation. Combining the limit velocity and the coefficient dispersion, a relationship is produced that regulates the dispersion sensitivity to flow velocity in a particular position. When velocities are high, the dispersion is probably negligible, and an advanced transport model should be suitable for any analytical need. A general rule for neglecting dispersion phenomena was not found.

Dong & Selvadurai (2008) developed a model capable of modelling the advection-dominated transport process; this process was based on an operator-splitting Euler-integration-based Taylor-Galerkin scheme, which models a problem that exhibits spherical symmetry. Williams & Tomasko (2008) presented the analytical solution of a one-dimensional contaminant transport undergoing advection, dispersion, sorption, and first-order decay, which was subject to a first-order decaying contaminant concentration at the source and a Type I Dirichlet boundary at infinity. Li & Cleall (2011) studied analytical solutions for advection and dispersion of a conservative solute in a one-dimensional, double-layered, finite porous media. The solutions have been applied to five different scenario-combinations of fixed, fixed flow and zero concentration gradient conditions at the input and output limits. These applications have been verified against numerical solutions from a finite-element-based approach and an existing closed-form solution for double-layered media, and an excellent correlation was found in both cases. Pérez Guerrero et al. (2013) have extended the Duhamel theorem, which was originally established for diffusion-type problems, to the case of advective-dispersive transport subject to transient (time-dependent) boundary conditions. As the analytical solution of the advection-dispersion solute transport equation remains useful for a large number of applications in science and engineering, several studies have concentrated in this direction. Berger et al. (2017) proposed to solve the advection-diffusion differential equation through two numerical schemes, the Scharfetter–Gummel and Crank–Nicolson approaches, whose efficiencies were investigated for both linear and nonlinear cases.

As the general solution of the advective-dispersive problem was largely discussed in the literature, applications in pipe networks are piecemeal due to the complexity of replicating analytical solutions for each pipe in the network with concatenating boundary conditions.

After a short discussion of the proposed modelling approach, the present study focuses on the following weak points of the state-of-the-art distribution network water quality models: the use of a simplified numerical transport model, the inability to consider dispersion, and the absence of experimental validation in low Reynolds regimes.

According to this aim, contamination experiments were performed using a conservative tracer in the laboratory water distribution network of University of Enna ‘Kore’ (Italy), and a field validation of the real water distribution network in Zandvoort (Netherlands)) was re-evaluated. The results of a state-of-the-art advective model and a diffusive-dispersive model was calibrated against laboratory and field data (Piazza et al. 2018).

The modelling analysis was carried out using the state-of-the-art advective EPANET model (standard advective model) and using an upgraded version of the model (diffusive-dispersive-advective model), including diffusion and dispersion equations proposed by Romero-Gomez & Choi (2011).

The advective phenomenon is generally linked only to the velocity of the fluid, as it involves the simple rigid translation in the space and time of the contaminant. Diffusion-dispersive phenomena are instead due to concentration gradients, which cause the migration of contaminants from high concentration areas to low concentration ones. Romero-Gomez & Choi (2011) realised that the presence of the solute trailing long after a tracer pulse has passed a fixed downstream position reveals that the velocity of dispersion towards the end of the pulse is stronger than the velocity near the front of the pulse. This result occurs because the low speed regions close to the wall strongly hinder the transport of the solute due to the non-slippery boundary condition and such condition differently applies to dispersion upstream and downstream from the contaminant injection.

• is the dimensionless pipe length that defines the location of solute migration, L, with respect to the pipe diameter, d;

• is the Reynolds number that accounts for the mean flow velocity's (um) geometric dimensions (d) and fluid conveying properties (kinematic viscosity, ν);

• is the Schmidt number that accommodates the solute properties (solute diffusion coefficient, D_AB);

• is the time, which is defined as the ratio between the location of solute migration, L, and the flow velocity (u_m).

The models were applied to the experimental network of Enna University – UKE (see (De Marchis et al. 2016; Piazza et al. 2017) for details), in which the network is operated as a single long loop (Figure 1) and to a real operational water network in the Netherlands, as described in Blokker et al. (2010).

The laboratory network is a closed water supply distribution network that is made up of three loops (only one was open during the experimental analysis), 10 nodes and 11 pipes, which have the characteristics of DN 63 mm, thickness 5.8 mm and length approximately 45 m and are arranged in nearly horizontal concentric circles with radii of 2.0 m; the network is supplied by four tanks, which can store up to 8 m³ of water. The supply tanks are connected to a group of four pumps and then to an air vessel in order to stabilise pressure. The pumping system behaves thus as a constant load tank, keeping the pressure constant and equal to a pre-set value between 1 and 6 bar and having a tolerance of 0.05 bar; the pumping system varies the speed of the pumps.

The system flows in the pipes are monitored by five electromagnetic flow meters installed in certain sections (trunks 4–5, 6–7, 7–8, 9–10 and 11–12). All the meters have sensitivity equal to 1 L/h, maximum flow equal to 5,000 L/h and measurement error equal 0.1% independent on flow. Pressure cells and multi-jet water meters are present in each node. Outflows in nodes are monitored with multi-jet turbine flow meters, compliant with EU MID directive, with nominal flow equal to 2,500 L/h (Q3 according to MID), starting flow equal to 4 L/h, errors lower than 5% under 40 L/h and lower than 2% in other cases (with flows between 40 L/h than the maximum 3,130 L/h). Additionally, WiFi real-time remote-controlled conductivity probes (Figure 2) were positioned at each node and connected with all the monitoring appliances to a central computer, which was also able to regulate flows supplied to the users by means of remote-controlled valves. The conductivity probes have an accuracy ranging from 1,300 μS to 40,000 μS and return as output the conductivity values in μs, total dissolved solids (TDS), and salinity. The salinity values are derived from the practical salinity scale (PSS-78). Further details about the laboratory network can be found in De Marchis et al. (2016). Experiments were carried out varying the outflows in nodes 6, 7, 9 and 11 (maintained constant in each experiment) with flows ranging between 80 L/h and 400 L/h in order to obtain different flow regimes in pipes. No leakages were applied in the present study.

The model was applied considering the different flow regimes and the maximum, minimum and average value of the dimensionless travel time (T) was calculated, assuming the conservative solute (NaCl) with D_AB = 1.2 × 10⁻⁹ m²s⁻¹. Table 1 shows the maximum, minimum and average values calculated for the UKE network.

The real network is situated in the town of Zandvoort, in the northwest of the Netherlands, along the sea (near Haarlem). The network was built in the 1950–1960s and consists of 5.7 km of Ø100 mm lined cast iron pipes and 3.5 km of Ø100 mm PVC pipes; it supplies 1,000 homes, two hotels and 30 beach clubs. The area is supplied from one point with a fixed head through a booster pump; there are no tanks in the network (Figure 3). Inflows are monitored with an electromagnetic flowmeter with similar characteristics to those used in laboratory experiments and supplied flows are monitored with turbine flow meters compliant with MID directive (maximum error lower than 5%), variable depending on water meter age, diameter and installation. The use of water in the network was determined by the historical flow patterns at the stimulation station, measured by the Provincial Water Company Noord-Holland (PWN), and the domestic water demand accounted for 70% of the total demand. Drinking water is distributed without any disinfectant, as is common in the Netherlands. A tracer study with NaCl was performed between 2 September and 20 October 2008.

Solute transport monitoring was enabled by dosing sodium chloride (NaCl) within a booster location nearby the network inlet, raising the electrical conductivity (EC) from about EC = 57 mS/m without dosage to about EC = 68 mS/m. Short intermittent contamination events (3 h) were performed with an inter-event time of 20 h. The contamination was carried out for seven weeks, but here we reported only a few days. EC values were measured at four locations in the system, and the water age was determined from the EC. Two models were constructed that are distinguished by demand allocation: Model_TD (top-down) and Model_BU (bottom-up). The first was allocated to all demand nodes with a correction factor DMP. This correction factor is the base demand and it was assigned according to the demand category. The bottom-up demand allocation was done with the use of the end-use model SIMDEUM that considers a stochastic water demand pattern. (Blokker et al. 2010). As leakage in the Netherlands is generally very low (2–4%) (Beuken et al. 2006; Geudens 2008), no leakage is assumed in this network.

In the upgraded version of EPANET model, the equations initially proposed by Romeo-Gomez & Choi (2011) were adopted concerning the dispersion/diffusion phenomenon; the results modelled were compared with the experimental data obtained by contaminating the water network in the Netherlands.

Initially, the simple advective approach (standard advective EPANET model) and the advective-diffusive approach (based on the Romero-Gomez and Choi formulation) were applied to the laboratory network and to the real network in order to compare experimental results and the numerical analysis. The results show the potential impact of diffusion/dispersion processes with respect to the water quality of the distribution networks (Figure 4).

The simplified advective model does not need calibration, as the quality module is a consequence of the hydraulic model parameters that were considered fixed in the present study. The diffusive-dispersive-advective approach requires the calibration of Equations (1)–(4) and specifically of the forward and backward dispersion coefficients, Ef and Eb.

The calibration process was based on a simple trial and error procedure that was performed with the aim of maximising of the Nash–Sutcliffe (N-S) convergence criterion over the other measured and simulated concentrations.

Figure 5 compares the experimental and numerical results (with and without dispersion) obtained by contaminating the UKE network (having opened every branch of the network) at node 6 with sodium chloride for a duration of 12 min and a mass of 370 g, leading to a constant concentration of 3,700 mg/L. In addition, the effect of the Reynolds number on the diffusive-dispersive phenomenon was evaluated. The results shown refer to nodes 6 (a), 7 (b) and 8 (c) of the UKE network.

The two most relevant parameters (the backward and forward dispersion parameters Eb and Ef) were calibrated at 2.6 and 3.1, respectively. The figure shows three cases in which the flows are laminar, in transition or turbulent between laminar and turbulent flow.

In Figure 5(a), the advective and dispersive behaviour coincide perfectly, since we are in the Turbulent Flow Regime and the Reynolds number is equal to 6,168. In this case, the dimensionless travel time assumes a value of 7.49 10⁵. It should be noted that in Figure 5(b), for the Transition Flow Regime with Re = 3,600, the agreement between the experimental and diffusive-dispersive-advective results is good. The advective model is not able to represent the process overestimating peaks and misses the delay of the concentration peak in time. The N-S criterion in this simulation was equal to 0.96 for the advective-dispersive-diffusive approach and 0.78 for the simple advective approach. In this case, considering a lower Reynolds number, we note an increase in the value of the dimensionless travel time equal to 1.285 10⁶.

Figure 5(c) and 5(d) show an experiment in which the Flow Regime is Laminar and Re is approximated to 1,500. The diffusion-dispersion phenomenon is more evident, and the differences between the two modelling approaches are large both in terms of peak estimation and peak delay. The N-S criterion in this simulation was equal to 0.74 for the advective-dispersive-diffusive approach and negative for the simple advective approach, demonstrating that this last method is unsuitable to analyse low-velocity flows in pipes. In Laminar Flow Regime, we note that the dimensionless travel time value assumes the maximum value equal to 2.998 10⁶.

For Re > 20,000, the advective approach overestimates a peak concentration of less than 10%; differences between the two models and the experiments become irrelevant for Re values higher than 35,000. The value of the dimensionless travel time, for Turbulent Motion Regime, assumes the minimum value equal to 2.31 10⁴.

The laboratory experimental investigation allowed for the identification of the values of the Reynolds numbers, for which dispersion and diffusion play a relevant role in contaminant propagation. The experimental campaign also demonstrated that advective simplification is acceptable when turbulent flow is well established. After the laboratory application, the two modelling approaches were applied to the real case study provided by a real water distribution network in the Netherlands.

Figure 6 shows the numerical results regarding water age, which were obtained by comparing the advective solutions determined in the study performed by Blokker et al. (2010) and the advective-dispersive-diffusive approach presented in this paper for the monitoring station located at Burg, Fennemaplein.

The comparison of the two approaches presents a substantial agreement in terms of water age. This is probably because the water age computation algorithm is based in both cases on a complete mixing approach at the nodes, and dispersion only plays a modest role in the evaluation of transit time in the pipes.

If the objective of the study is the evaluation of water detention in the system, the use of the complete advective-dispersive-diffusive approach is not justified, as the difference is limited to an average of −12% and has peaks around −20%. However, if the analysis aims to the investigation of contaminant concentrations, the two approaches agree in term of peak concentration and diverge significantly in terms of persistence of water contamination.

Figure 7 shows the results of the simulation and the electro-conductivity (EC) measurements during a contamination experiment on 3–4 September 2008. The NaCl-pulse at the inlet location into the network was from 3:40 to 6:40 on the first day. The pulse reached the different locations, depending on the transfer time, between the morning of 3 September and the morning of the next day. In particular, this pulse arrived on 4 September at the location Sterflat Friedhoffplein after 9:00 in the morning. For this reason, the figure shows two peaks (the second due to the next contamination experiment) in all locations apart from Sterflat Friedhoffplein. The figures clearly show that the advective approach tends to provide a much shorter contamination with a much smaller mass of the contaminant reaching the user. The diffusive approach provides a more realistic distribution of contaminant concentrations, respecting peaks and providing a better estimation of contaminant masses. Once again, calibration provided better results for the advective-dispersive-diffusive approach, providing higher values of the N-S criterion (respectively, 0.79 and 0.88 in the two time frames presented in Figure 7(a) and 7(b)) obtained for a backward dispersion coefficient equal to 2.42 and a forward dispersion coefficient equal to 2.60. Using the same parameter values in the other two locations, the N-S criterion was equal to 0.68 and 0.71 in the two other locations, represented in Figure 7(c) and 7(d).

Lastly, the Reynolds number for the entire network in Zandvoort was determined. It is observed that within 24 h, the speeds are variable depending on the diameters present in the network and the flow rates taken by the users. This variability makes the global network flow regime oscillate in a range that assumes values of the Reynolds number with a minimum equal to 198 and a maximum equal to 99,684 (respectively, laminar and fully turbulent). Applying the numerical model in this network, we have calculated the minimum, maximum and average value of the dimensionless travel time equal to 2.42 10³, 4.815 10⁶ and 1.21 10⁴, respectively.

The previous application showed the reliability of the advective-dispersive approach when dealing with low-velocity water distribution networks; the advective approach showed a general overestimation of peak concentrations and an underestimation of contaminant persistence in pipes. As presented in this paragraph, the two modelling approaches were compared in the identification of optimal water quality sensor locations. The approach was as follows.

Considering the UKE laboratory network, the optimisation problem was solved using the NSGA-II algorithm for positioning a set of three sensors. Three objective functions were considered, namely Detection likelihood (F_1), Detection time (F_2) and Redundancy (F_3):

• F_1: Detection likelihood, i.e. the probability of a sensor configuration to detect the contamination;

• F_2: Detection time, i.e. the average time passed between contamination and detection in the 200 simulations constituting each individual;

• F_3: Detection redundancy, i.e. the probability that the contamination is detected by two sensors within 20 min.

The objective functions were slightly adapted from those presented in Preis & Ostfeld (2008) in order to comply with the smaller dimensions of the analysed network and they were equally weighted in the selection of the optimal sensor location. Contamination is randomly setup in nodes (random contaminant mass, contamination duration and contamination node). User demands in all nodes were fixed and equal to 2.5 L/min. The inlet head was fixed as well to 3.5 bar. At the same time it was an experiment, in which the contamination was performed in node 6 with a duration of 12 min and a mass of 370 grams (leading to a constant concentration of 3,700 mg/L).

The optimisation problem was run with 200 generations, each made by 50 individuals, and the mutation and crossover probability are equally set to 20%.

The optimisation problem using the advective model provides two possible configurations for the sensors positioning, as presented in Figure 8. Table 1 shows the characteristics and performance of the optimal configurations. A single configuration (based on sensors in nodes 6, 7 and 10) is able to maximise two objective functions (F_1 and F_2). To maximise function F_3, node 8 should be included in the monitoring campaign. With the use of three sensors in the optimal positions, 92% of the contamination episodes may be averagely detected (maximum value of function F_1), and 56% may be detected by at least two sensors within 20 min (maximum value of function F_3). The average detection time is optimally equal to more than 10 min. Figure 8 shows the Pareto fronts obtained for the three objective functions.

The optimal solutions considering dispersion differ from the simple advective case by considering both the location of optimal sensors and the objective function values (Table 1). Each objective function is optimised by a different sensor configuration, and only node 6 is present in all cases Table 2. Generally, the optimal configurations do not privilege central nodes of the system (like in the advective case), probably because concentrations decrease rapidly and move from the contamination node to the others, thus requiring the sensors to be widely distributed to increase detection likelihood. Pareto fronts are presented in Figure 9. The maximum detection likelihood (F_1) is provided by a configuration containing one external node and two internal ones (Figure 10).

The optimal locations are not contiguous; this distribution is probably related to the sharp attenuation of concentrations moving from the contamination node, requiring the sensors to be more distributed in the network. The detection time (F_2) does not change significantly in the two optimisation exercises, and the two optimal configurations are superimposed with the exception of one node. Additionally, the optimisation of redundancy function (F_3) shows different results in the two exercises: thus, including dispersion, the nodes providing the best result are widely distributed in the network, revealing the importance of contaminant attenuation in its detection.

The analysis showed how dispersive and diffusive processes are relevant in the simulation of solute propagation in water networks. The importance of such processes decreases with the presence of turbulence in pipes, and these processes can be considered negligible once the Reynolds number reaches values higher than 35,000 and turbulence is stable and established.

The laboratory analysis was able to highlight the impact of dispersion and diffusion in controlled conditions, but the tests on a real network added some interesting additional elements, supporting the use of more complex modelling approaches:

• In real networks, Reynolds numbers may vary significantly during the day, and diffusion-dispersion may result by having variable impact on contamination.

• Even if peaks are better represented by the complete approach, such a model is still incapable of representing the complexity of the recession limb of the pollutograph.

• The need for calibration of the diffusion-dispersion coefficients requires the availability of real data for calibration, even if the parameter values are not characterised by significant variations in the two analysed case studies.

Considering that one of the most frequent applications of numerical network models is the identification of optimal water quality sensor locations, some general considerations can be drawn, as follows:

• Detection likelihood decreases with the Reynolds number because contamination more often passes the sensor with concentration lower than the sensor sensitivity.

• Redundancy decreases and detection time increases because of the temporal lag that is clearly visible in Figure 5 depending on the importance of dispersive processes.

• When the network is characterised by low velocities, neglecting dispersive processes takes as a consequence an incorrect estimation of concentration peaks and persistence.

• Using advective-only approaches, the optimal sensor locations are mainly concentrated in the central part of the network; however, using the dispersive-diffusive approach, optimal locations are equally distributed between the central part of the network and the periphery.