The accurate modeling and simulation of the spread of contaminants within water distribution networks (WDNs) is an important task for drinking water security. In commonly used water quality simulation platforms the mixing of concentrations at junctions is assumed to be a complete mixing behavior. Experimental investigations have shown that this assumption is only true for certain flow conditions. Indeed, mixing depends on the geometry and the load configuration and often it is incomplete or there is some shortcut and preferential flow pathway with no mixing at all. In this paper, we present a more representative model of mixing at nodes that is valid for double T- and cross-junctions. Statistics on real WDNs are used to define realistic junction scenarios. From the latter, two-dimensional and three-dimensional computational fluid dynamics (CFD) simulations for the mixing process at different types of junctions were made and compared to experiments. Both the simulations and experiments were in agreement and show a difference in mixing of 10% and more compared to perfect mixing models. CFD results were used to build a lookup table and Kriging interpolation was applied for entries not in the table.

## INTRODUCTION

The security of drinking water is increasingly recognized as a major challenge for municipalities and water utilities. The safety and/or security of drinking water can be threatened by natural disasters, accidents or malevolent attacks. In the event of contamination, water spreads rapidly and hence extensively before the problem is detected. Contaminated drinking water can induce major epidemics, disrupt economic life and create mass panic. To protecting this critical infrastructure, water utilities install sensors to monitor the quality of the water delivered. To make the right decisions, we need to know where the source of the contamination is located, what the extent of the pollution is and where it goes. Without this crucial information we would not be capable of undertaking the appropriate countermeasures. Also, what to do with the continuous and large stream of data? Experiments have been conducted on a small scale but not at larger scale (see, e.g., Whittle *et al.* 2010). There is an urgent need for a decision-making tool for security issues regarding water distribution networks (WDNs). The SMaRT-Online^{WDN} project (SMaRT-Online^{WDN} 2015) intends to contribute to this need. Its main objective is the development of an online security management and reliability toolkit for WDNs.

The SMaRT-Online^{WDN} toolkit is based on a multi-probe sensor network and simulation models that we aim to be as reliable and accurate as possible. The software solution relies on data treatment and assimilation from a sensor network of water quantity values (pressure, flow rate) and water quality values (e.g., free chlorine, pH, conductivity, turbidity, temperature). The core of the online security management toolkit consists of a grid of smart sensors in combination with an online simulation model. The boundary conditions of the network model are regularly updated by measurement data guaranteeing the compliance of the model with the observations. The demand is also calibrated online. With this information, the online security management toolkit is able to reflect the current hydraulic state of the entire system. In addition, monitoring of water quality parameters supports the detection of biochemical contamination of the drinking water. For identification of the source and prediction of the pollution extent currently and for the following hours, the transport simulation model is central.

Existing transport models in software package (e.g., EPANET 2015) possess some limitations regarding solute propagation in the network. In the context of hazardous waterborne products, it is important to improve the models in order to take appropriate countermeasures. In many cases the three-dimensional (3D) geometry inside the network has an impact on the spread of substances and the 1D (one-dimensional) transport simplified model may not be accurate enough (e.g., hydrodynamic dispersion in pipes for laminar flow and mixing at junctions). Recently, it was demonstrated that a more realistic model of the axial dispersion for laminar and transitional flow should be introduced (Romero-Gomez *et al.* 2008a; Romero-Gomez 2010; Romero-Gomez & Choi 2011). Likewise, experimental investigations have shown that for certain flow conditions at cross-junctions and double-T-junctions ideal perfect mixing is not fulfilled (Choi *et al.* 2008). Two models that complete EPANET were proposed: the AZRED model by Choi *et al.* (2008) and the Bulk Advective Mixing (BAM) model by Ho (2009). In practice, both models have certain shortcomings. The BAM model contains a parameter which has to be determined by several experiments for each individual junction. This is not applicable for real-world networks, which contain hundreds or thousands of junctions. The AZRED model on the other hand is strictly applicable for turbulent flow regime and extrapolates the experimental results for values in that range.

This paper is an extension of a conference paper presented at the 12th International Conference on ‘Computing and Control for the Water Industry – CCWI2013’ by Braun *et al.* (2014) and presents the work on the project SMaRT-Online^{WDN} concerning the development of a 1D mixing model for mixing of transient contamination pulses in WDN, based on two-dimensional (2D) and 3D computational fluid dynamics (CFD) simulations. The paper is structured as follows. First, the recent non-perfect mixing models (AZRED and BAM) are detailed. Next, statistics on the junctions of the WDN in Strasbourg are presented. Then, the CFD-based methods are introduced for the derivation of enhanced 1D mixing models, simulation results of mixing at T- and N-Junctions are then presented, and finally, a conclusion and an overview of future work are provided.

### Mixing models

This section presents the three most popular mixing models, namely the Perfect Mixing Model, the BAM model and AZRED.

#### Perfect Mixing Model

The concentrations ‘*c*’ are in mg L^{−1}, the flows rates ‘*Q*’ are in m^{3}s^{−1} and the sum of flow going in is equal to the sum of flow going out. This model is the one that is generally used (such as in EPANET) because it is the simplest and is easy to implement.

#### BAM model

This model completes the preceding model by modeling imperfect mixing at crosses, also called X-junctions. Indeed, Ho (2009) presents results of simulation and experiments that show that, in some cases, the mixing is not perfect. They mostly focus their work on cross-junctions with two adjacent entries, two adjacent outputs and equal diameters.

#### AZRED

The third model has been developed at the University of Arizona by Choi *et al.* (2008). This model relies on a lookup table filed with experimental results. It is valid for cross and double-T junctions composed of pipes of equal diameters. It was concluded that the degree of mixing depends mostly on the velocity (or Reynolds number) repartition in and out as well as the type of junction studied. Also, as this was an experiment-only study, results of mixing for extreme values of flow rates can only be extrapolated.

### Exemplary statistics on junctions of a WDN

To provide a good frame of reference for the flow conditions and junction configurations in a real network, exemplary statistics have been performed on the WDN of the Urban Community of Strasbourg (CUS). The results are based on a 24-hour simulation in Porteau (2015) using the reduced network graph with 12,730 links and 10,460 nodes; the pipe length is about 964 km. Graph simplification reduces the number of links and nodes by approximately 60% without affecting the hydraulics in the system and the type and distribution of junctions. This is done in five steps, which are as follows: (1) remove the dead end pipes or antennae without consumption; (2) merge adjacent pipes in series without consumption; (3) check for connected components without hydraulic source isolated by closed valves and remove them (manually operated); (4) Remove antenna with consumption if diameter is less than 95 mm and length is less than 200 m; and (5) merge adjacent pipes in series with small consumption at internal nodes. The simplification tools are available with the latest version of the Porteau software solution. The values of interest for the CFD model are the dominant junction types and the relevant flow conditions. The following conclusions can be drawn from the statistics: (1) the pipe diameter of *D*=100 mm is the dominant diameter in the network; (2) junctions with equal pipe diameters account for 44% of all T-junctions; (3) the most common dimension less distance for double T-junctions is *L*=5*D*; and (4) laminar flow regime is relevant for wide areas of the network. More detailed information on the statistics is given in Braun *et al.* (2014).

## CFD-BASED METHOD FOR THE DERIVATION OF ENHANCED 1D MIXING MODELS

Our aim is to develop enhanced mixing models which reflect the mixing of flows with different concentrations in a more realistic way than the existing approaches described above. These enhanced mixing models will be based on results gained from 2D and 3D CFD simulations. The aim is the implementation of the resulting solution in a 1D transport model. This approach is described in the following sections which cover in order: the governing equations for CFD models for mixing at T- and N-junctions, the model parameters and the method for the evaluation and extrapolation of the simulation results.

### CFD model for laminar mixing at T- and N-junctions

In this and the following section, the governing equations for CFD modeling for laminar and turbulent mixing will be presented, respectively. These models will be used for the investigation of mixing at T- and N-junctions.

*et al.*2002). The following simplifications can be applied: the fluid is incompressible, body forces and also inertia terms (due to the small scale and the slow variation of the flow) can be neglected. Navier-Stokes equations can be stated by the conservation Equation (3) and the momentum Equation (4). where

*ρ*denotes the density (constant parameter),

*p*the pressure and

**the velocity field. The spatio-temporal distribution of a single phase conservative contaminant (i.e. no biological degradation, no chemical reactions considered) is described by the convection-diffusion partial differential Equation (5), where**

*u**c*is the concentration of the contamination, and

*D*is the diffusion coefficient (which is the diffusivity of the contaminant). The velocity field

**for the convective transport is determined by Navier-Stokes Equations (3) and (4). The right hand side of (5) is equal to zero as we do not consider any sources inside the domain. The inlet of a contamination source is considered by a time-dependent boundary condition.**

*u*Figure 2(a) explains the nomenclature for T-junctions. The main inlet pipe is referred to by , the main outlet pipe is referred to by and the secondary pipe is referred to by in Figure 2(a). Figure 2(b) defines the nomenclature for N-junctions. This nomenclature can also be applied for the U-junction with equilateral secondary pipes. In an N-junction, the main inlet and the main outlet are defined the same way as at the T-junction. In addition, the secondary inlet pipe is referred to by and the secondary outlet pipe is referred to by . The pipe diameter is defined as which corresponds to typical value of the water utility CUS (Strasbourg).

#### Boundary conditions for T-junction

Flow rates at the inlet pipe and at the main outlet pipe are given, and a free pressure outlet boundary condition is set at the secondary outlet pipe . For the main inlet pipe a contamination impulse is set as concentration boundary conditions. The injection of the contaminant into the pipe is modeled as a Gaussian distribution focused at the pipe center. The temporal profile of the injection pulse is also chosen as a Gaussian distribution.

#### Boundary conditions for N-junction

The boundary conditions for the N-junction are chosen similarly. The only difference is that the flow rate secondary inlet pipe is also a boundary condition with zero concentration. The laminar entrance length, that gives the distance after which the laminar flow profile is fully developed, is set as (diameter).

### CFD model for turbulent mixing at T- and N-junctions

This section describes the steps taken for considering turbulence in the CFD model. Turbulent flow in general is described by the Navier-Stokes equations. This can be done in direct numerical simulations (DNS). DNS resolves the flow for eddies on all length scales encountered on the domain. Since the characteristic length scale gets smaller for higher Reynolds numbers we would have to use very fine meshes for the simulation. Because of the high computational cost in practical cases, we use turbulence modeling. Two popular methods are the Reynolds-Averaged Navier-Stokes equations (RANS) and the Large Eddy Simulation (LES).

*p*in the cross-section of a pipe can be represented by the sum of a time-averaged component and and an oscillating component and with an average value of zero. Averaging the fluctuations relates them to the mean flow. If we include these definitions into the Navier-Stokes equations we can write the momentum equation as follows:

The flow is calculated with respect to the time-averaged quantities. Since we cannot calculate the fluctuating properties they have to be modeled through the factor . This factor is also called the turbulent viscosity. For the RANS approach a number of models have been developed based on experimental and numerical observations. In cases that have to model the viscous sub-layer near surfaces as well as the outer layer like the pipe junctions presented in our studies the shear stress transport (SST) k-ω model is an appropriate choice (Menter 1994).

This diffusion coefficient is defined in a phenomenological sense, by analogy with the molecular diffusivity. In the gradient approach to modeling turbulent diffusion, is added to the molecular diffusion coefficient (Romero-Gomez *et al.* 2008b).

In contrast to the RANS model that averages in the time-domain, the LES model filters the spatial domain. This means that in LES we resolve large eddies and use a Sub-Grid-Scale (SGS) model to determine the effects of smaller vortices. With the Smagorinsky-Lilly Model, one of the first SGS models in LES was introduced in 1963 (Smagorinsky 1963). It is based on the equilibrium of the energy produced by the large eddies and the energy dissipated by the smaller vortices. A detailed description of the LES model for the interested reader is given by Ferziger & Perić (2002).

The boundary conditions for the turbulent flow simulations are chosen in the same way as for the laminar N-junction. For the solute transport, results presented in this paper we used stationary boundary conditions for the concentration and constant concentration over the cross-section of the inlet. For more details see Ansys (2011).

### Model parameters

The parameter values of the physical flow and mixing model are given by the density and the dynamic viscosity of water (Table 1). For the dilution model the diffusivity is the characterizing parameter. In the simulations presented in this paper the diffusivity of sodium chloride in water was chosen in order to have a better comparability with the experimental results in the project. The temperature for both the flow model and the dilution model is assumed to be room temperature.

Density | |

Dynamic viscosity | |

Diffusivity | |

Temperature |

Density | |

Dynamic viscosity | |

Diffusivity | |

Temperature |

In the case of the T-junction, only the relative outlet flow rate has to be defined. The indication of the pipes follows the one shown in Figure 2(a). The value of both the relative inflow and the relative outflow velocity are chosen to be . A third factor is introduced for N-junctions. It is defined by the dimensionless distance *L* between two T-junctions . Herestands for the distance between the two T-junctions. It is set in relation to the pipe diameter. The area of interest for the non-dimensional distance is defined as.

## Method for the evaluation and interpolation

In this section, the method of evaluation is described with a focus on the comparability of the transient CFD results and stationary 1D mixing models like BAM or AZRED. Further, interpolation methods are evaluated with respect to the use of the CFD simulation data in an online context.

### Method for the evaluation of the transient mixing behavior

### Interpolation approaches

The results created by the simulations described in the previous sections depend on a number of variables like the Reynolds number, the relative inflow and outflow velocities and the distance between two T-junctions. For the derived 1D mixing model we need to interpolate on an irregular grid. The accuracy of the method we use will determine the accuracy of the final model. There exist a number of interpolation approaches of which the most prominent range from simple linear interpolation to high dimensional spline interpolation. In the course of our research we took a close look at two specific methods for multivariate interpolation. They are the radial basis functions (RBF) interpolation scheme and the Gaussian Process Regression (GPR). Both methods can cope with irregular grids and are not as dependent on the location of the sample as a linear or spline interpolation. For the final implementation of the algorithm, these two methods have the advantage that they are very computationally efficient compared to the more classical approaches. This applies especially for calculating the value of a high dimensional point .

The RBF are usually applied to approximate functions that are only known at a finite number of points. It makes it possible to evaluate the approximated function often and efficiently.

Radial basis functions . | Φ(r)
. | Parameters . |
---|---|---|

Gaussian | c > 0 | |

Multiquadrics | c > 0 | |

Inverse multiquadrics | c > 0 | |

Inverse quadratics | c > 0 | |

Polyharmonic splines | r^{2k}log(r) | |

r^{2k−1} |

Radial basis functions . | Φ(r)
. | Parameters . |
---|---|---|

Gaussian | c > 0 | |

Multiquadrics | c > 0 | |

Inverse multiquadrics | c > 0 | |

Inverse quadratics | c > 0 | |

Polyharmonic splines | r^{2k}log(r) | |

r^{2k−1} |

*Kriging*or GPR (Cressie 1993) is introduced for multi-dimensional inference. The method is a stochastic model based on the definition of a probability distribution at each of the

*N*sample points . A spatial inference can be calculated as a linear combination of those for any point of the domain. The weights are computed by assuming stationary state for first and second moment. This results in a constant mean for all variables and a correlation between two variables that only depends on the distance between them. The usual chosen form of correlation is the following: with

*p*between 0 and 2. When

*p*is chosen as 2, the correlation is named

*Gaussian*and the mean-square of the distribution is infinitely differentiable. In a more general case,

*θ*and

*p*can be fitted using maximum likelihood estimation (MLE) (Biles

*et al.*2007). If the prior information is well chosen, this method is expected to be the best unbiased linear predictor (Cressie 1993).

To calculate the weights in Equation (16) an optimization problem is formulated that minimizes the variance taking into account the fact that the first moment is stationary. This method can also be used in the engineering field as a meta-modeling tool (as is done in our case) for building a black box model over a designed set of computer experiments.

## SIMULATION RESULTS OF MIXING AT T- AND N-JUNCTIONS

The aim of this section is to present, in order, results from 2D laminar simulations, the approach and results of the calibration of the turbulent Schmidt number for the RANS mixing models and a comparison between 2D and 3D simulations. In the final subsection the results of 3D simulations and experiments are briefly discussed.

The 2D simulations have been performed on a COMSOL platform (version 4.3) while the 3D model is implemented in ANSYS (version 14.0). By means of the 2D simulations, it should be evaluated if the results are comparable to the 3D results below. For the 3D model, we performed laminar simulations and turbulent LESs. In the simulations discussed below, involving 2D and 3D laminar and turbulent RANS simulations, we evaluate the accuracy of the reduction in dimension and the degree of abstraction in the turbulence model.

### Results from 2D simulations

The 2D laminar simulations performed on a COMSOL platform concentrate on three way T-junctions and different configurations of three-way N-junctions. A major conclusion for the junction in a T modification is, that it can be seen as an X-Junction under the condition that the flow at one of the pipes is equal to zero. The simulations have also shown that the complete mixing assumption can only be applied to N-junctions if the relative pipe length between the two T intersections is larger than . In all other cases including the X-junction with the two consequent junctions have to be treated as a single unit.

Presented in this section are the X-junction and the N-junction with a distance of . The results are shown in Figure 3. The three axes of the diagrams show the relative inflow velocity , the relative outflow velocity and the difference in mass fraction for the complete mixing model and the CFD leaving the main pipe of the junction. For the extreme values of and both models are in good agreement since the difference for both approaches is equal to zero. For the area of , Figure 3(a) shows the plane produced by the complete mixing model. This shows that, in contrast to the complete mixing assumption, the CFD simulation shows a negligible amount of mixing. For this effect can be seen as well but not as prominently. Figure 3(b) shows the difference between the complete mixing model and the CFD simulations for . It shows that the complete mixing model still overestimates the amount of mixing but that both models are closer to each other. This reinforces the fact that for long distances between the two junctions the complete mixing model is valid. Figure 4 compares the BAM model and the CFD simulation by the difference in mass fraction . It is obvious that both models are in relatively good agreement despite the area of . This can be explained by the enhanced mixing due to the diffusion neglected by the BAM model.

### Calibrating the turbulent Schmidt number for the RANS mixing models

When modeling the mixing process with RANS simulations we have to calibrate the turbulent Schmidt number to predict the appropriate amount of mixing. In the simulations presented in this paper, the benchmark for the calibration is the full 3D LES. This calibration has been done for a standard N-junction with the dimensionless distances . Further the Reynolds number has been considered for the values of and . Table 3 shows the mass fraction at outlet pipe (3) for the N-junction with and and for different values of the turbulent Schmidt number. For the cases described in this passage, it has been found that based on the LES gives reasonable results over a broad range of boundary conditions.

. | Sc_{T} = 0.2
. | Sc_{T} = 0.3
. | Sc_{T} = 0.4
. | Sc_{T} = 0.71
. | LES . |
---|---|---|---|---|---|

Re = 5,000 | 0.53 | 0.57 | 0.55 | 0.58 | 0.53 |

Re = 10,000 | 0.63 | 0.60 | 0.53 | 0.53 | 0.63 |

. | Sc_{T} = 0.2
. | Sc_{T} = 0.3
. | Sc_{T} = 0.4
. | Sc_{T} = 0.71
. | LES . |
---|---|---|---|---|---|

Re = 5,000 | 0.53 | 0.57 | 0.55 | 0.58 | 0.53 |

Re = 10,000 | 0.63 | 0.60 | 0.53 | 0.53 | 0.63 |

### Comparison between 2D and 3D simulations

After our initial simulation on the 2D model, the question poses itself how well the results of the simplified 2D model correspond to a more complex model that pays respect to all three spatial dimensions. To investigate this, we compared a 2D and a 3D model of an N-junction with a dimensionless distance between the junctions of . This is done for laminar as well as turbulent flow regime (RANS). The relative inflow and outflow are chosen as and. Table 4 shows the results of these simulations in term of the mass fraction . The results show that for the laminar flow regime, there is a considerable difference between the 2D and the 3D model. For the turbulent model no such difference can be observed. Looking at the concentration distribution on the domain we get a clue as to why there is such a large difference for the laminar model. Figure 5(a) shows that the flow at the first junction on the inlet side is forced into a stratified flow due to the inability to spread in any other direction. This constraint is directly imposed through the reduction in dimension. The cross-section of the symmetry plane of the 3D model shown in Figure 5(b) in contrast, reveals an area of low concentration at the upper wall of pipe . Unlike the 2D model, the flow from pipe 2 moves along the bottom of the pipe connecting the two junctions, resulting in a much lower mass fraction *w*_{3}. The observation that the mixing in laminar flow cannot be simplified as easily as in turbulent flow is supported by the experimental results presented by Choi *et al.* (2008) and seems to be the main reason for their model to be restricted to turbulent flows.

Mass fraction . | 2D . | 3D . |
---|---|---|

Laminar | 0.75 | 0.63 |

Turbulent | 0.7 | 0.7 |

Mass fraction . | 2D . | 3D . |
---|---|---|

Laminar | 0.75 | 0.63 |

Turbulent | 0.7 | 0.7 |

As stated before, in contrast to mixing in the laminar case, the effect of modeling in 3D in the studied case has no significant influence for turbulent mixing. In Figure 6(a), we can see that the 2D model develops a stratified flow similar to the laminar case, but we can see as well that the difference to the 3D model in Figure 6(b) is not nearly as severe as it is for the laminar model. An explanation for this behavior could be given by the fact that the enhanced mixing by the turbulent motion dominates the process.

Following the experience from the comparison of the 2D and 3D RANS simulations as well as the calibration of the turbulent Schmidt number, Table 5 shows the mass fraction for two different Reynolds numbers for turbulent flow regime and for a relative distance between the junctions of , , and . They confirm the observation that for a well-calibrated turbulent Schmidt number a simplified 2D RANS model can replicate the results of the full 3D LES.

w_{3,RANS}|w_{3,LES}
. | L = 5D
. | L = 8D
. | L = 10D
. | L = 20D
. |
---|---|---|---|---|

Re = 5,000 | 0.64|0.53 | 0.57|0.57 | 0.55|0.55 | 0.51|0.58 |

Re = 10,000 | 0.63|0.63 | 0.60|0.60 | 0.53|0.53 | 0.53|0.53 |

w_{3,RANS}|w_{3,LES}
. | L = 5D
. | L = 8D
. | L = 10D
. | L = 20D
. |
---|---|---|---|---|

Re = 5,000 | 0.64|0.53 | 0.57|0.57 | 0.55|0.55 | 0.51|0.58 |

Re = 10,000 | 0.63|0.63 | 0.60|0.60 | 0.53|0.53 | 0.53|0.53 |

### Comparison between 3D simulations and experiments

In Table 6, a comparison between some simulations and experiments is presented for the case Re = 5,000 and distance *L* = 5D. The difference between the experiments is less than 5% of total mass of contaminant and therefore they are in good agreement. The differences may be either due to difficulties in the measurements or lack of details on the geometry simulated. In the last column is the result for the case of perfect mixing, it can be observed there is a significant difference of about 10% and more between simulations and experiments.

Q_{in,rel}
. | Q_{out,rel}
. | Fluent LES simulations (%) . | Experiments (%) . | Perfect mixing (%) . |
---|---|---|---|---|

0.7 | 0.3 | 44–66 | 39.9–60.1 | 30–70 |

0.7 | 0.5 | 66–44 | 64.2–35.8 | 50–50 |

0.7 | 0.7 | 82–18 | 84.3–15.7 | 70–30 |

Q_{in,rel}
. | Q_{out,rel}
. | Fluent LES simulations (%) . | Experiments (%) . | Perfect mixing (%) . |
---|---|---|---|---|

0.7 | 0.3 | 44–66 | 39.9–60.1 | 30–70 |

0.7 | 0.5 | 66–44 | 64.2–35.8 | 50–50 |

0.7 | 0.7 | 82–18 | 84.3–15.7 | 70–30 |

## CONCLUSION AND FUTURE WORK

In this paper, we have presented investigations regarding the mixing behavior of three- and four-way junctions (T-, X- and N-junctions). CFD simulations have been performed in laminar, transitional and turbulent flow regimes. 2D and 3D domains have been used for the simulations. From the CFD simulation results, an enhanced 1D mixing model has been derived by a special interpolation method (interpolation on an irregular grid). The new 1D-mixing model has been integrated in the WDN simulation platform Porteau. The junction configurations and hydraulic conditions have been derived from a real-world, large network (water utility CUS, Strasbourg). A plausible contamination scenario has been defined (pulse-like inflow of contamination).

The performance of our enhanced mixing model has been compared with BAM and AZRED model. We have shown that for the X-junctions, the CFD-based approach is comparable to the BAM model and AZRED. Owing to consideration of diffusion in the CFD-based approach, slight improvements compared to BAM models have been achieved. Regarding N-junctions (which are not covered by BAM), the CFD-based approach in the laminar regime is more accurate than AZRED. As with AZRED, these values are extrapolated from turbulent regime. The T-junctions are not considered in the AZRED model, hence a comparison with AZRED is not possible.

A general advantage of our approach compared to AZRED is that even extreme values of the relative inflow and outflow are considered accurately. However, in the AZRED for extreme inflow and outflow conditions the results are extrapolated which may lead to inaccurate results.

The results of the 2D and 3D simulations have been compared. It has been shown that laminar mixing cannot be fully apprehended by 2D simulations. However, for turbulent flow this approach is feasible. It has been shown that with the Reynolds-averaged model, similar results are achieved to those using LESs (which are computationally expensive). Condition for this is the accurate calibration of the turbulent Schmidt number. For the scenarios investigated in this paper, the differences between the turbulent 2D RANS model and 3D LES models are not significant, i.e. the mass fraction at the two outlets and concentration distribution in the 1D mixing model were quite similar. The resulting 1D model has been implemented in a hydraulic solver and was validated in experiments on a medium-sized test network.

Future work will focus on laminar 3D simulations in order to improve the performance of the developed mixing module for laminar flow conditions. Furthermore, a mixing model for further four-way junctions (e.g., U-junctions) has to be derived. An interesting topic for further research is also the investigation of the effect of slowly varying flow rates and the comparison of the results with detailed 3D simulations and the simplified 1D model.

## ACKNOWLEDGEMENTS

The project is supported by the German Federal Ministry of Education and Research (BMBF project: 13N12180) and by the French Agence Nationale de la Recherche (ANR project: ANR-11-SECU-006). Further we thank Reik Nitsche for conducting the experiments at the Water Technology Center (TZW).