Improving incomplete mixing modeling for junctions of water distribution networks

The mixing phenomenon is increasingly recognized as a serious concern in drinking water quality models of water distribution networks (WDNs). Indeed, in most existing water quality models (such as EPANET (Environmental Protection Agency Network Evaluation Tool) (Rossman 2000) or WaterCAD (Water Computer-Aided Design)/WaterGEMS (WATER Geographic Engineering Modeling System) (Bentley Systems 2020)), the mixing of water quality constituents at junctions is simplified and considered to be complete and instantaneous, which can greatly reduce the accuracy of the predicted concentration values.

A complete mixing assumption was adopted in drinking water quality models to reduce simulation time and model complexity. The degree of mixing at junctions depends on a wide range of parameters, such as the junction type (cross or tee junction), junction configuration (e.g., junctions with different pipe diameters in their legs), and hydraulic conditions. Due to the importance of water quality models in urban WDNs and technological advances provided by high-speed computers, more accurate models have been developed within the last decades, allowing for the simulation of mixing at junctions instead of considering the simplistic assumption of complete mixing.

Fowler & Jones (1991) and Ashgriz et al. (2001) were among the first researchers to find that mixing should not be considered complete at the junctions of WDNs, unlike the assumption employed in most existing water quality software such as EPANET (Rossman 2000) or WaterCAD (Water Computer-Aided Design)/WaterGEMS (WATER Geographic Engineering Modeling System) (Bentley Systems 2020). Many researchers conducted laboratory studies on such a mixing phenomenon to investigate the details of mixing in cross junctions (Orear et al. 2005; Ho et al. 2006; McKenna et al. 2007; McKenna et al. 2008; Song et al. 2009; van Summeren et al. 2017; Song et al. 2018; Shao et al. 2019; Yousefian & Duchesne 2022a) or double-tee junctions (Orear et al. 2005; Ho et al. 2006; Ung et al. 2014; Cotes et al. 2016; Ung 2016; Gilbert et al. 2017; Song et al. 2018; Shao et al. 2019; Grbčić et al. 2020b). Likewise, many authors used experimental results of the mixing phenomenon at junctions to calibrate different computational fluid dynamics (CFD) tools (Orear et al. 2005; Webb & van Bloemen Waanders 2006; Ho et al. 2007; Webb 2007; Romero-Gomez et al. 2008c; Braun et al. 2015; Cotes et al. 2016; Hernández Cervantes et al. 2018; Grbčić et al. 2020b; Grbčić et al. 2020a). Alternatively, other researchers focused on the impact of the mixing phenomenon on sensor network designs, optimal sensor placements, and the distribution of disinfectants or residual chlorine in WDNs (Romero-Gomez et al. 2008a; Romero-Gomez & Choi 2011; Mompremier et al. 2015; Mompremier et al. 2017).

Several attempts have been made to develop different models for incomplete mixing (IM). AZRED (AZ refers to the University of Arizona, where the model was developed, and RED is a Spanish word for network) is one of the first models to project the concentrations of solute at junction outlets (Austin et al. 2008; Choi et al. 2008). In this model, the ratio of Reynolds (Re) number in one inlet to the other inlet (inlet Re numbers ratio) and the ratio of Re number in one outlet to the other outlet (outlet Re numbers ratio) are independent variables used to compute the concentration in outlets.

EPANET-BAM (Bulk Advective Mixing) is another model for IM that originally considered only the advection term of the transport equation and neglected the effect of diffusion in mixing (Ho 2008). In this model, the solute concentration in the outlet adjacent to the inlet with a higher flow rate was considered equal to the solute concentration in the inlet with a higher flow rate. The solute concentration of the other outlet was then calculated based on the mass balance equation. Later on, a scaling parameter was added in EPANET-BAM to take into account the effect of the diffusion term in the mixing phenomenon (McKenna et al. 2007; Ho & Khalsa 2008; Romero-Gomez et al. 2008b).

Shao et al. (2014) developed an empirical model for IM, whereby the independent variables are the flow momentum of inlets and outlets, defined as the ratio of the squared flow rate to the pipe section area. The dependent variable of this model is a newly defined flow distribution factor employed to project the solute concentration in outlets.

Besides AZRED, EPANET-BAM, and the analytical model by Shao et al. (2014), others were developed specifically for cross junctions. Hernández Cervantes et al. (2021) proposed a new model for IM in cross junctions where the inlet flow ratio, the outlet flow ratio, and the inlet concentration ratio are used as input variables, while the ratio of concentrations in the outlets is the output. In addition, Hammoudi (2021) suggested a regression equation to estimate the mass fraction at the outlets of a cross junction based on the ratio of flows in the inlets, the ratio of flows in the outlets, the density and diffusivity of the contaminant, the pipe diameter of the cross junction, and the Re number of contaminated inlets for the mixing in WDNs.

Among the existing models, two models can be employed for junctions with two pipe sizes where the opposite legs have the same pipe size. The first of these, EPANET-BAM-WRAP, is an IM model based on EPANET that estimates the solute concentration in the outlets of the junctions with two pipe sizes (Ho & O'Rear 2009). In this model, like in EPANET-BAM, a scaling parameter is defined to consider the effect of diffusion in mixing, for each junction. The second is an analytical model by Yu et al. (2014b), where it is assumed that the main pipes (the inlet with the higher flow rate and the outlet opposing that inlet) have a larger pipe diameter. This model can be used for three ratios of pipe diameter (mm), namely 25/25, 25/32, and 25/50, within cross and double-tee junctions (Yu et al. 2014b, 2014a, 2016). In this empirical model, a new mixing index was proposed based on the conservation of flow and solute at junctions. Finally, a regression equation was employed in this model to project the mixing index as a function of the inlet pipe diameter ratio, the inlet Re number ratio, and the outlet Re number ratio.

Further details regarding the existing IM models and studies carried out on the mixing phenomenon in the junctions of WDNs can be found in Yousefian & Duchesne (2022b). This review paper revealed that most existing IM models were developed based on laboratory experiments with small pipe diameters (25–50 mm) and low pressure.

Yousefian & Duchesne (2022a) studied the impact of pressure on the mixing phenomenon and found that for pressures higher than 140 kPa, the pressure does not have a significant impact on the mixing phenomenon. Later, Yousefian (2023) investigated the accuracy of the existing IM models within laboratory experiments and showed that using the existing IM models for real-world pipe diameters (100 mm and more) can lead to erroneous results. Besides, the IM models that can work for junctions with two pipe sizes either need to be calibrated based on experimental results (e.g., EPANET-BAM-WRAP) or can only be used for three specific pipe diameter ratios (e.g., analytical model by Yu et al. (2014b)).

Furthermore, Yousefian & Duchesne (2022b) found that no model exists to consider IM for common junctions, such as those cross junctions with pipe diameters of 150 × 100 × 150 × 100 mm, 100 × 150 × 150 × 150 mm, or 100 × 100 × 50 × 100 mm. Therefore, the objective of this study was to develop an IM model that can be employed for cross junctions with real-world pipe diameters (e.g., 100 and 150 mm) and also for cross junctions with different pipe diameters in their legs.

Among the characteristics of drinking WDNs, pipe size, and pressure are important parameters that are different between the experiments performed to develop the existing IM models and real-world conditions. Yousefian & Duchesne (2022a) recently showed that pressure cannot significantly affect mixing. However, the same authors also showed, in Yousefian (2023), that the junction size has a significant impact on the mixing behavior, when the junction size is smaller than 100 mm. In order to determine the characteristics of the experimental setup in the laboratory, the pipe and flow characteristics of a North American WDN serving a population of about one million people were investigated. To do so, the number and type of junctions with IM, the size of commonly used pipes, the configuration of junctions, and the most frequent flow rates were investigated for that WDN. The mini WDN Laboratory of Institut National de la Recherche Scientifique (INRS; see description below) was then redesigned and rebuilt based on these conditions to study the mixing behavior under real-world conditions.

The most frequent conditions of real-world WDNs were obtained from the analysis of 24-h average water demand on a normal day in the selected North American city. In this analysis, all junctions have been examined to determine if the junction is within the complete mixing conditions or not in normal flow conditions. To do so, the junctions were first divided into cross junctions and double-tee junctions. Among each type, those junctions with two inlets and two outlets were then selected. Next, those junctions with adjacent inlets were identified as junctions with IM conditions. Moreover, for double-tee junctions, those with a distance between tee junctions less than 10 diameters of the pipes were considered double-tee junctions with IM. According to such criteria, approximately 65% of the cross junctions in this WDN were within the conditions of IM, while only about 1% of its double-tee junctions had characteristics that led to IM. Besides, about 60% of the cross junctions with IM characteristics have different pipe sizes in their legs. It was also found that the most frequently used pipe sizes in this WDN were 100 and 150 mm. Therefore, based on these observations, six cross junction configurations with pipe sizes of 50, 100, and 150 mm were selected to study the mixing behavior within them. Those configurations are listed in Table 1. For cross junctions with different pipe sizes, two different types of reducers were used: a reducer coupling that gradually changes the pipe diameter from 150 to 100 mm in a length of 10 cm (configurations B and C) and a reducer bushing in which the pipe diameter suddenly reduces from 100 to 50 mm (configurations E and F).

It was observed that the flow rates between 1.50 and 3.00 LPS were the most common during a normal day in the studied WDN. Accordingly, for laboratory experiments, two series of flow scenarios were tested: the first series was used to develop a new IM model for cross junctions with different pipe sizes (Table 2) and the second series was employed to validate the new IM model (Table 3).

In addition to the flow rate scenarios tested in the laboratory (which were selected through the analysis of a real-world WDN), others were tested through CFD simulations to fill the gap between the laboratory flow scenarios and some extreme flow rate scenarios observed in the real-world WDN. The flow rate scenarios selected for the CFD experiments are shown in Table 3.

Tap water was used in the whole network and the injection of a soluble tracer (salt; NaCl) into the network was carried out through a pulse (or an injection) pump. The distance between the salt injection point and the cross junction was about 3 m in all configurations to ensure that longitudinal mixing occurred entirely in water before reaching the cross junction (Figure 2(b)). In all laboratory experiments, the solution of salt was injected into the southern pipe (red arrow in Figure 2(a)), tap water was flowing in the western pipe (blue arrow in Figure 2(a)), and the northern and eastern pipes were the outlets (orange arrows in Figure 2(a)). In each cross junction configuration (Table 1) and for all flow scenarios of the laboratory experiments (Tables 2 and 3), the mixing phenomenon was investigated for two pressure values (320 and 430 kPa). Also, for each pressure and flow rate scenario within the cross junction configurations, 10 different concentrations of salt were injected into the southern pipe.

Numerical (CFD) experiments were also carried out to study the mixing phenomenon within the flow scenarios listed in Table 4, which complemented the laboratory experiment flow scenarios (listed in Tables 2 and 3). To do so, the k-epsilon model of OpenFOAM-v2210 within the solver of twoLiquidMixingFoam was used. OpenFOAM-v2210 is an open-source CFD tool using the finite volume method. The solver, twoLiquidMixingFoam, can consider the transport and mixing of two liquids with different specifications. OpenFOAM was selected for CFD simulations because it is an open-source tool. Since OpenFOAM works with the k-epsilon model, this closure model was employed for turbulent modeling. More details about the closure model and the governing equations of this model, as used in OpenFOAM, are presented in Goodarzi et al. (2020).

The CFD model was calibrated and validated based on the experimental results of all scenarios displayed in Table 2. The calibration was conducted by adjusting the turbulent Schmidt number value in order to minimize the difference between the dimensionless eastern outlet concentrations obtained from CFD simulations and laboratory observations, for all flow scenarios in Table 2 and all junction configurations in Table 1. Four turbulent Schmidt number values were tested: 0.01, 0.135 (as suggested by Romero-Gomez et al. (2008b)), 1.00, and a value adopted from the contour of the fitted turbulent Schmidt number provided by Romero-Gomez et al. (2008c). This last value is a function of both the inlet Re number ratio and the outlet Re number ratio.

The same fluid specification was considered for contaminated and clean water, namely a density of 1,000 kg/m³ and a kinematic viscosity of 1.00 × 10⁻⁶ m²/s. The geometry of the junctions was created through SolidWorks (Version 2021), and the snappyHexMesh tool of OpenFOAM was used to generate the mesh. Mesh analysis was carried out, and the final mesh was selected for each configuration from different meshes with a total number of cells in the range of 440,000 to 108,160,000. The boundary conditions applied for these simulations in all CFD experiments are presented in Table 5.

The laboratory flow scenarios that were used to develop the model (Table 2) were also used to calibrate the transport equation of the twoLiquidMixingFoam solver (the turbulent Schmidt number). Consistent with the research by Romero-Gomez et al. (2008b), it was observed in the current study that the CFD simulation results with a turbulent Schmidt number of 0.135 lead to the best agreement with the experimental results. Indeed, the root-mean-squared error (RMSE) between the simulated and observed dimensionless concentrations in the eastern pipe was 0.214, 0.018, 0.326, and 0.029, respectively, for turbulent Schmidt numbers 0.01, 0.135, and 1.00, and the turbulent Schmidt number values adopted from the contour provided by Romero-Gomez et al. (2008c). Therefore, a value of 0.135 was assigned to the turbulent Schmidt number for all CFD simulations carried out in this research.

In order to fit a single regression equation for all studied configurations, an analysis was carried out to find the principal combination of independent variables. The principal combination of independent variables was defined as the one based on which the outlet concentrations could be sorted. Different ratios of inlet and outlet velocities, ratios of inlets and outlet flow rates, and ratios of inlet and outlet Re numbers were taken into account to find the principal combination.

In some IM models, the inlet Re number ratio (⁠⁠) and the outlet Re number ratio (⁠⁠) were taken into account as independent variables (Austin et al. 2008; Yu et al. 2014b). However, it was observed from our experimental results that these ratios were not the principal combination explaining the mixing behavior (concentrations at the outlets). Indeed, in some cases, the outlet concentrations were different when the inlet Re number ratio and the outlet Re number ratio were the same, as shown in Table 6.

Table 6

Counter example for considering inlet and outlet Re number ratios as the principal combination

Configuration .	Flow scenario .	.	.	.	.
A	Lab-3	0.50	1.00	0.70	0.03
E	Lab-5	0.50	1.00	0.49	0.01
F	Lab-1	0.50	1.00	0.88	0.07

Configuration .	Flow scenario .	.	.	.	.
A	Lab-3	0.50	1.00	0.70	0.03
E	Lab-5	0.50	1.00	0.49	0.01
F	Lab-1	0.50	1.00	0.88	0.07

View Large

Four possibilities of sorting the eastern outlet concentrations as a function of the inlet and outlet flow ratios are illustrated in Figure 4. The other four possibilities for sorting the data are the exact mirrors of these figures. Accordingly, as it can be seen in Figure 4, it is not possible to sort the eastern concentrations in ascending or descending order based on the inlet flow ratio and the outlet flow ratio.

It is observed in Figure 5 that for all studied configurations, when the ratios of (⁠⁠) and (⁠⁠) are the same, the dimensionless concentrations in the eastern pipe (mixing behavior) are similar. Therefore, those two flow ratios of the eastern pipe to the southern pipe and of the eastern pipe to the western pipe were used as independent variables in the IM model developed in this study.

These last equations have some limitations which should be taken into account before employing them. First, the equations work only for cross junctions within the six studied configurations. The pipe diameters of cross junctions in laboratory experiments were 150 mm or less and, therefore, the proposed equations can be used for cross junctions within the pipe diameter range of 50 to 150 mm. The flow rates tested to develop this equation were between 0.50 and 80.00 LPS and the ratios (⁠⁠) and (⁠⁠) were between 0.50 and 40.00; it is recommended that the proposed equations are used for the same range of flow rates.

Figure 7 shows that the newly developed model leads to the least differences between the observed and computed concentrations. It is also observed that, for the studied flow scenarios, the model developed in this study does not have any outlier, while some are present in the results of the AZRED model and the IM model by Shao et al. (2014). The RMSE between the simulated and observed concentrations of the eastern pipe, for configurations A and D, is 0.02 for the newly developed model, while it is of 0.05 and 0.07, respectively, for the AZRED model and the IM model by Shao et al. (2014).

This newly developed equation can be integrated into EPANET to improve water quality modeling in real-world water distribution systems. The current EPANET model with complete mixing assumption leads, in some cases, to errors in applications of water quality simulations such as water age modeling, disinfectant residual modeling, and water quality sensor placement. Although there are some existing IM models that can improve the accuracy of water quality modeling, they were developed and applied only on a limited number of junction types. However, the newly developed equation could reduce this limitation by considering more types of junctions, such as cross junctions with real-world pipe diameters and cross junctions with different pipe diameters.

The mixing phenomenon was studied at cross junctions with equal and unequal pipe sizes. Laboratory experiments were conducted under real-world conditions, such as 100 and 150 mm pipe sizes. Besides, OpenFOAM was calibrated and employed to simulate mixing under flow scenarios complementing the laboratory-tested flow scenarios.

Most existing models with IM were developed within laboratory conditions (i.e., pipe sizes smaller than in real-world networks) and for junctions with equal pipe sizes. Some models developed for junctions with unequal pipe sizes either depend on predetermined parameters, like the scaling parameter of EPANET-BAM-WRAP, or work with only limited pipe size ratios, like 25/32 and 25/50 in the analytical model by Yu et al. (2014b). This study investigated real-world cross junction sizes, such as 150 × 100 × 150 × 100 mm and 100 × 150 × 150 × 150 mm, which no existing IM model can simulate. A single equation was proposed to compute outlets' concentrations for all considered junctions' configurations. Based on the obtained results, it can be concluded that

1. For the studied junction configurations and flow rate scenarios, a Schmidt number of 0.135 showed the best agreement between outlets' concentrations obtained from laboratory experiments and numerical simulations.

2. The outlet concentrations could not be sorted in ascending or descending order for any of the following ratios: inlet Re number ratio and outlet Re number ratio, inlet flow rate ratio and outlet flow rate ratio, and the momentum ratio (defined by Shao et al. (2014)).

3. The outlet concentrations could be sorted in ascending or descending order, as a function of the ratio of the flow rate of the eastern pipe to the flow rate of the southern pipe and of the ratio of the flow rate of the eastern pipe to the flow rate of the western pipe; those two ratios were thus used to compute the outlet concentrations in the developed model.

4. The outlet concentrations computed by the developed model are closer to observations than those computed by two other IM models, namely AZRED and the IM model by Shao et al. (2014).

Although the developed model in this study could work for a broader range of junction configurations, it is recommended to be used with flow rates between 0.50 and 80.00 LPS and flow rate ratios between 0.5 and 40. Besides, it is important to emphasize that the findings of this study are also subject to some limitations. First, the pipe sizes were limited to 150 mm or less. Second, cross junctions with unequal pipe sizes were installed using fittings for laboratory experiments. Third, the flow rates in laboratory experiments were limited to 1.50 and 3.00 LPS. Therefore, some issues need to be investigated in the future. Indeed, the principal combinations of independent variables identified in this research should be evaluated for other junctions' types and configurations (e.g., double-tee junctions, other cross junction sizes, etc.). Nevertheless, the developed model and the future proposed research of this study could improve water quality models in order to increase the use in disinfectant residual management, water quality simulations in WDNs with different water sources, water quality sensor placement, and contaminant source identification.

The authors would like to thank Mr Pierre-Olivier Schwarz and Mr Laurent Émond for their precious help in conducting this research. This research was funded by the Natural Sciences and Engineering Research Council of Canada, Grant No. RGPIN-2018-05473.

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

The authors declare there is no conflict.