Abstract

Contaminant intrusion in pipelines during transients is a remarkable mechanism, which leads to a decline in the quality of the contained water. The negative pressure of water hammer pressure waves is the trigger for the suction of pollution from the surrounding leak area, and hence deteriorating water quality. The volume of contamination intruded into the pipeline is investigated using mathematical and numerical modeling of the phenomenon. To elucidate this phenomenon in real pipe systems, the intrusion amount is estimated for 72 different scenarios including: two lengths of pipeline (i.e. short and long), three different leak locations, three different fluid velocities in the pipe, two leak diameters and two pipeline materials (elastic and viscoelastic). The results showed that the amount of intrusion in viscoelastic pipes was clearly less than that in elastic pipes, especially in long pipelines. The critical zone of high intrusion risk is identified close to the downstream valve for small leak sizes, nevertheless, it is difficult to estimate this zone in the case of large leaks due to significant interactions between nodal components (valve, leak, reservoir).

INTRODUCTION

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.

Figure 1

Pictures showing air valves as a pathway of contaminant intrusion: (a) Semi-submerged air valve; (b) completely submerged air valve (Besner et al. 2011).

Figure 1

Pictures showing air valves as a pathway of contaminant intrusion: (a) Semi-submerged air valve; (b) completely submerged air valve (Besner et al. 2011).

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.

Contaminant intrusion modeling

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.

Laboratory evidence of intrusion

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 aims of the study

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.

METHODOLOGY

Mathematical and numerical water hammer simulation

Governing equations

The following mass and momentum conservation equations are used for simulating fluid transients in elastic pipes (Chaudhry 2014): 
formula
(1)
 
formula
(2)
where x = distance along the tube, t = time, g = gravitational acceleration, D = pipe internal diameter, A = pipe cross-section area, Q = discharge, H = pressure head, f = friction coefficient, and c = pressure wave speed. Figure 2 depicts a schematic of a reservoir-pipe-valve system with a leak to investigate the contaminant intrusion phenomenon. The extrusion of fresh water and intrusion of contaminated liquid takes place through the leak on the pipeline during transients generated by the valve maneuver. The resistive effect of the external soil is neglected.
Figure 2

Schematic of the reservoir-pipe-valve system with a leak to investigate the contamination intrusion phenomenon.

Figure 2

Schematic of the reservoir-pipe-valve system with a leak to investigate the contamination intrusion phenomenon.

The leak discharge can be calculated by the orifice equation (Massey 1998): 
formula
(3)
where are the pressure head outside the leak (given due to the hydraulic conditions of the leak surroundings) and head inside the pipe (determined by water hammer solution), respectively, is the leak effective area, and is the discharge through the leak.

Initial conditions

The steady-state velocity and pressure head in the leaky pipeline constitute the initial conditions of the water hammer solution. Use of the energy equation between the reservoir and the leak and then the leak and the valve, including the Darcy–Weisbach head loss, allows for evaluating the flow rate at either side of the leak (Swamee & Sharma 2008): 
formula
(4)
 
formula
(5)
in which is the pressure head at the reservoir, is the elevation of reservoir, and are the lengths of the pipes before and after the leak, and are discharges before and after the leak, is the leak elevation, is the valve elevation, and is the pressure head at the valve.

Numerical simulation

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).

According to Figure 3, by combining the and equations: 
formula
(6)
 
formula
(7)
with the continuity equation at the leak: 
formula
(8)
and the orifice (Equation (3)), the following relationship is obtained: 
formula
(9)
where parameters , , and are evaluated at the previous time step (Chaudhry 2014).
Figure 3

Schematic of the leak, intrusion and the computational transient model.

Figure 3

Schematic of the leak, intrusion and the computational transient model.

Contaminant intrusion quantification

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.

Contaminant advection in Lagrangian approach

The contaminant transport is governed by the advection-diffusion equation, which is based on mass and Fick's law (Jaynes & Rogowski 1983). In the intrusion problem, the diffusion term is negligible compared to the convection term so that (Rossman & Boulos 1996): 
formula
(10)
in which u is fluid velocity, and is the concentration of the contaminant fluid. Assuming that and by using the chain rule, Equation (10) reduces to (Basha & Malaeb 2007): 
formula
(11)

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.

Figure 4

Various scenarios of suction at the leak site during transients.

Figure 4

Various scenarios of suction at the leak site during transients.

Concentration at the leak control volume

Figure 4 shows three possible scenarios of intrusion in terms of flow direction at either side of the leak node during the water hammer. With the assumption of complete and sudden mixing at the site of the leak node (Al-Omari & Chaudhry 2001), the following formula is proposed to compute the concentration of input contamination due to the water hammer: 
formula
(12)
in which is the concentration of the contaminated water at the leak node, is the concentration of the intruded flow from leak site, is the concentration of flow at the upstream of leak, is the concentration of flow at the downstream of leak. The derivation of this equation is shown in Appendix A in Supplementary Materials.

Intrusion volume

Considering Equation (12) and Figure 4, the concentration of the contaminant parcel (CCP) in the Lagrangian model is calculated by means of the following relationship (Mansour-Rezaei & Nasser 2013): 
formula
(13)
in which is given by Equation (12) and N is the number of time steps in one suction cycle of flow contamination due to the effects of transients. The volume of the contaminated parcel can be calculated as follows: 
formula
(14)
where represents times at which contaminants intrude from the leak in a suction cycle of contamination and is the time step of the hydraulic model.

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.

Equation (14) estimates intrusion during one suction cycle of the water hammer, therefore the total intruded volume (VCPt) is given by: 
formula
(15)
in which is Heaviside function defined by: 
formula
(16)
Equation (15) may be simplified as follows: 
formula
(17)
where M is the total number of time steps in the transient event. In Equation (17), the non-zero terms in the summation correspond to the time steps during which the flow rate at the leak (QL) is negative, meaning that contaminated water intrudes into the main pipe flow from the leak. The time duration for which the flow rate at the leak point is negative is one of the principal quantities in the estimation of intrusion which will be further discussed in the case studies.

Comparison with experimental data

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.

Table 1

Pipe system specifications of Jones et al.’s experiment (2018)

Internal diameter of pipe (mm) 50 
Pipe length (m) 140 
Orifice diameter (mm) 1.55 
Orifice distance from reservoir (m) 75 
Wave speed (m/s) 570 
Flow rates (L/s) 1, 2, 3, 4 
Reservoir pressure head (m) 10, 20, 30, 40 
Internal diameter of pipe (mm) 50 
Pipe length (m) 140 
Orifice diameter (mm) 1.55 
Orifice distance from reservoir (m) 75 
Wave speed (m/s) 570 
Flow rates (L/s) 1, 2, 3, 4 
Reservoir pressure head (m) 10, 20, 30, 40 
Figure 5

The laboratory facility of Jones et al.’s (2018) intrusion experiment. The figure is extracted from their paper.

Figure 5

The laboratory facility of Jones et al.’s (2018) intrusion experiment. The figure is extracted from their paper.

Assessment of numerical solution

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 Q0 = 1 L/s, reservoir head = H0 = 20 m; (b) Q0 = 2 L/s, H0 = 10 m; (c) Q0 = 2 L/s, H0 = 30 m; (d) Q0 = 4 L/s, H0 = 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.

Figure 6

The numerical simulation results at the upstream valve for the experimental facility (Jones et al. 2018) with various initial flow rates (Q0) and initial pressure heads (H0) indicated above each figure.

Figure 6

The numerical simulation results at the upstream valve for the experimental facility (Jones et al. 2018) with various initial flow rates (Q0) and initial pressure heads (H0) indicated above each figure.

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 Q0 = 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 H0 + 81.43) with a favorable coefficient of determination (R2 = 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 (H0 = 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.

Figure 7

Comparison of the volume of intrusion for a range of different transient conditions: (a) varying initial head and initial flow Q0 = 2 L/s; (b) varying initial flow rate and initial head H0 = 20 m. The dashed line represents the linear relation fitted to the experimental data reported by Jones et al. (2018).

Figure 7

Comparison of the volume of intrusion for a range of different transient conditions: (a) varying initial head and initial flow Q0 = 2 L/s; (b) varying initial flow rate and initial head H0 = 20 m. The dashed line represents the linear relation fitted to the experimental data reported by Jones et al. (2018).

Viscoelastic behavior of pipes

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.

Figure 8

The effect of viscoelastic behavior of the pipe wall on the volume of intrusion for a range of different transient conditions: (a) varying initial head and initial flow Q0 = 2 L/s; (b) varying initial flow rate and initial head H0 = 20 m. The continuous and dash-dot line represent the result of the viscoelastic and elastic models, respectively.

Figure 8

The effect of viscoelastic behavior of the pipe wall on the volume of intrusion for a range of different transient conditions: (a) varying initial head and initial flow Q0 = 2 L/s; (b) varying initial flow rate and initial head H0 = 20 m. The continuous and dash-dot line represent the result of the viscoelastic and elastic models, respectively.

CASE STUDIES AND DISCUSSION

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 (xl) 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 U0 = 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.

Table 2

Reservoir and pipeline specifications

Test (1)Test (2)Common specifications
D (mm) 600 108  (m) 
L (m) 2,300 540  (m) 
(m) 40 14  (m) 
xl 800; 1,100; 1,500 180; 270; 360 f 0.02 
 0.01; 0.001 0.01; 0.001  0.67 
 0.14; 0.28; 0.56 0.0045; 0.0091; 0.018 csteel (m/s) 1,000 
H0 (m) 41; 44; 57 15; 19; 35 cpvc (m/s) 390 
U0 (m/s) 0.5; 1; 2 0.5; 1; 2 Valve closure time (s) 
Test (1)Test (2)Common specifications
D (mm) 600 108  (m) 
L (m) 2,300 540  (m) 
(m) 40 14  (m) 
xl 800; 1,100; 1,500 180; 270; 360 f 0.02 
 0.01; 0.001 0.01; 0.001  0.67 
 0.14; 0.28; 0.56 0.0045; 0.0091; 0.018 csteel (m/s) 1,000 
H0 (m) 41; 44; 57 15; 19; 35 cpvc (m/s) 390 
U0 (m/s) 0.5; 1; 2 0.5; 1; 2 Valve closure time (s) 
Table 3

Total volume of contaminant parcel (VCPt) for various pipe systems and transient conditions

Test 1
Test 2
xL (m)U0 (m/s)VCPt, ELVCPt, VExL (m)U0 (m/s)VCPt, ELVCPt, VE
0.001 800 0.50 4.73 0.00 0.001 180 0.50 0.33 0.06 
1.04 44.42 0.00 1.00 0.63 0.21 
2.09 82.59 8.08 2.03 0.73 0.18 
0.001 1,100 0.50 6.72 0.00 0.001 270 0.51 0.46 0.10 
1.04 60.47 0.00 1.00 0.88 0.32 
2.08 113.68 10.99 2.03 1.09 0.27 
0.001 1,500 0.50 5.82 0.00 0.001 360 0.51 0.43 0.10 
1.04 64.02 0.00 1.00 0.88 0.40 
2.08 130.20 14.06 2.02 1.18 0.34 
0.01 800 0.50 0.00 0.00 0.01 180 0.51 0.27 0.00 
1.01 45.55 0.00 1.07 0.97 0.63 
2.08 175.79 0.00 2.12 1.85 0.92 
0.01 1,100 0.51 0.00 0.00 0.01 270 0.52 0.38 0.00 
1.01 59.64 0.00 1.04 1.16 0.92 
2.1 221.30 0.92 2.1 2.20 1.36 
0.01 1,500 0.5 0.00 0.00 0.01 360 0.5 0.17 0.00 
1.0 31.65 0.00 1.0 0.95 0.85 
2.0 186.13 0.00 2.1 1.83 1.48 
Test 1
Test 2
xL (m)U0 (m/s)VCPt, ELVCPt, VExL (m)U0 (m/s)VCPt, ELVCPt, VE
0.001 800 0.50 4.73 0.00 0.001 180 0.50 0.33 0.06 
1.04 44.42 0.00 1.00 0.63 0.21 
2.09 82.59 8.08 2.03 0.73 0.18 
0.001 1,100 0.50 6.72 0.00 0.001 270 0.51 0.46 0.10 
1.04 60.47 0.00 1.00 0.88 0.32 
2.08 113.68 10.99 2.03 1.09 0.27 
0.001 1,500 0.50 5.82 0.00 0.001 360 0.51 0.43 0.10 
1.04 64.02 0.00 1.00 0.88 0.40 
2.08 130.20 14.06 2.02 1.18 0.34 
0.01 800 0.50 0.00 0.00 0.01 180 0.51 0.27 0.00 
1.01 45.55 0.00 1.07 0.97 0.63 
2.08 175.79 0.00 2.12 1.85 0.92 
0.01 1,100 0.51 0.00 0.00 0.01 270 0.52 0.38 0.00 
1.01 59.64 0.00 1.04 1.16 0.92 
2.1 221.30 0.92 2.1 2.20 1.36 
0.01 1,500 0.5 0.00 0.00 0.01 360 0.5 0.17 0.00 
1.0 31.65 0.00 1.0 0.95 0.85 
2.0 186.13 0.00 2.1 1.83 1.48 

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 U0 = 0.27 ms–1, xL = 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 QL 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).

Figure 9

Time series of (a) pressure and (b) discharge at the leak location in test 2, elastic pipe, U0 = 0.27 m/s, xl = 270, and δ = 0.001.

Figure 9

Time series of (a) pressure and (b) discharge at the leak location in test 2, elastic pipe, U0 = 0.27 m/s, xl = 270, and δ = 0.001.

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.

Effect of pipe material

Since the wave speed in viscoelastic pipes is smaller than elastic pipes, less volume of intrusion is expected. Based on several scenarios implemented, an index to quantify the effect of pipe material on the amount of intrusion can be defined. For example, in the first test, among all 36 cases, half are run for elastic (EL) and half for viscoelastic (VE) pipe. The average of VCPt for each pipe material denoted by VCPtm can be defined by: 
formula
(18)
where Nexp denotes the number of tests whose average are computed for each material; Nexp = 18 herein. The following index, which is defined by the ratio of mean intrusion in viscoelastic pipes to that in elastic pipes, illustrates the reduced extent of intrusion in viscoelastic materials: 
formula
(19)

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.

Pressure–velocity interdependence

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 Vpipe 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 U0 = 2 m/s is less than that of U0 = 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).

Figure 10

Pressure–velocity interdependence for several tests with = 14 m, (a) elastic pipe (a = 1,000 ms–1), δ = 0.001; (b) elastic pipe (a = 1,000 ms–1), δ = 0.01; (c) viscoelastic pipe (a = 390 ms–1), δ = 0.001; (d) viscoelastic pipe (a = 390 ms–1), δ = 0.01.

Figure 10

Pressure–velocity interdependence for several tests with = 14 m, (a) elastic pipe (a = 1,000 ms–1), δ = 0.001; (b) elastic pipe (a = 1,000 ms–1), δ = 0.01; (c) viscoelastic pipe (a = 390 ms–1), δ = 0.001; (d) viscoelastic pipe (a = 390 ms–1), δ = 0.01.

Leak size and leak location

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.

Figure 11

Dimensionless time histories of flow rate (inflow/outflow) at the leak location in elastic pipe for U0=2 ms–1 and different δ, (a) test 1, xl= 1,150 m; (b) test 2, xl= 360 m.

Figure 11

Dimensionless time histories of flow rate (inflow/outflow) at the leak location in elastic pipe for U0=2 ms–1 and different δ, (a) test 1, xl= 1,150 m; (b) test 2, xl= 360 m.

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 (xL/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)).

Figure 12

Dimensionless leak discharge for different leak locations, elastic pipe (a = 1,000 ms–1), U0 = 1 ms–1, δ = 0.01, (a) test 1; (b) magnification of the indicated square in (a); (c) test 2; (d) magnification.

Figure 12

Dimensionless leak discharge for different leak locations, elastic pipe (a = 1,000 ms–1), U0 = 1 ms–1, δ = 0.01, (a) test 1; (b) magnification of the indicated square in (a); (c) test 2; (d) magnification.

Leak size–fluid velocity interdependence

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.

Figure 13

Leak size–fluid velocity interdependence for test 1, elastic pipes (a = 1,000 ms–1), Δh = 40 m, (a) xl = 800 m; (b) xl = 1,100 m; (c) xl = 1,500 m.

Figure 13

Leak size–fluid velocity interdependence for test 1, elastic pipes (a = 1,000 ms–1), Δh = 40 m, (a) xl = 800 m; (b) xl = 1,100 m; (c) xl = 1,500 m.

DISCUSSION

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 (xL), whereas the inflow extent is proportional to the effective leak size (Ae) and the Joukowsky pressure (Hres ± V0c/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 Hres and V0 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.

CONCLUSIONS

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).

ACKNOWLEDGEMENTS

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.

SUPPLEMENTARY MATERIAL

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

REFERENCES

REFERENCES
Al-Omari
A. S.
Chaudhry
M. H.
2001
Unsteady-state inverse chlorine modeling in pipe networks
.
J. Hydraul. Eng.
127
(
8
),
669
677
.
Bergant
A.
Simpson
A. R.
Tijsseling
A. S.
2006
Water hammer with column separation: a historical review
.
J. Fluids Struct.
22
(
2
),
135
171
.
Besner
M. C.
Ebacher
G.
Jung
B. S.
Karney
B.
Lavoie
J.
Payment
P.
Prévost
M.
2010
Negative pressures in full-scale distribution system: field investigation, modelling, estimation of intrusion volumes and risk for public health
.
Drinking Water Eng. Sci.
3
(
2
),
101
106
.
Brunone
B.
1999
Transient test-based technique for leak detection in outfall pipes
.
J. Water Resour. Plann. Manage.
125
(
5
),
302
306
.
Chaudhry
M. H.
2014
Applied Hydraulic Transients
.
Springer
,
New York
,
USA
.
Collins
R. P.
Boxall
J. B.
Karney
B. W.
Brunone
B.
Meniconi
S.
2012
How severe can transients be after a sudden depressurization?
J. Am. Water Works Assn.
104
(
4
),
E243
E251
.
Fernandes
C.
Karney
B.
2004
Modelling the advection equation under water hammer conditions
.
Urban Water J.
1
(
2
),
97
112
.
Fontanazza
C. M.
Notaro
V.
Puleo
V.
Nicolosi
P.
Freni
G.
2015
Contaminant intrusion through leaks in water distribution system: experimental analysis
.
Proc. Eng.
119
,
426
433
.
Fox
S.
Collins
R. P.
Boxall
J. B.
2017
Experimental study exploring the interaction of structural and leakage dynamics
.
J. Hydraul. Eng.
143
(
2
),
04016080
.
Gullick
R. W.
LeChevallier
M. W.
Case
J.
Wood
D. J.
Funk
J. E.
Friedman
M. J.
2005
Application of pressure monitoring and modelling to detect and minimize low pressure events in distribution systems
.
J. Water Supply Res. Technol. Aqua
54
(
2
),
65
81
.
Jaynes
D. B.
Rogowski
A. S.
1983
Applicability of Fick's law to gas diffusion 1
.
Soil Sci. Soc. Am. J.
47
(
3
),
425
430
.
Jones
S.
Shepherd
W.
Collins
R.
Boxall
J.
2018
Experimental quantification of intrusion volumes due to transients in drinking water distribution systems
.
J. Pipeline Syst. Eng. Pract.
10
(
1
),
04018026
.
Keramat
A.
Tijsseling
A. S.
2012
Waterhammer with column separation, fluid-structure interaction and unsteady friction in a viscoelastic pipe
. In: (
Anderson
S.
, ed.).
Proceedings of the 11th International Conference on Pressure Surges, Lisbon, Portugal
,
October 24–26
.
BHR Group Limited
,
Cranfield
,
UK
, pp.
443
460
.
Keramat
A.
Tijsseling
A. S.
Ahmadi
A.
2010
Investigation of transient cavitating flow in viscoelastic pipes
. In:
Proceedings of the 25th IAHR Symposium on Hydraulic Machinery and Systems, Timişoara, Romania
. Vol.
12
,
September 20–24
,
IOP Publishing
,
UK
.
Keramat
A.
Tijsseling
A. S.
Hou
Q.
Ahmadi
A.
2012
Fluid-structure interaction with pipe-wall viscoelasticity during water hammer
.
J. Fluids Struct.
28
,
434
455
.
Keramat
A.
Gaffarian Kolahi
A.
Ahmadi
A.
2013
Waterhammer modelling of viscoelastic pipes with a time-dependent Poisson's ratio
.
J. Fluids Struct.
43
,
164
178
.
Kirmeyer
G. J.
Friedman
M.
Martel
K.
Howie
D.
LeChevallier
M.
Abbaszadegan
M.
Karim
M.
Funk
J.
Harbour
J.
2001
Pathogen Intrusion Into the Distribution System
.
AWWA Research Foundation
,
Denver, CO
,
USA
.
Lindley
T. R.
2001
A Framework to Protect Water Distribution Systems Against Potential Intrusions
.
Doctoral dissertation
,
University of Cincinnati
,
Cincinnati, OH
,
USA
.
Mansour-Rezaei
S.
2013
Contaminant Intrusion in Water Distribution Systems: Advanced Modelling Approaches
.
PhD Thesis
,
University of British Columbia
,
Kelowna, BC
,
Canada
.
Mansour-Rezaei
S.
Naser
G.
2013
Contaminant intrusion in water distribution systems: an ingress model
.
J. Am. Water Works Assn.
105
(
1
),
E29
E39
.
Massey
B. S.
1998
Mechanics of Fluids
.
Taylor & Francis
,
Cheltenham, UK
.
Mays
L. M.
2000
Water Distribution Systems Handbook
.
McGraw Hill and American Water Works Association
,
New York
,
USA
.
Meniconi
S.
Brunone
B.
Ferrante
M.
Berni
A.
Massari
C.
2011
Experimental evidence of backflow phenomenon in a pressurised pipe
. In:
Proc. 11th Intl. Conference on Computing and Control for the Water Industry – CCWI
.
American Water Works Association
,
Washington, DC
,
USA
.
Meniconi
S.
Brunone
B.
Ferrante
M.
Capponi
C.
Carrettini
C. A.
Chiesa
C.
Lanfranchi
E. A.
2015
Anomaly pre-localization in distribution–transmission mains by pump trip: preliminary field tests in the Milan pipe system
.
J. Hydroinf.
17
(
3
),
377
389
.
Mora-Rodríguez
J.
López-Jiménez
P. A.
Ramos
H. M.
2012
Intrusion and leakage in drinking systems induced by pressure variation
.
J. Water Supply Res. Technol. Aqua
61
(
7
),
387
402
.
Naser
G.
Karney
B. W.
2008
A transient 2-D water quality model for pipeline systems
.
J. Hydraul. Res.
46
(
4
),
516
525
.
Payesteh
M.
Keramat
A.
2017
Sensitivity analysis of hydraulic parameters on contaminant intrusion in transient conditions
.
AUT J. Civil Eng.
In press
.
doi:10.22060/ceej.2018.14230.5595
.
Rossman
L. A.
Boulos
P. F.
1996
Numerical methods for modeling water quality in distribution systems: a comparison
.
J. Water Resour. Plan. Manage.
122
(
2
),
137
146
.
Schuster
C. J.
Ellis
A. G.
Robertson
W. J.
Charron
D. F.
Aramini
J. J.
Marshall
B. J.
Medeiros
D. T.
2005
Infectious disease outbreaks related to drinking water in Canada, 1974–2001
.
Can. J. Public Health/Rev. Can. Sante'e Publique
96
(
4
),
254
258
.
Soares
A. K.
Covas
D. I. C.
Ramos
H. M.
Rei
L. F. R.
2009
Unsteady flow with cavitation in viscoelastic pipes
.
Int. J. Fluid Mach. Syst.
2
(
4
),
269
277
.
Starczewska
D.
Collins
R.
Boxall
J.
2015
Occurrence of transients in water distribution networks
.
Proc. Eng.
119
,
1473
1482
.
Swamee
P. K.
Sharma
A. K.
2008
Design of Water Supply Pipe Networks
.
John Wiley & Sons, Inc.
,
Hoboken, NJ
,
USA
.

Supplementary data