Interdependence of flow and pipe characteristics in transient induced contamination intrusion: numerical analysis

It is highly probable that water quality declines during water hammer since the negative pressure of the fluid is able to suck pollutants from possible leaks to the distribution system. Contaminant intrusion is usually the main cause of water quality loss that may endanger clients' health. In the UK, Hunter et al. (2005) reported results from a postal questionnaire, finding a strong association between loss of water pressure at the home tap and the incidence of diarrhea, and they estimated that the cost of such illness could exceed £100 million/year in England. Three issues to be discussed in more detail provoke this phenomenon: contaminant sources, driving forces (low/negative pressure), and pathways (Fox et al. 2015).

Contaminant source is usually available in urban areas such as sewer lines, garbage disposal areas, leaking underground petroleum storage units, etc. Consequently, buried water mains are often exposed to a wide range of chemical and microbial contaminants. Schuster et al. (2005) analyzed 288 waterborne outbreaks in Canada occurring during 1974–2001, and identified a variety of pathogens in the outbreaks.

Studies show that contaminants often intrude into pipe systems through a variety of pathways such as submerged air valves, leak points, repair and installations, faulty seals, joints, and service connections (Starczewska et al. 2015). Besner et al. (2010) showed that bacterial contamination has been found in soil and water surrounding pipes, and hence leaks, which cannot be easily managed, could provide a pathway for contaminant intrusion. Figure 1 displays the images of semi- and completely submerged air valves provided by Besner et al. (2011) which provide examples of contaminant pathways.

Negative gauge (or relative) pressure of pipe flow is the driving force which sucks pollutants to the main flow. Negative pressure is likely to occur due to pump operation, pipe replacement, valve closure/opening, demand variations, etc. (Collins et al. 2012; Meniconi et al. 2015). Lindley (2001) indicated an adverse pressure gradient as a susceptibility condition for contaminant intrusion in 62.8% of the outbreaks in USA during 1971–1998. Kirmeyer et al. (2001) and Mays (2000) provided a detailed list of causes for low/negative pressure events in the normal operation of a water distribution system (WDS) including: variation of system demand, firefighting, pump failure, valve closure, flushing and ground water level increase. Gullick et al. (2005) also carried out pressure monitoring for 43 sites in eight WDSs. Over a period of 4,640 days, they reported 21 negative pressures that lasted less than 3 minutes, mainly caused by sudden pump shutdowns.

Most previous studies have suggested the Eulerian–Lagrangian method to investigate the propagation of the intruded contaminants during transient conditions. Some details on the literature of the method are provided here.

Fernandes & Karney (2004) were among the first who predicted the behavior of contamination intrusion from the leak point caused by water hammer. They examined how the concentration of trapped chlorine is altered, and they also further investigated the amount of contamination in the vicinity of the leak site as a result of both inflow and outflow through the leak. Basha & Malaeb (2007) introduced and emphasized the application of a combined Eulerian–Lagrangian method whose basis is a constant velocity for the movement of contaminant parcel. Mora-Rodríguez et al. (2012) proposed a combined 1D-3D (respectively used for the transient flow simulation and for the intrusion quantification) procedure to compute total volume of intrusion which is of great importance if the external soil affects the discharge through the orifice. Mansour-Rezaei & Nasser (2013) studied the suction and release of leak contamination during transient pressure oscillations. They utilized the method of characteristics (MOC) to solve the hydraulic equations in conjunction with a Lagrangian method to simulate the emission of pollutants. Considering the possible existence of soil around the leak site, a model for estimating the incidence of contamination from the surroundings of the leak site into the tube is provided by Mansour-Rezaei & Nasser (2013). Payesteh & Keramat (2017) recently exploited the Eulerian–Lagrangian method to conduct a sensitivity analysis of the hydraulic parameters affecting the amount of intrusion in a reference reservoir-pipe-valve system with a leak. They concluded that the magnitude of the negative pressure at the leak point is the most important factor.

Meniconi et al. (2011), Mora-rodríguez et al. (2012) and Fox et al. (2014, 2015) experimentally revealed that when a water hammer occurs in a network, the negative pressure wave sucks the contamination around the leak site into the distribution system. Subsequently, the intruded volume travels towards the downstream of the leak site. The laboratory results of Fontanazza et al. (2015) showed that the amount of contamination due to the permeation mechanism in semi-filled pipes is relatively higher than the contamination induced by transient flows. They also showed that in pressurized pipes the amount of intrusion during transient flows is directly governed by the magnitude and duration of the negative pressure at the leak position. Jones et al. (2018) measured the contaminant intrusion rate under specified initial conditions in terms of discharge and steady-state pressure.

The effect of each parameter on intrusion was separately studied by Payesteh & Keramat (2017), nevertheless, it draws no conclusion regarding interconnection between parameters. Therefore, in this research, the simultaneous effect of two or more parameters during transient flows in pressurized pipes is assessed by means of numerical experiments. In fact, the assessment of the effective factors (pipe length, main flow rate, leak size and location, and pipe material) in the amount of contaminant intrusion in either the laboratory or the field requires extensive time and costs, hence this study aims to numerically address the parameters' contributions. The Eulerian method of characteristics model transient pipe flow and the Lagrangian solution of the advection equation capture total Volume of Contaminant Parcel (VCPt) penetrating through a leak. The VCPt establishes as a criterion to compare various transient scenarios and the interconnection between key parameters. The numerical investigations carried out in a short and long transmission line comprise 72 different hypothetical transient scenarios to quantify the amount of intrusion due to negative pressure wave at the leak point.

Equations (1) and (2) are solved using the Method of Characteristics (MOC). The transient state at the leak node contains four unknowns comprising flow rate upstream, downstream and through the orifice (leak) and pressure head at the leak node (Brunone 1999).

The constitutive equations to estimate the concentration of the contaminated water and its volume in the vicinity of the leak node and total volume of intrusion are illustrated in this section. The Lagrangian approach of the advection equation, which resolves complexities of a non-uniform advection speed and hence time step mismatch between advection and transients models (Fernandes & Karney 2004), is further elaborated.

In the widely established Lagrangian method, as revealed in Equation (11), the contamination is modeled as a parcel with a varying volume but constant concentration that moves with the velocity of the main flow which itself is governed by the continuity and momentum equations during transients. Consequently, the transient nature of the water hammer may change the direction of flow in the vicinity of the leak location, thus leading to intrusion and/or extrusion of contamination.

As seen in the schematic of a leaky zone in Figure 2, during the steady state flow the water leaves the system due to the positive pressure head difference between the pipe and the outside. The leaking flow can provide favorable conditions for the growth of microbes around the leak site. The resulting moisture can also dissolve the soil chemicals around the leak site and facilitate the suction of hazardous materials into the pipe. When the negative pressure wave of the water hammer arrives, the contamination around the leak area is sucked into the pipe, as seen in Figure 4. With the onset of the next cycle of the water hammer, the positive pressure wave at the leak pushes part of the contamination out of the pipe and also to the upstream of the leak. By evacuating the water hammer wave, the reciprocating movement of the contamination ends around the leak site and ultimately a part of the contamination is trapped in the pipe (Fernandes & Karney 2004; Meniconi et al. 2011). The number of trapped contaminant parcels in the pipe depends on the number of back and forth cycles of the water hammer.

Because the main aim of this study is to compare the volume of intrusion in various hydraulic scenarios, possible (partial) extrusion of the contaminant parcel is neglected during a transient cycle. This simplifying assumption is obviously conservative and in fact corresponds to a worse case as it estimates the contamination as larger than the actual one.

To examine the illustrated numerical method for the estimation of the intrusion volume, the laboratory experiment conducted by Jones et al. (2018) was numerically simulated and then the computations were verified against experimental results. In another survey, the transient induced intrusion in viscoelastic and elastic pipes were compared in order to assess the significance of the retarded behavior of viscoelastic pipes.

The large-scale laboratory facility at the University of Sheffield (Fox et al. 2017) comprises 140 m of 50 mm internal diameter medium-density polyethylene (MDPE) pipeline with the specifications provided in Table 1. A circular orifice of diameter 1.55 mm is situated about 75 m downstream of the supply reservoir. The leakage orifice is enclosed within a 400-mm length of larger 380-mm diameter outer pipe creating a volume surrounding the orifice, as sketched in Figure 5. To quantify the intrusion of fluid contained in the outer pipe into the main pipeline, a clear riser pipe is connected to the outer pipe such that the level in this riser allows intrusion to be measured. In the following, the experimental data of intrusion in this facility (Jones et al. 2018) are incorporated to verify the explained numerical formulations.

The VCPt quantity in Equation (17) is directly determined by the leak discharge history, which is governed by the water hammer model. Therefore, the verification herein consists of water hammer modeling in the leaky pipe system of the Sheffield experiment. In the numerical solution, use is made of a Discrete Vapor Cavity Model (DVCM) in a viscoelastic pipe including unsteady friction (Soares et al. 2009; Brunone & Berni 2010; Keramat et al. 2010). Almost all transient tests described in Jones et al.’s (2018) paper experience column separation, so the DVCM which is of acceptable accuracy (Bergant et al. 2006) is adopted herein. The generalized Kelvin–Voight model (Keramat et al. 2012) with some typical coefficients provided in Keramat & Haghighi (2014), and the Zielke–Vardy model for UF are included in the solver to improve its performance; in Keramat & Tijsseling (2012) the complete form of the cited mathematical model and corresponding numerical solution are illustrated.

Figure 6(a)–6(d) depicts the results for four executed transients generated by the upstream valve closure including: (a) initial flow rate Q₀ = 1 L/s, reservoir head = H₀ = 20 m; (b) Q₀ = 2 L/s, H₀ = 10 m; (c) Q₀ = 2 L/s, H₀ = 30 m; (d) Q₀ = 4 L/s, H₀ = 20 m. These results are compared with the measured pressure heads reported by Jones et al. (2018) (Figure 3 of their paper). The agreement between the two sets of figures is acceptable, and the slight discrepancy between the two sets of graphs may be attributed to possible uncertainties of input parameters or measurements.

In the other investigation conducted by Jones et al. (2018), the variations of intrusion volume subject to different values of reservoir pressure head for the same steady state flow rate Q₀ = 2 L/s was studied and the results were compared with the experimental measurements. In the experimental study of Jones et al. (2018), the collected data points are accompanied by the regression line (VCPt = –1.23 H₀ + 81.43) with a favorable coefficient of determination (R² = 0.947), which means that the data points agree reasonably with the fitted line. Figure 7(a) depicts the comparison where the continuous line indicates the numerical results and the dashed line represents the aforementioned regression line (Jones et al. 2018). Figure 7(b) provides another comparison in which the reservoir pressure head is kept constant (H₀ = 20 m) while the flow rate alters to capture the variations of the intrusion volume. Likewise, the computations of intrusion (continuous line) are compared with the experimental regression line (dashed line) reported by Jones et al. (2018). As Figure 7(a) and 7(b) demonstrates, the fitted lines of experimental observations and the numerical calculations of intrusion volume are in acceptable agreement.

This numerical survey serves as an investigation into the significance of the retarded behavior of the viscoelastic pipes. With respect to elastic pipes, the viscoelastic ones have two crucial characteristics: their smaller pressure wave speed and a lag in their response to an applied load. The former is easy to study by means of regular water hammer models in elastic pipes but the latter requires a stress-strain analysis to include the retarded behavior of these materials in their dynamic response (Covas et al. 2004, 2005; Keramat et al. 2012, 2013; Pezzinga et al. 2016), which is addressed in this section. To this aim, the simulations in the previous section that consider the viscoelastic effects are now repeated using the elastic model and the results are provided in Figure 8. The dash-dot lines correspond to the results of the elastic water hammer model and the continuous lines indicate those of the viscoelastic model. The results support the deduction that consideration of the retarded behavior contributes to lessen the intrusion. The reduction on intrusion amounts can be justified by the function of viscoelastic materials in attenuating the water hammer pressures drops, which form the main driving force of intrusion in this study. Because neglecting the creep behavior of viscoelastic pipes corresponds to the worst case of transient-induced intrusion, the evaluations of intrusion in the next section neglect this property and consequently, the viscoelasticity is accounted for solely by its low wave speed.

In order to provide a comprehensive understanding of the contamination intrusion phenomenon in various conditions, 72 different cases were studied. In each one, different leak locations and sizes, different fluid velocities due to different reservoir pressures in elastic and viscoelastic pipe for short and long pipeline were investigated (Table 2). In test 1, pipeline length L = 2,300 m, pipe diameter D = 600 mm and valve pressure head (a typical model of suburban water transmission mains). In test 2, L= 540 m, D = 108 mm, and the rest is as previous (as a sample of urban pipelines). For the two tests: (i) three positions for the leak (x_l) including the middle, one third and two third from upstream, (ii) two leak size ratios being leak area over the cross-sectional area of flow = 0.01 and 0.001, (iii) three fluid velocities representative of conventional speeds in urban water transmission systems including U₀ = 0.5, 1 and 2 m/s; (iv) two pipe materials being steel and PVC (Poly Vinyl Chloride), altogether making 72 scenarios studied (Table 2). Note that for viscoelastic pipes, only their small wave speed (Table 2) is applied in the simulations and their creep behavior is neglected as discussed in the previous section.

In each case, the numerical model of the water hammer allows for evaluation of transient flow rate at the leak which is then exploited to compute VCPt based on Equation (17). For example, Figure 9 depicts the time series of pressure and discharge at the leak location for test 2 with U₀ = 0.27 ms^–1, x_L = 270 and ⁠. The dashed line in Figure 9(b) corresponds to the zero leakage so that the flow rate below this line accounts for intrusion which occurs in four cycles in a 10 s time span. It is worthwhile to note that according to Figure 9(b), the whole area between the plotted line for Q_L and the zero line (the quantity ⁠) is apparently a positive value meaning that, overall, a volume of water exits the system during the water hammer. This volume may consist of entirely fresh water, entirely contaminated water or a partial combination of both. If the outflow is either fully or partially fresh water, it is most likely that the contaminated parcel is trapped in the pipe (Fernandes & Karney 2004).

The procedure to compute VCPt is implemented for the 72 scenarios and the results are summarized in Table 3. The discussions associated with the findings of this table are provided in the following.

In test 1, which corresponds to a long pipeline, the effect of viscoelasticity in declining intrusion is prominent, making (neglecting the creep behaviour) which represents a very small amount of contaminant intrusion. In the second test, however, this ratio reaches ⁠, meaning that in short pipes the viscoelastic property of the pipe wall is less effective on the reduction of intrusion.

As seen in Equation (17), the amount of intrusion mainly depends on the magnitude of the negative pressure and its duration at the leak location. Consequently, viscoelastic pipes have smaller intrusion due to their smaller wave speed. The significance of duration of down-surge is clearer when the intrusions of long and short pipes are compared; the results of which were provided earlier for the two tests.

In order to investigate the interaction between velocity and pressure of fluid on the entry of contamination in all scenarios, two non-dimensional quantities are used to plot the numerical results. One adopts the steady state fluid velocity and pressure head difference and the other allows for total volume of contaminant parcel ⁠, where V_pipe is the volume of the pipeline. The non-dimensional numbers are defined on account of the Joukowsky pressure and the fact that intrusion volume increases by the pipe length and flow cross-section.

The general trend regarding the amount of intrusion is expected to be according to: (i) Joukowsky's pressure rise, so that higher wave speed and/or initial velocity bring about higher intrusion (Figure 10(a) and 10(b)), and (ii) steady state pressure head at the leak (dictated by the upstream reservoir head), whose higher quantities produces less intrusion simply because, in Equation (3), decreases. However, several counterexamples have been found among the 72 cases. For instance, in test 2 for and viscoelastic pipe, the intrusion of the case U₀ = 2 m/s is less than that of U₀ = 1 m/s (Figure 10(c)). Another example of unusual interdependence among affecting parameters is observed with changing reservoir pressure (while the same pressure at the valve is maintained). In a velocity range from 0.5–1 m/s, increasing reservoir pressure leads to rising contaminant intrusion (Figure 10(a), 10(b) and 10(d)), while for velocities higher than 2 m/s, the pressure of the reservoir has a suppressing effect on negative pressure at the leak, thus reducing intrusion (Figure 10(c)). These findings are in agreement with the laboratory results of Jones et al. (2018).

Considering Equation (17), two aspects clearly contribute to increase VCPt: the magnitude of the negative leak discharge (or the magnitude of negative pressure at the leak) and the duration (denoted by d) for which negative discharge at leak occurs. These two components are evident in Figure 11(a) and 11(b) which shows how they alter by the leak size variation in both tests. In Figure 11, the horizontal axis is the dimensionless time, which is scaled by the transient wave travel time from the valve to the reservoir, and the vertical axis represents the dimensionless discharge given by transient leak flow rate over its steady state quantity. Regarding the inflow duration d to the main pipe, the role of leak size ratio δ is of great significance. Specifically, as seen in Figure 11(a) (test 1) and Figure 11(b) (test 2), for δ = 0.001 (blue curve), d is the summation of several inflows occurring at a number of water hammer periods, while for δ = 0.01 (red curve) only two rarefaction waves at the leak location determines d.

According to Figure 12, for studying the leak location (in the large leak case δ = 0.01), long duration times are accompanied by lower magnitudes of negative discharge and short duration times correspond to higher magnitudes (black lines) in both tests. The opposing trend of these two key quantities reveals the maximum at which the highest amount of intrusion occurs. This pattern is in fact due to the interaction between the main water hammer wave (determined by the valve maneuver and the reservoir head) and leak induced waves. This complicated interaction eventually causes the maximum intrusion to be formed at the middle of the pipeline in both cases of leak location (x_L/L = 0.5, Figure 12(b) and 12(d)). For the small leak case (δ = 0.001), the interaction between the valve and the leak-induced waves are less dominant in d, that is to say wave reflections from the leak are negligible so that the amount of intrusion is simply dictated by the main water hammer waves (generated by the valve). Therefore, the duration d is mainly governed by the leak distance from the upstream reservoir meaning that distant leaks (from reservoir) result in higher d and hence more intrusion (see Figure 10(a) and 10(c)).

According to Equation (3), it is expected that the increase in the leak area will increase the contamination intrusion to the pipe, but the results manifest a number of counterexamples. For instance, Figure 13 shows that at low fluid velocity (e.g. 0.5 m/s) with increasing leak diameter, the amount of contamination entered to the pipe decreases, whereas at high velocities (say 1 and 2 m/s) with increasing leak diameter, contamination sucked into the pipe rises. This behavior is justified with the significant co-relation between the leak size and steady state velocity, which form the influential components of intrusion volume. As a result, one cannot deduce a general trend for intrusion quantification in terms of the leak size or fluid velocity separately. In fact, at a low fluid velocity, the increase in leak diameter leads to increased outflow from the pipe during the first positive pressure wave, hence the transient pressure wave is rapidly damped. This results in a slight negative pressure in later cycles of the water hammer, which ultimately reduces the suction of the contamination into the pipe. Conversely, at higher flow rates, the effective parameter is leak diameter with which the amount of contamination entered to the pipe is altered.

This section serves as an overview of the quantification of the intrusion volumes in the provided case studies, and expresses how they relate or what they add to the literature.

The numerical and experimental studies in the literature of the intrusion problems fall into two research categories: (i) quantification of both constituent volume and its concentration (e.g. Fernandes & Karney 2004; Mansour-Rezaei 2013; Mansour-Rezaei & Naser 2013), and (ii) only the quantification of the intrusion volume (e.g. Fox et al. 2014, 2015, 2017; Fontanazza et al. 2015; Jones et al. 2018). The first branch often associates with the chlorine advection, dissipation and decay in pipes while the second commonly addresses the intrusion of the contaminated liquid through leaks in which a precise modeling of concentration is of less significance. In the latter group, diffusion and decay mechanisms are usually disregarded as they are taken into account when the estimation of concentration is of interest. In this research, diffusion is neglected but conversely the following two prudent assumptions apply: (i) possible extrusion of contaminants, which may occur when positive pressure prevails at the leak point, is neglected, and (ii) among the three scenarios of contaminant suction at the orifice site illustrated in Figure 4, the worst case in which the concentration after the intrusion is equal to that of the outside (Figure 4(c)) is assumed throughout the intrusion duration. Consequently, the Lagrangian approach (e.g. Basha & Malaeb 2007) estimates the length of the contaminant parcel which is VCPt/A where A is the cross-sectional area of flow. Considering the two types of intrusion problems, the corresponding case studies addressed in the literature are also different. For instance, Fernandes & Karney (2004) and later on Naser & Karney (2008) studied the variations of the water quality due to valve opening or closing in reservoir-pipe-valve systems (without leak) whereas Fox et al. (2017) and Jones et al. (2018) focused on the quantification of intrusion volumes sucked through a leak. The numerical investigations of the current research aim to give insight into the prediction of the intrusion amounts occurring due to valve closing in a reservoir-leaky pipe-valve system subject to various system parameters.

In the studied reservoir-leaky pipe-valve system, the pipe material, the initial velocity (or the reservoir pressure, assuming a constant head at the downstream valve), the leak size and the leak location, are treated as independent quantities whose variations affect the intrusion volume. The results manifest that the intrusion trends cannot directly be predicted by variations of these quantities. More specifically, a precise evaluation of the intrusion volume requires solving the transient hyperbolic equations in conjunction with the orifice relation. However, the intrusion duration (d) and the extent of the inflow rate are two important determinants to predict the variations of the intrusion volume with respect to the independent quantities. Fortunately, these two quantities can be estimated by approximate relations: the duration is evaluated by the water hammer period (4L/c) and the leak location (x_L), whereas the inflow extent is proportional to the effective leak size (A_e) and the Joukowsky pressure (H_res ± V₀c/g). In view of the effect of the fundamental water hammer period, one can infer that long pipes are prone to higher amounts of intrusion. The two terms in the Joukowsky relation have opposite effects on the intrusion volume meaning that increasing the first term enhances and the second term declines the intrusion, noting the challenge of the interdependence between H_res and V₀ in this relation addressed in Figure 10. The wave speed c also has contradictory impacts on the duration and inflow magnitude, thus making it difficult to predict whether its variation rises or drops the intrusion volume. Although the numerical case studies showed that, on average, pipes with lower wave speed (viscoelastic pipes) have less intrusion, a couple of cases showing higher intrusion for lower wave speed can be found in Table 3.

In this study, the Eulerian method of characteristics was employed to model transient flow and determine total Volume of Contaminant Parcel (VCPt) entering from the leak site. Seventy-two different hypothetical transient scenarios were considered to study the amount of intrusion due to negative pressure wave at the leak point. The key findings of the executed numerical experiments are:

• Viscoelastic pipes are greatly advantageous in intrusion reduction, especially in long transmission lines which have higher intrusion (due to their higher fundamental water hammer period). The ratio of intrusion (defined by intrusion volume in viscoelastic pipes to that in elastic pipes on average) for the pipe length 2,300 m is 0.027, while the ratio for pipe length 540 m is 0.496.

• It may seem that higher Joukowsky pressure or leak size produces more intrusion volume, but this is not always true and the numerical results demonstrate several counter-examples.

• Large leak sizes are prone to significant leak induced wave reflections, thus making the leak position zone of high intrusions quite unpredictable (case dependent).

• Small leaks cause less significant valve-leak-reservoir interactions of transient waves. This means that the amount of intrusion is dominated by valve-closure waves so that high intrusions are more likely when the leak is close to the downstream valve owing to longer duration of negative pressure.

• There is a positive correlation between the leak size and the intrusion amount if the velocity of flow is high (considering the defined dimensionless number), and a negative correlation is found for low velocities. Considering the orifice relation, higher velocity establishes a considerable negative transient pressure at the orifice, hence large orifices produce more intrusion. However, low velocity generates a small negative gauge pressure which may be suppressed by reservoir pressure (in the case of a large leak).

Italian MIUR funded this research, and the University of Perugia is acknowledged within the program Dipartimenti di Eccellenza 2018–2022. The financial support from Jundi-Shapur University of Technology is also appreciated.

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/hydro.2020.069.