**1. Introduction**

Urban flooding events can cause significant economic and societal disruption. Numerous studies [1–3] have suggested that the occurrence of flooding in urban areas is likely to increase in the future due to increased urbanisation and changes in precipitation patterns, making intense rainfall events and the inundation of local drainage systems more common. The majority of urban flooding hazard studies focus on the economic damage, or direct risks to the public derived from hydraulic modelling of the depth and velocity of floodwaters resulting from historic or design rainfall events (see, e.g., in [4,5]). However, an increasing number of studies have also considered the public health risks of exposure to flood water, which may take the form of long term mental impacts [6], or illness from direct exposure of the public to contaminated flood water. Urban floodwater may contain a mix of rainwater, stormwater runoff and waste/foul water from surcharging urban drainage systems and therefore may contain harmful bacteria [7,8]. For example, ten Veldhuis et al. [9] sampled and analysed flood water from three urban flooding incidents in the Hague, the Netherlands in areas served by combined sewers. In the study, values of intestinal enterococci and *E. coli* were found to be 1 to 3 orders of magnitude higher than values for good bathing water quality according to the EU Directive 2006/7/EC.

Understanding the concentrations, transport and fate of harmful contaminants in urban floodwaters for effective health risk assessment is challenging [10]. Current state of the art urban flood risk models consider urban hydrological processes and utilise hydrodynamic principles to route resulting flows in both piped drainage and surface overland systems, with interaction (i.e., mass transfer) nodes such as manholes or gullies, which are commonly represented by weir or orifice equations [11–13]. Although flood model calibration and validation is often difficult due to a paucity of full scale data, such tools are generally considered to give tolerable predictions of flood depths and are widely used for risk evaluation and asset management [14]. Recently, Mark et al. [15] developed an approach to integrate an understanding of contaminant transport and health risk into flood models, utilising the 2D Advection Dispersion Equation to simulate the mixing and transport of wastewater surcharging from drainage systems within overland surface flow (assuming a constant pathogen level within the surcharging flow). However, such approaches can significantly simplify a number of processes concerning sources, transport, survival and transformations of harmful contaminants (e.g., see in [10]). The number of additional terms and associated parameters required to account for transport and fate processes exacerbate non-identifiability and equifinality issues which are a common problem for complex integrated models [16]. To develop a more robust understanding of health risks posed by urban flood waters, detailed information is required concerning individual processes associated with sources, transport pathways and life cycles of pathogens from sewer/drainage networks to surface flows and on urban surfaces. For example, recent studies have considered the behaviour of waterborne pathogens on different urban surfaces [17] and evaluation of pathogen levels in urban rainfall runoff flows [18].

However, as far as the authors are aware, no studies to date have considered the exchange of contaminated material (in soluble or particulate form) from drainage/sewer networks to surface flows during flood events via interaction structures such as gullies and/or manholes. Flows in and around surcharging hydraulic structures are highly complex and three-dimensional, especially during interactions with surface flood flows [19]. It is also likely that contamination concentrations within urban drainage/sewer networks will vary significantly as the proportion of stormwater and quantity and nature of contaminated material (i.e., dissolved vs. entrained solids) within the network varies during flood events. Numerous studies have considered the mixing of soluble material in manhole structures in the absence of interacting surface flood flows, demonstrating that mixing/transport (and thus mass exchange) processes are sensitive to geometrical characteristics and poorly described using commonly used simple models such as the 1D ADE which are commonly used to model pollutant transport and mixing in piped networks (see, e.g., in [20,21]). More complex 3D CFD based approaches have been shown to be able to quantify hydraulic and solute mixing processes in hydraulic structures such as manholes [12,22–25]. However, to date such models have not been experimentally validated in urban flood situations which include complex interactions between piped and surface flows [19]. While such 3D models are too computationally expensive to be used in direct design or network simulation, validated CFD models can be used to conduct experiments which may elucidate relationships and timescales describing the transport of materials to surface flows, understand the influence of geometric or hydraulic variables on mixing and mass transport characteristics or be used to calibrate simpler models.

Understanding how contaminants move from sewer networks to surface flows is a key aspect for understanding health risks posed by urban floods and possibly to foster the design of techniques to mitigate negative effects. This study conducts a detailed 3D numerical simulation of flow and soluble mass transport through a manhole during surface flooding conditions where net sewer to surface exchange flows are simulated. Whilst the focus of this study is limited to soluble pollutants only (i.e., those fully dissolved in the flow), it is recognised that the transport of contaminated solid material (e.g., fine sewer sediments) is also relevant in this context. The aims of the paper are to (1) compare the model outputs to new hydraulic and solute transport experimental datasets collected in a scale model facility designed to study interactions between pipe and surface flows. (2) Conduct numerical experiments to provide a more complete understanding of mass exchange to surface flows via hydraulic structures, including characteristic timescales associated with the occurrence of steady mass flow rate conditions.

#### **2. Materials and Methods**

Section 2.1 presents details of the setup used to gather experimental data to evaluate the numerical model. Section 2.2 provides a definition of key timescales and processes to be explored using CFD modelling and Section 2.3 describes numerical model and tests undertaken.

#### *2.1. Experimental Setup*

To collect data required for evaluating the numerical model, an experimental testing campaign was conducted using a physical 1:6 scaled model of a linked sewer/surface system, constructed at the University of Sheffield (Figure 1) [11,19,26–31]. The model is composed of a surface "floodplain" 8.2 m long, 4 m wide, constructed from acrylic (slope of 0.001 m/m). This floodplain is connected to a piped sewer system via a manhole with a diameter of 0.240 m (simulating a 1.440 m manhole at full scale, a size typical of UK urban drainage systems for pipes diameters up to 900 mm [32]). The sewer comprises a 0.075 m (internal) diameter clear acrylic pipe (simulating a 0.450 m pipe at full scale). To simulate flooding conditions, a series of steady flows were passed into the inlets at the upstream boundary of the sewer system and the floodplain. During each test, a portion of the flow within the piped network passed into the surface system via the manhole structure, with the remaining flow passing to the pipe outlet tank via the downstream boundary. The scheme of the facility is displayed in Figure 1. – ed of a surface " oodplain"

**Figure 1.** Scheme of the exchange structure showing the floodplain, sewer pipe and the manhole.

— — —QOutS re monitored and controlled via Labview™ The experimental facility was equipped with three electromagnetic flowmeters (two of them at the sewer and surface flow inlet—*QInS* and *QInF*—and one in the outlet of the sewer—*QOutS*) of 0.075 m internal diameter. The accuracy of the flow meters was validated using volumetric discharge readings at the laboratory measurement tank. A butterfly flow control valve was fitted to the pipe that feeds the sewer and the floodplain, calibrated such that steady inflows from 1 to 11 L/s can be set. Electromagnetic flowmeters and butterfly valves were monitored and controlled via Labview™ software. For all the tests conducted, flows were first established and allowed to stabilise before data values were recorded. Once established, data were collected for a period of 3 min to define reliable temporally averaged values for each flow mater. Mean steady state flow exchange rate through the

= −

manhole structure (*Qe*) was quantified based on mass conservation principles (i.e., *Q<sup>e</sup>* = *QInS* − *QOutS*). During the experimental campaign, water column pressure at the sewer outlet point was measured using pressure sensor (Figure 1). This sensor was calibrated to directly convert the output signal (mA) to gauge pressure and this procedure was conducted using a pointer gauge. The measure values were compared against defined calibration outcomes and errors were quantified to be ±0.69 mm within the water depth range of 0 to 600 mm. Values recorded with the pressure sensor were then gathered in real time by using the same Labview™ software described previously.

Experiments to understand solute transport and mass exchange were undertaken by injection of a neutrally buoyant soluble fluorescing dye (Rhodamine WT) into the sewer pipe >8 m upstream of the first measurement point (Cyclops 1 in Figure 1). The distance between the location of the injection and the measurement areas was higher than 10D (D = sewer pipe diameter) to allow cross sectional mixing [33]. Measurement of concentration vs. time profiles upstream and downstream of the manhole was conducted using Cyclops-7F™ fluorimeters. For this experiment, dye of concentration 10−<sup>3</sup> mg/L was fed into a constant head tank, from where injection into the sewer pipe was controlled by a manual open/close valve. For each test conducted, a 15 s duration pulse of dye of was introduced into the inflow pipe, and the resulting in-pipe concentrations monitored using the fluorimeters. The electrical sensor output was converted to concentration using experimentally predetermined calibration equations.

Experimental tests were conducted under steady state hydraulic conditions over a range of sewer inflows (*QInS*) and surface inflows (*QInF*), producing different flow exchange rates (*Qe*). Reynolds number (Re) for these tests ranged from 1.37 × 10<sup>6</sup> to 1.72 × 10<sup>6</sup> in the sewer inlet, which indicates a fully turbulent flow condition. Surface flow depths measured 350 mm upstream of the centreline of the manhole ranged between 5 to 17 mm over the tests conducted. Full details of these test conditions along with their numerical replication in CFD are presented in Table 1 in Section 3.

#### *2.2. Timescales and Mass Exchange Processes*

For a given pulse of soluble contaminant passing within a pipe network entering an exchange structure (e.g., a manhole) during sewer-to-surface flow exchange conditions (i.e., *Q<sup>e</sup>* > 0), a proportion will pass through the structure remaining within the pipe network and a proportion will exit to the surface flow. The change in total solute mass within the exchange structure at a given point in time can be expressed as

$$\frac{dM\_m}{dt} = \dot{M}\_{\rm PI} - \dot{M}\_e - \dot{M}\_{\rm PO} \tag{1}$$

where *dM<sup>m</sup> dt* is the rate of change in mass of solute within the exchange structure (mg/s), . *MPI* is the solute mass flow rate entering the exchange structure via the pipe network (mg/s), . *Me* is solute mass flow rate (mg/s) leaving the exchange structure to the surface flow and . *MPO* is the solute mass flow rate leaving the exchange structure via the pipe network (mg/s).

Considering the arrival of a well-mixed solute of concentration (CPI) at the inlet to the exchange structure under steady inflow conditions and Equation (1) above, a number of characteristic timescales can be defined.


mass flow through an inlet/outlet is a product of the rate of hydraulic flow rate and mean solute concentration, the proportion of solute mass exchanging though each outlet will become equivalent to the flow partition through the structure (Equations (2) and (3)).

$$\frac{\dot{M}\_{\mathcal{E}}}{\dot{M}\_{\text{PI}}} = \frac{\mathbb{Q}\_{\mathcal{E}}}{\mathbb{Q}\_{\text{PI}}} \tag{2}$$

$$\frac{\dot{M}\_{\rm PO}}{\dot{M}\_{\rm PI}} = \frac{Q\_{\rm PO}}{Q\_{\rm PI}}\tag{3}$$

Further experiments on specific hydraulic structures are required to understand the characteristic timescales (*t*1,2,3) and how these are affected by local flow characteristics. Flow structures and mixing processes in tanks and urban drainage structures have been studied previously but in the absence of surface flow interaction. For example, a general description of flow structures inmanholes under surcharged pipe conditions is given in [12]. When the sewer pipe inflow enters the manhole, three distinguished flow zone can be commonly observed [12,23,34]. A part of the pressurised flow, known as the diffusion zone, expands inside the manhole at a ratio of 1:5 towards the manhole diameter length. The remaining strong velocity zone forms a conical shape which has the same central axis as the inlet pipe. The slope of this cone is generally 1:6.2 towards the manhole length and travels through the manhole diameter towards the outlet. This conical form may create different distinctive scenarios based on the manhole to sewer pipe diameter ratio (φ*m*/φ*m*) and available surcharge depth (s). For 3.0 < φ*m*/φ*<sup>p</sup>* < 4.5 and with s > 0.2φ*m*, the core velocity region travels out of the manhole without contributing to the mixing process [23,35]. This is the most conventional size and surcharge depth characteristics for an overflowing manhole commonly seen the drainage systems and corresponds to the present study. In these cases, the diffusion zone is mainly responsible for solute mixing inside a manhole [20,23]. Part of the diffusion zone interacts with the manhole wall and travels upward. Later, this upward moving flow further divides in two components of which the first part exits through the surface and the last part recirculates within the manhole. However, how these structures interact with a surface flow, how effective they are in transporting solute mass to the surface and key timescales for well-mixed conditions (i.e., Equations (2) and (3)) are currently unclear and will be analysed in the current work using CFD techniques.

#### *2.3. CFD Modelling*

The hydraulics of the experimental model was reproduced using three-dimensional CFD modelling tools OpenFOAM® v.18.12 within interFoam solver [36–38], which considers the two-fluid system as isothermal, incompressible and immiscible utilising a Volume of Fluid (VOF) model [39]. Despite Larger Eddy Simulation (LES) models being known to model the turbulence structures of the flow more effectively, LES models are significantly more computationally expensive than those of RANS models. Moreover, RANS models are also reported in the literature for their accuracy in replicating manhole hydraulics properly and efficiently [29,34,40]. The model uses a single set of Navier–Stokes equations (Equations (4) and (5)) for both fluids with additional equations to describe the free-surface (Equation (6)). The interFoam within RNG k-ε Reynolds-averaged Navier–Stokes equations also requires Equations (7) and (8).

$$\nabla \cdot \mathfrak{u} = 0 \tag{4}$$

$$\frac{\partial \rho \mathbf{u}}{\partial t} + \nabla . (\rho \mathbf{u} \mathbf{u}) = -\nabla p^\* + \nabla . \pi - \mathbf{g}. \mathbf{x} \nabla \rho + f\_{\sigma} \tag{5}$$

$$
\frac{\partial \alpha}{\partial t} + \nabla \cdot (\alpha \mathfrak{u}) + \nabla. [\mathfrak{u}\_{\mathbb{C}} \alpha (1 - \alpha)] = 0 \tag{6}
$$

$$\frac{\partial \rho k}{\partial t} + \nabla \cdot (\rho k \mu) = \nabla \cdot (\Gamma\_k \nabla k) + P\_k - Y\_k \tag{7}$$

*Water* **2020**, *12*, 2514

$$\frac{\partial \rho \varepsilon}{\partial \mathbf{t}} + \nabla \cdot (\rho \varepsilon \mathbf{u}) = \nabla \cdot (\Gamma\_{\varepsilon} \nabla \varepsilon) + P\_{\varepsilon} - Y\_{\varepsilon} + D\_{\varepsilon} \tag{8}$$

Where *u* is the mean velocity vector in the Cartesian coordinate; ρ is the density of the fluid mix; *g* is the acceleration due to gravity; *t* is the time; τ is the shear stress tensor; *p* ∗ is the modified pressure adapted by removing the hydrostatic pressure from the total pressure; *f*<sup>σ</sup> is the volumetric surface tension force (where CSF and interface curvature are included); α is the VOF function; *k* is the turbulent kinetic energy; ε is the energy dissipation; Γ*<sup>k</sup>* and Γ<sup>ε</sup> are the diffusion for *k* and ε, respectively; and *<sup>P</sup>*, *<sup>Y</sup>* and *D* are the Production, Dissipation and Additional term for RNG, respectively.

In this work, an additional solute transport model was added to the interFoam VOF model. The main advection–dispersion equation used in the model is

$$\frac{\partial \mathbf{c}}{\partial \mathbf{t}} + \nabla \cdot \left(\mathbf{u} - v\_s \frac{\mathbf{g}}{|\mathbf{g}|}\right) \mathbf{c} = \nabla. (\alpha v\_t \nabla \mathbf{c}) \tag{9}$$

where *c* is the solute concentration of the flow, *v<sup>s</sup>* is the terminal velocity due to gravity (which is zero for a neutrally buoyant solute) and ν*<sup>t</sup>* is the turbulent kinematic viscosity of water, which is a function of the turbulence of the flow [41] and taken to be equivalent to diffusivity [42,43]. The multiplication of ν*<sup>t</sup>* by α prevents solute particles from entering the air phase [42].

Earlier model validation works by the authors presenting measured velocities using PIV within the same experimental facility [29] showed that RNG k-ε model is a suitable RANS modelling choice for predicting water elevation and velocity profiles and hence is chosen for this work. This turbulence model can also capture complex flow and is known for better performance for separating flow [22,23,29]. Apart from wall boundary condition, five open boundaries were prescribed in the model: two inlet and two outlet boundaries at the sewer pipe and floodplain, respectively, and an atmosphere boundary at the floodplain (Figure 1). The inlet boundaries were prescribed as fixed velocities, while the outlets were applied as fixed pressure. This measured temporal mean pressure data was used for the sewer outlet pressure boundary condition (measured at *POutS*). The atmosphere boundary was set as equal to atmospheric pressure and zero gradient for velocity to have free airflow if required. All the wall boundaries were prescribed as noSlip condition. The sewer pipe walls were considered as rough wall applying equivalent sand roughness height (*ks*). Further details of measured head losses within the experimental facility can be found in [26].

Cfmesh v1.1 [44] was used to generate the hexahedral computational meshes, keeping the maximum mesh size as 10 mm towards all three Cartesian directions. The boundary meshes were kept small in such a way that 30 < y <sup>+</sup> < 300, keeping three boundary layers at the all wall boundaries. A standard wall function was applied to all the walls, which has been shown to be appropriate for the application of boundary turbulence effects for such mesh sizes [36], eliminating the necessity of fine layered boundary meshes. Figure 2 shows part of the computational mesh created for this work. The rest of the CFD model such as the choice of different meshes, different solvers parameters and solution schemes were obtained from another CFD model validated in an earlier work depicting the same experimental set-up [23,29]. The maximum Courant–Friedrichs–Lewy (CFL) number was kept as 0.9. The cluster computing system at the University of Coimbra was used to run the simulations using MPI mode. Each simulation was run for 300 s to reach steady state conditions. For comparisons with experimental datasets, the measured solute concentration for each test condition was applied through the sewer inlet pipe at Cyclopes 1 when the hydraulic model reached a steady state. Unsteady model results were saved at every 0.01 s interval. Model solute concentrations were extracted at different sections and compared with the experimental measurements.

**Figure 2.** Computational mesh for the study showing the manhole and its connections to sewer pipes and surface.

– Experimental tests 1–6 (including repeats) listed in Table 1 were replicated to perform calibration and testing of the CFD model. In these tests, solute concentration in the inflow pipe was taken as the measured value recorded at Cyclopes 1. Measured and predicted concentration curves are compared at the downstream measurement point (Cyclopes 2) along with measured and modelled flow rates in the pipe and exchanging to the surface (i.e., flow partition).

− – *—* To isolate and understand the effects of the manhole separately from the pipe network, the calibrated model was further applied to Test 1-A, 2-A, 3-A and 4-A under identical hydraulic conditions to tests 1, 2, 3 and 4, respectively. However, in these tests, solute concentration was applied uniformly at the sewer pipe/manhole inlet boundary at a steady concentration of 1.2 × 10 <sup>−</sup><sup>6</sup> mg/L. Resulting concentration time series were extracted at the manhole–pipe junctions (section A and B in Figure 1) and at the exit to the surface flow (section C in Figure 1). Solute mass flow rates to the surface flow and the downstream pipe as well as characteristic time scales of the manhole as described in Section 2.1 were calculated from these tests (described further in Section 3.3). Test 1 was further extended numerically by changing the downstream boundary pressure; by decreasing 9.5% (Test 1-B) and by increasing 15% (Test 1-C) to enable lower and higher pipe to surface exchange flows (1.21 L/s and 2.43 L/s), respectively.

#### **3. Results**

#### *3.1. Model Calibration and Validation—Hydrodynamics*

− – – Calibration of wall roughness (*ks*) was performed based on experimental test results of flow exchange through the manhole to the surface. Applying a higher *ks* in the sewer pipe leads to lower flow through the outlet pipe with higher flow exchange from the manhole to the surface, and vice versa. The experimental values from Test 4 were used for calibration purposes as it had a sewer inlet flow which was median to all the sewer flows tested herein. Modelled *ks* values ranging from 1 × 10 −6 to 1 × 10 <sup>−</sup><sup>3</sup> mm were simulated in the CFD model. Results showed that *ks* = 0.0005 mm gives a comparable modelled value of the flow partition to the experimentally observed values (*Q<sup>e</sup>* within 1.7%). This value of *ks* is valid for smooth surfaces such as acrylic which is appropriate to the experimental setup used here. The same *ks* value was applied to the rest of the hydraulic simulations (Tests 1–3 and 5–6) for model validation. Table 1 compares experimentally measured and modelled steady state flow rates in the pipe and exchanged to the surface (*QOutS* and *Qe*) for each test, along with measured boundary conditions and calculated Reynolds numbers. Modelled and measured flow rates are found to be

−

within 1.7% in all cases. Figure 3 presents resulting calculated velocity streamlines and vectors within the manhole during Test 4.


**Table 1.** Experimentally observed and numerical flow rates for each test case. Solute injections for tests 4–6 were repeated 3 times. –

\* Test 4 data was used for hydraulic calibration.

**Figure 3.** Hydraulic conditions inside the manhole during Test 4. (**a**) Streamline of the flow indicating a general circulation pattern, (**b**) mean velocity vectors at the top horizontal plane of the manhole and (**c**) mean velocity vectors at the horizontal plane passing through the sewer pipe axis of the manhole. At all cases, main flow direction is from right to the left.
