Definition and application of the Péclet number threshold for water quality analysis in water distribution networks

Within the water distribution systems, velocity in pipes is significantly variable during the day. This causes laminar, transition and turbulent flow regimes to co-occur in the network which determines different advective and diffusive–dispersive contaminant transport mechanisms.

In which, is the solute concentration, D is the constant diffusion coefficient, u the mean flow velocity, k the coefficient of first order reaction, while t and x representing the variables of time and space, respectively. Although in many chemistry books, D is defined as diffusion coefficient (it determines how far molecules move on average in a given period of time; the diffusion coefficient depends on the size and shape of the molecule, and the temperature and viscosity of the solvent), in ADR equation it is mostly called ‘dispersion coefficient’. Obviously, in Equation (2), D is the ‘longitudinal dispersion coefficient’ or ‘axial dispersion coefficient’, since the dispersive process develops only longitudinally upon x.

Currently, to model the water quality within water networks, the EPANET (Rossman 2000) model is used which solves the contaminant transport equation by introducing simplifications based on solving the mass balance equation of the plug flow substance and on a simplified reaction kinetics. This is advantageous from a computational point of view as it does not produce particularly significant errors in purely turbulent flow regimes.

Other advective–reactive models have subsequently been developed to consider several species of contaminants acting simultaneously within the water distribution networks, to evaluate the WDS response to diverse contamination events (Eliades et al. 2016; Abhijith & Ostfeld 2022), neglecting the diffusive–dispersive phenomena.

Nevertheless, diffusive–dispersive processes assume great importance when low velocities and laminar/transition flow regimes take part and affect the actual propagation of contaminants within the system, as highlighted in Piazza et al. (2020). Furthermore, in Abokifa et al. (2020b) it is observed that neglecting dispersive transport and spatial aggregation demands can overestimate residual chlorine concentrations in dead ends and this generates more exacerbated errors in reduced demand scenarios. To overcome this error, in developing the WUDESIM model (open-source C/C ++ toolkit for modeling water quality), the authors determined three closed-form correction factors depending on the Péclet number and Damköhler number, to model the advective–dispersive transport and decay of the constituents (Abokifa et al. 2020a).

Subsequently, Romero-Gomez & Choi (2011) proposed a new formulation for the dispersion coefficient, as they realized 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 a condition differently applies to dispersion upstream and downstream from the contaminant injection.

For this reason, they have specified the dispersion coefficient considering the effect that the direction of mass flow has on the dispersion, backward and forward (Romero-Gomez & Choi 2011).

All the cited studies demonstrated that the use of diffusive–dispersive models produces results that are remarkably representative of the actual behaviour of the contaminant within the network, managing to be effective in the representation of the transport mechanism; on the other hand, extensive use of the diffusive–dispersive approach produces an increase in the needed computational effort that may be necessary for a small portion of the network only.

Waeytens et al. (2013) proposed an inverse modelling method to reconstruct the concentration field in a bounded domain, by solving the adjoint problem to quantify the sensitivity of the data misfit functional to the boundary controls. This allowed us to reduce the computational costs obtaining qualitatively and quantitatively similar results to those observed.

For this reason, using the appropriate numerical model according to the case to be analysed is of fundamental importance.

Several studies (Rapp 2017) state that, if the Péclet number is small, diffusive mass transport predominates over advective mass transport; otherwise, if the Péclet number is high, the advection is strongly dominant and diffusive processes can be considered negligible. This is generally true if we adopt a 1D cross-section averaged simplification is applied, but as pointed out in Chhabra & Shankar (2018), near the solid boundaries, diffusive processes are dominant thus giving rise to concentration boundary layers, which is essential to studying the interaction of both transport mechanisms.

In the literature, numerous definitions of the Péclet number (Huysmans & Dassargues 2005) are available that do not allow identifying a single threshold related to the dominancy between the two contaminant transport mechanisms.

This study aims to verify whether the validity of the Péclet number application to the 1D Advection–Dispersion equation is still valid when solving the Advection–Dispersion equation in the 2D case and, consequently, if such a numerical threshold can be used in real networks, where cross-sectional variability is relevant (Hart et al. 2021), to discriminate pipes in which 1D advective simplification can be used.

The determination of the Péclet number was carried out using three models: the advective EPANET model (Rossman 2000), the 1D advective–diffusive model (Shang et al. 2023) and the new dynamic model proposed by the authors EPANET-DD (Dynamic-Dispersion) model (Piazza et al. 2022), on a linear pipeline, contaminated using a Gaussian-type concentration pattern of sodium chloride (⁠ m²/s). Three scenarios have been considered in which the pipe length varies in a range of 0.018–107 m, the diameter in a range of 0.024–0.50 m and the flow rate is variable in a range between 0.0000397 and 0.56 m³/s randomly. These ranges have been chosen in order to take into account the pipes' variability that may be found in the distribution networks.

Similar to the NSE coefficient, KGE = 1 indicates perfect agreement between the simulations and observations. For KGE values ≤0, analogous to what the authors observed for NSE values, all negative values below the threshold KGE = 0 indicate results with poor model performance.

The root mean square error (RMSE) is the quadratic mean of the differences between the observed and predicted values. RMSE value equal to 0 would indicate a perfect fit to the data (Tyagi 2022).

The classical advective EPANET model (Rossman 1994) solves the advective transport equation using the finite volume method (FVM). The model, by assigning the mass of the substance to discrete volume elements of the pipeline network, solves a fundamental mass balance of the plug flow that considers advective transport and any kinetic reaction processes. In this case, the effect due to longitudinal dispersion is neglected, as it is not considered necessary in most operating conditions.

The 1D advective–diffusive model allows you to simulate the longitudinal diffusion of a solute in a water pipe. The approach used in the present study refers to that developed by Shang et al. (2023), in which the dispersion coefficient was calculated using a 2D stochastic approach called RWPT (Random Walk Particle Tracking) to model the motion of an instantaneous pulse of solute particles under both laminar and turbulent conditions.

The new model EPANET-DD solves the equations under quasi-steady flow conditions, solving the hydraulic problem under steady flow conditions with the EPANET-MATLAB-Toolkit (Eliades et al. 2016) and the advection–diffusion–dispersion equation under dynamic flow conditions in the 2D case with the classical random walk method (Delay et al. 2005), implementing the diffusion and dispersion equations proposed by (Romero-Gomez & Choi 2011).

Initially, the results from the simplified 1D advective approach and the 1D diffusive–dispersive approach were compared with the 2D diffusive–dispersive detailed approach to show the impact of diffusion processes and of the cross-sectional variability of contaminant distribution. A specific case, with a flow rate equal to 0.001291 m³/s, a pipe diameter equal to 0.036 m, and three sections with a length (denoted Sez) equal to 8.08, 55.18, and 80.37 m, was used as an example for a critical flow regime.

It is observed that in the case of a flow regime with well-established turbulence (Re > 50,000), the difference in the application of the 1D and 2D advective–dispersive model is minimal in correspondence with the first section (Figure 1(a)), in terms of concentration peak value, which is lower by applying the 2D model. By increasing the pipe length, a difference between the two models is noticed. The 2D advective–dispersive model tends to anticipate the contamination event compared to the one-dimensional case (Figure 1(b)). This deviation is more evident considering longer pipe lengths (Figure 1(c)). Furthermore, it is observed that the use of the 2D model results in a lower peak concentration than the 1D modelling.

A similar behaviour is observed considering the critical flow regime (3,000 < Re < 50,000), in which the event simulated using the 2D advective–dispersive model is anticipated and this difference increases as the pipe length increases (Figures 2 and 3). The anticipated event given by the 2D model is probably linked to the different velocity profiles. In fact, in the case of the 1D advective–dispersive model, the velocity profile is averaged along the cross-section of the tube. If the 2D model is used, the velocity profile changes to a logarithmic profile if turbulent shear prevails (Salama 2021).

Tables 1–3 show the geometric and hydraulic characteristics of the linear pipeline used in the study for the purely turbulent and critical flow regimes. Tables 4–6 show KGE, NSE, and RMSE coefficients determined by correlating the EPANET-DD (2D-ADR) and 1D ADR, EPANET-DD and EPANET (1D Advective), and 1D ADR and EPANET models.

Considering the objective of the paper, i.e. verifying if Péclet number can be used to discriminate the conditions (pipe lengths, diameters and flow regimes) for which 1D advective simplification can be accepted, we should expect higher correlation between EPANET-DD and EPANET and between 1D ADR and EPANET when the Péclet number is high. The correlation between EPANET-DD and 1D ADR shows the role of cross-sectional averaging simplification in the analysis.

It is observed that considering different pipeline lengths and the two different flow regimes, higher values of the axial Péclet number are related to lower correlations both between ADR approaches and advective simplification and between 2D approach and 1D simplification. Furthermore, this is specially true for the transitional flow regime (Table 2 and Figure 2) and for low values of Re (Table 3 and Figure 3).

Concerning Figure 2(c), it is noted that although the Reynolds number is high (Re = 45,652) the effect linked to dispersion-diffusion is evident and the simple advective model cannot describe the real transport mechanisms.

It is observed that for low values of the axial Péclet number (Pe < 5,000), the KGE coefficient has a notable variability (between 0 and 1). This is mainly due to two aspects related to the length of the pipeline and the flow regime that develops inside it. In fact, since the Péclet number is directly proportional to the length of the pipe, in the presence of pipes with limited lengths and purely turbulent flow regimes the KGE coefficient presents quite high efficiencies. If the flow regime inside the pipe is critical or laminar, the KGE coefficient values approach zero or are negative.

In Figure 4, we can see the non-applicability of the Péclet number to both the 2D and 1D advection–dispersion equation. In fact, the 1D advection–dispersion equation (Equation (5)) states that very high values of the Péclet number mean that the contribution linked to diffusion–dispersion can be neglected and that the process is dominated by advection. However, what has been stated is not confirmed by the results obtained, as despite high values of the Péclet number, there is no dominance of the advective processes over the diffusive–dispersive ones.

The present study allowed us to evaluate the applicability of the Péclet number to the 2D advective–dispersive–diffusive model EPANET-DD, in order to define a priori the dominance of one process over the other.

The research was developed using three models: the EPANET advective model, the 1D advective–dispersive model, and the 2D advective–dispersive–diffusive EPANET-DD model to simulate the motion of a conservative tracer inside a linear pipe.

The KGE, NSE and RMSE efficiencies of the three models used in the study were evaluated as a function of the axial and radial Péclet number.

The study highlighted that the Péclet number cannot be used as a parameter to identify the dominance of the advective or dispersive process in real systems, as it was observed that for low values in the axial Péclet number (Pe < 5,000), it has a notable variability of the KGE statistical parameter, varying between 0 and 1 when comparing the 1D and 2D advective–dispersive models, the 1D advective–dispersive model and advective model and the 2D advective–dispersive model and advective model. This is mainly due to the relevant role of pipe length when diffusive–dispersive processes are significant (such as in transitional and laminar flows).

A similar result was obtained by considering the relationship of the KGE parameter as a function of the radial Péclet number.

By comparing the models used in the present study to model a contamination event, it was observed that diffusive–dispersive processes take on considerable relevance, not only in the presence of laminar or critical flows, but also for turbulent flows. In fact, the results showed that in the presence of a purely turbulent flow regime, in which the literature states diffusive–dispersive processes can be neglected, in reality they are present and this is also confirmed by the low values of the axial Péclet number determined in this study.

This is linked to the definition of the Péclet number which is directly proportional to the pipe length. Therefore, for pipelines with short lengths, despite the flow regime occurring on the pipeline, the diffusive–dispersive processes are still relevant and cannot be neglected.

This study was carried out within the RETURN Extended Partnership and received funding from the European Union Next-GenerationEU (National Recovery and Resilience Plan – NRRP, Mission 4, Component 2, Investment 1.3 – D.D. 1243 2/8/2022, PE0000005).

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.