Effect of surcharge on gully-manhole flow

Abstract In an urban drainage system, surface system and buried system are connected through several linking elements like gullies and manholes. To model an urban flood through a 1D or 2D routing model, proper representation of these linking elements is mandatory. This paper focuses on the effects on inline manhole head loss coefficient and gully flow due to change in manhole surcharges to get knowledge of their characteristics, flow patterns and influence in the overall flow. The solver interFoam with VOF capability under open source CFD toolbox OpenFOAM® is used for the numerical simulations. The CFD model created is a real scale model following the experimental setup installed at the University of Coimbra. The CFD model results were compared with discharge and water depth data at the manhole and velocity data at the gully, measured from the physical model. The manhole head loss coefficient was found very high at low surcharge level and reaching a constant value of 0.3 for higher surcharge levels. The gully flow was found to follow typical orifice flow equation with different discharge coefficients for different water depth conditions of the manhole. Three different surcharged zones were identified at the manhole on which discharge coefficients are dependent. The shear stress pattern at the manhole floor is also found varying due to surcharge depth as well as the gully flow. The variation showed distinctive configuration at different surcharge zones.


Introduction
The Urban drainage system is termed as the main conveyance route for draining extensive rainfall and flood in a busy city. The system is composed of many linking elements that convey the flow during a rainfall event from the major or surface system to the minor or buried system. Manhole and gully are two very common drainage structures where the flow is normally complex, highly turbulent and possibly multiphase. When entering to a manhole or gully, the drainage flow undergoes a sudden expansion and/or contraction with several hydraulic losses, which are important in flood routing simulations. Most urban flood routing models are one or two dimensional and therefore cannot represent the complex flow pattern of these structures (Leandro et al., 2009b). These structures usually are not modelled but connected as a point entity and translated by discharge and head loss coefficients calculated using empirical equations (Leandro and Martins, 2016). State-of-the-art routing models in urbanised areas are Dual Drainage (DD) models as they simulate both surface flow and flow in buried pipes simultaneously among other urban key features (Djordjević et al., 2004;. These models use discharge coefficients to connect the two systems through linking elements. However, they also have weaknesses in considering linking elements as very few existing models are available to calibrate these coefficients (Djordjević et al., 2005). The hydraulic losses and efficiencies are very much dependent on the flow conditions and should be studied in order to allow the translation of the different structures behaviours throughout the possible flow conditions in a model. At low flow conditions, drainage pipes are partially filled and the estimation of the head loss in a manhole can be calculated considering an open channel flow. But in the case of flood events, the system becomes pressurised, uses its storage potential and creates surcharged flow conditions. In those cases, estimation of head loss is complex, involving more structural and hydraulic factors of the system. The flow through the gully can also vary due to the surcharge heights of the manhole. This paper aims to determine manhole head loss coefficients and gully discharge coefficients at different surcharge levels of a manhole.
Researchers pointed out many factors that influence the head loss in a manhole using experimental model facilities. The shape of the manhole is an important factor in determining the flow conditions. Four different types of manholes are mainly found in different countries (Asztely, 1995), which are known as Type I, II, III and IV. Type I has the simplest design with no guided channel. Type II and III has guided http://dx.doi.org/10. 1016/j.jher.2017.08.003 channel where the channel is half circular invert and U-type invert respectively. Type IV is a modified form of type II, where the manhole bottom is sloped towards the guided channel. Out of different types of the manhole, Type I is the most investigated. Sangster et al. (1958) studied Type I surcharged manhole where the inlet pipe invert was at the same level with the manhole bottom. The head loss coefficient varied between 0 and 0.3 with different inlet-outlet pipe diameters ranging from 76 to 145 mm. A unsystematic variation on energy losses at small manholes was reported by Sangster et al. (1958) and Lindval (1984). Lindval (1984) also found the head loss coefficients varying between 0.08 and 0.88 by using different manhole and inlet pipes with manhole to inlet diameter ratios ranging from 1.7 to 4.1. Many authors reported an oscillatory surface with high head loss coefficients at very low surcharge conditions (Howarth and Saul, 1984;Lindval, 1984;Pedersen and Mark, 1990). The coefficient of head loss was as high as 1.3 for reduced surcharge flow, as reported by Arao and Kusuda (1999). Pedersen and Mark (1990) investigated Type I, II and III manholes with scaled physical models with a pipe to manhole diameter ratios of 1. 22, 2.11, 3.20 and 4.94. They reported the loss coefficient to increase with manhole to pipe diameter ratios. They also explained this phenomenon through submerged jet theory; originally described by Albertson et al. (1950). Guymer et al. (2005) and Lau et al. (2007) checked the velocity distribution of surcharged manholes and described their results in line with submerged jet theory.
Several researchers studied surcharged manholes using CFD. Stovin et al. (2008) explained the approach to validate a CFD manhole model. One major issue in CFD modelling of drainage structure is locating the position of the free surface. Two different techniques are commonly followed for modelling the free surface; namely: Rigid lid approximation and Volume of Fluid (VOF) method. In the Rigid lid approximation, the free surface is assumed as a frictionless lid at a fixed positioned. Lau et al. (2007) and Stovin et al. (2013) used Rigid lid approximation using CFD modelling tool Fluent (ANSYS Ins, 2009) to simulate scaled surcharged manhole behaviour and to predict solute transport through it. However, rigid lid approximation requires prior knowledge of the exact location of the free surface. In Volume of Fluid (VOF) method (Hirt and Nichols, 1981), the position of the free surface is not approximated but modelled by introducing the concept of volume fraction. This is considered as the recommended procedure when the position of the free surface is not known (Bennett, 2012;Lau, 2008). However, not much work has been done to model a manhole flow using VOF method. Moreover, all studies mentioned here used scaled manhole models. Uses of prototype scale manholes are rare.
In compared to manholes, fewer works have been done on the hydraulics of a gully. A typical gully collects runoff from the roadside curbs and delivers it to a nearby manhole. The gully outlet is normally placed at a higher location than the manhole floor. So at a low surcharge condition, the gully outflow enters the manhole as a plunging jet. This could be related to the typical hydraulics of a drop manhole structure. Some researchers focused on the hydraulics of a drop manhole but mainly discussing the energy dissipation of the drop (Arao et al., 2012;Carvalho et al., 2013;Carvalho and Leandro, 2012;Chanson, 2004;Granata et al., 2011Granata et al., , 2014. Martins et al. (2014) described the hydraulic performance of the gully at drainage condition. They used weir and orifice equations to explain the flow through the gully outlet. Carvalho et al. (2011) investigated flow through a gully using VOF based 2DV numerical model. They compared water depths, streamlines, velocities, pressures and the inlet discharge coefficients at different outlet sizes and positions. Gómez and Russo (2009) and Lopes et al. (2016) described flow efficiency of transverse gully grate using experimental and CFD models respectively. Although the main task of a gully is to drain surface flow to the underground system, the opposite directional flow may occur due to any exceptional event when gully becomes pressurised and surcharged water flood the surface area. Leandro et al. (2014) and Lopes et al. (2015) described the reverse flow in a surcharged gully. Romagnoli et al. (2013) measured the velocity of a typical gully at reverse flow condition using Acoustic Doppler Velocimetry (ADV). Few research has been reported on effects of manhole surcharge on gully flow. Noh et al. (2016) used scaled urban drainage model along with Monte-Carlo simulations to estimate discharge coefficients between different sewer network components. But the research did not consider other structural elements such as manhole types and surcharge heights corresponding to the height of gully as these are responsible for different coefficient values.
In this work, the effect of manhole surcharge on manhole head loss coefficients and manhole-gully discharge coefficients have been studied. The analysis has been done using open source CFD tools OpenFOAM®. The CFD model replicated the real scale experimental facility installed at the hydraulic lab of the University of Coimbra. The work uses VOF model to examine the flow behaviour of the drainage system. Numerical model data was compared with discharge, water depth or pressure data measured at the manhole and ADV velocity data measured at the gully. This work describes the experimental and the numerical model at the 'Methodology' section. In the 'Results and Discussion' section, comparisons were drawn between numerical and experimental works. Attempts were also made to describe the manhole and gully flow using basic hydraulics. The last section concludes the work.

Experimental setup
The Multiple Linking Element (MLE) installation at the hydraulic laboratory of University of Coimbra was used to validate the CFD modelling procedure. The facility has been described in other literature such as Leandro et al. (2009a) and Carvalho et al. (2013). A part of the facility has been used in this work containing a rectangular surface drain (0.48 m wide) with a slope of 1:1000, a rectangular gully pot (0.6 m × 0.24 m × 0.30 m) and a circular manhole (1 m diameter) (Fig. 1). The manhole is a "Type I" (Asztely, 1995) manhole with no guided channel. A 300 mm horizontal inlet-outlet pipe is connected to the manhole, whose invert level is at a height of 0.10 m from the manhole floor. The manhole to inlet-outlet pipe diameter ratio (Φ m /Φ p ) is 3.33. The manhole does not receive any flow from its top. A lid is loosely connected that ensures equal pressure between inside and outside of the manhole. The gully is connected to the manhole with an angular pipe of 80 mm diameter, at an angle of 63°in plan and 90°in vertical.
The gully-manhole element replicates a typical drainage system in Portugal. Data from three electromagnetic flow meters were used to define boundary conditions and to verify simulations. The first one is located at the upstream of Drain Inlet, the second one is the referred at the pipe between manhole and gully and the third is at the downstream of Pipe Outlet; from which different discharges at the drain, gully and manhole can be measured (Fig. 1). The first and the third flow meters are just outside the computational domain. The second flow meter is within the domain and it is represented by a small contraction in the numerical model. The contraction zone has a diameter of 60 mm. Three pressure sensors are also installed in the setup located at the outlet pipe shown in Fig. 1 as P18, P22 and the manhole centre, shown as P23.

Experimental test cases
The grate of the gully was removed throughout the experiment. Two experimental data sets have been used for the numerical model validation. In the first set (SE1), flow through the drain and gully was observed and measured giving 19.8 l/s flow at upstream of the Drain Inlet. The flow through the rectangular drain upstream of the gully was a free surface flow with a depth of 0.068 m; which yielded a Froude number of 0.6. A part of the incoming drain flow passed through the gully outlet and entered the manhole, the rest being overflown through the drain and falling as a free fall to the reservoir tank through the Drain Outlet ( Fig. 1). The intercepted flow by gully enters the manhole as a free fall plunging jet generating a recirculation zone in the manhole. This flow mixes with the inflow coming through the Pipe Inlet and drains out through the Pipe Outlet. The gully was given special attention. The velocity fields at the gully were measured using Acoustic Doppler Velocimetry (ADV) at three vertical planes parallel to the longitudinal flow direction of the drain. The first and the third planes (Plane 1 and Plane 2 of Fig. 2 left panel) are at 0.05 m distance from the two longitudinal walls of the gully; which made each of the planes 0.07 m apart from the central line of the gully. The second plane is the central plane (Plane C). Each plane had 121 point velocity measurements, at an interval of 0.050 m and 0.025 m towards horizontal and vertical directions respectively ( Fig. 2 right panel). Since the drainage capacity of the gully was lower than the flow in the drain, the water level was higher than the top level of the gully. The velocity measurement was taken inside the gully only and was not extended to the water surface. Detail description of the measurement can be found at Beg et al. (2016b).
In the second set of experimental works (SE2), only flow in the pipe passing through the manhole inlet was used and recorded using a discharge flow meter. The three pressure sensors (P23, P22 and P18 shown in Fig. 1) were also used to check the pressures of the flow at different sections. The inflows and the manhole surcharges were controlled using two valves at the inlet and outlet pipe. The pressure sensors measured piezometric pressures for both free surface and pressure flow conditions, which were converted to piezometric head considering the bottom of the manhole as the zero datum. Eighteen experimental runs were performed, keeping the outlet valve opened at 30% (low flow), 40% (medium flow), and 60% (high flow) to observe different flow velocities inside the pipe and the manhole. All discharges and corresponding water pressures throughout the system were recorded at a steady flow condition. This setup was utilised at different combinations of inflows ranging from 12 l/s to 145 l/s and manhole water level up to 1.40 m (Beg et al., 2016a).

Numerical model description
The objective of the numerical modelling validation with experimental measurements was to characterise the incoming flow through the gully and check the flow path in the manhole during drainage condition. Open source three dimensional CFD model tools OpenFOAM® v2.3.0 was used in this study. The solver interFoam was chosen which includes Volume of Fluid (VOF) method (Hirt and Nichols, 1981) to track the free surface or interface location between two fluids. This method uses volume fraction indicator function α to determine the amount of liquid present in each cell. In the case of = α 0 or 1, the cell volume is considered filled with air or water respectively; while 0 < α < 1 represents that the cell contains the free surface as it is partially filled with water.
The interFoam solver uses a single set of Navier-Stokes equations for the two fluids independently for water and air, and additional equations  Journal of Hydro-environment Research 19 (2018) 224-236 to describe free-surface where the velocity at free-surface is shared by both phases. The solver considers a system of isothermal, incompressible and immiscible two-phase flow. The model deals with Reynolds averaged conservation of mass and momentum expressed according to Rusche (2002) as: where g is the acceleration due to gravity, u is the velocity vector in the Cartesian coordinate, τ is the shear stress tensor, * p is the modified pressure adapted by removing the hydrostatic pressure (ρg·x) from the total pressure, ρ is the density and f σ is the volumetric surface tension force.
The viscous stress term is defined by the Newton's law of viscosity for incompressible fluid, The advection equation to describe free-surface in VOF method uses an interfacial compressive term to keep the interface region confined in a small space (Rusche, 2002;Weller, 2002) and is described as: The last term of Eq. (4) is the compressive term. The term − α α (1 ) ensures that the compressive term or compressive velocity u c is calculated only at the interphase (when 0 < α < 1). This velocity acts at the perpendicular direction to the interface and defined as: The volumetric surface tension f σ , shown at Eq. (2) is calculated by the Continuum Surface Force model (Brackbill et al., 1992).
where κ is referred as the surface curvature.
To model the turbulence phenomena considering RANS approach, two different − k ε model families were considered for the gully and manhole flow. For the simulations considering gully outlet flow, importance was given to the pipe flow at the gully outlet and Standard − k ε turbulent modelling approach was used as this is considered to give good prediction of fully turbulent flow (ANSYS Inc, 2013; Ghoma, 2011). For the manhole head losses simulations, resolving the complex flow inside the manhole was considered essential. Therefore, the Renormalization group of − k ε turbulent model (Yakhot et al., 1992) was used for these simulations. This approach is described to be more accurate to predict complex shear flow with separation  and more responsive to the effects of rapid strain and streamline curvature (ANSYS Inc, 2013;Bennett, 2012;Carvalho et al., 2008). Both of these turbulence calculation approaches use two closure equations for k (turbulent kinetic energy) and ε (Energy dissipation).
The dynamic viscosity (μ) is calculated as: where ν 0 and ν t are molecular viscosity and turbulent viscosity respectively.

Mesh generation
Three different computational meshes were simulated with different grid sizes to check the mesh grid convergence using RNG − k ε turbulent model. The three meshes have dx = 6.6 mm (Mesh 1), 20 mm (Mesh 2)  and 25 mm (Mesh 3) respectively. The number of cells in these meshes were 8.64 million, 310 thousand and 156 thousand respectively. Richardson extrapolation was considered in the mesh convergence study following Celik et al. (2008). The global refinement ratio was 3.8; which is bigger than the recommended minimum value of 1.3. The three meshes were simulated with 30 l/s inflow and 0.6 m of surcharge depth in the manhole. The axial velocity profile at the centre of the manhole and the centre of the outlet pipe were compared for the mesh analysis. A total of 38 point velocities were compared. The longitudinal velocities from the three meshes are shown in Fig. 3.
To calculate the mesh convergence study, different parameters of the three meshes were calculated. The apparent order p was calculated for each 38 velocity points by the iteration of the following equations: where e a 21 is the approximate relative error between Mesh 2 and Mesh 1 and calculated as: where Φ 1 and Φ 2 are the data from Mesh 1 and Mesh 2 respectively. In the analysis, oscillatory convergence was found at 22 points (57%). The approximate relative error close to the wall was found 53% comparing mesh 1 and 2, and 30% when comparing mesh 2 and mesh 3. However, at the centre where the jet stream is located, the approximate relative error was found 7.24% comparing mesh 1 and 2, and 6.2% comparing mesh 2 and 3. Large uncertainty close to the wall was found due to high gradient of velocity. The average GCI for mesh 2 was found 16%. From Fig. 3 (left panel), it can be seen that the fine mesh of 6.6 mm creates different flow structure close to the manhole bottom. While the remaining two meshes (20 mm and 25 mm) creates similar flow in the manhole. The location of vortex formation was changed at different meshes. A high percentage of oscillatory convergence, as well as change in vortex formation, resulted in high GCI values. As for this work, the focus is given to manhole head loss coefficient k, the value was checked at all the three meshes. The k value in mesh 1, mesh 2 and mesh 3 were found 0.321, 0.324 and 0.371 respectively. It can be said that although mesh 2 showed different flow structure on the small scale compared to mesh 1, considering the flow at the large scale, mesh 2 gives similar results for head loss coefficient, k. Considering the simulation time required, Mesh 1 demands very high computational cost as the number of cells is significantly high. Considering the accuracy level and computation costs, the cell size of mesh 2 is chosen (cell  size = 20 mm) for this work. The maximum cell size is kept 20 mm towards all the three directions. The mesh was further refined at the walls and joints of different geometrical shapes. The created final computational mesh has 821 500 computational cells with a little more than 1.01 million nodes. Some of the mesh properties can be seen in Fig. 1 and Table 1.

Definition of boundary conditions and simulation control parameters
Six open boundaries were prescribed to the numerical model ( Fig. 1): Drain Inlet, Drain Atmosphere, Drain Outlet, Pipe Inlet, Manhole Atmosphere and Pipe Outlet. The closed boundaries were prescribed as walls. The Drain Inlet was further divided into two parts for the incoming water and air phases respectively. The boundary data were calculated from the experimental model completed beforehand.
As the simulation used − k ε turbulent modelling and RNG − k ε turbulent modelling, OpenFOAM® requires six types of Boundary Conditions (BC) for each boundary. They are alpha.water (water fraction in each cell volume), U (velocity vector in Cartesian domain), p_rgh (relative bottom pressure corresponding to datum), k (turbulent kinetic energy), ε (energy dissipation) and nut (turbulent viscosity). The first three BC's are required for hydraulic modelling while the last three are required for turbulence calculation.
For both inlets, fixed velocities/discharges were applied using alpha.water and U. Pressure data (p_rgh) were prescribed at the outlet boundaries. Both of the atmosphere boundaries were kept as zeroGradient velocity and relative air pressure as zero; so that air could be exchanged if necessary. All the wall BC's were kept as no-slip condition (i.e. velocity = 0). For the turbulent approach, values of k, ε and nut were calculated using the equations in FLUENT manual (ANSYS Ins, 2009), considering medium turbulence at the gully and manhole. All the walls are prescribed as wallFunction as this eliminates the necessity of fine layered boundary mesh and hence reduce the computational time (Greenshields, 2015). The numerical model represents the experimental facility, which was made of acrylic. The roughness was considered not important at the acrylic surface and thus was not taken into account at the wall BC's.
The model was ready to run after the boundary setup. During the simulation, adjustableRunTime was used keeping maximum CFL number to 0.8. Cluster computing system at the University of Coimbra was used to run the simulations at MPI mode using 16 processors. Each simulation took 40 s to reach a steady state; which took 138 h of effective computational time. For each case scenario, results were saved for 5 s of simulation run (considering the simulation reached steady state already) at an interval of 0.5 s, making a total of 11 time steps. Each analysis was done with the averaged result of these 11 time steps.

Comparison with experiments
During the experimental model run, pressure levels were recorded at different lengths of the pipe as well as the bottom centre of the manhole. The numerical model was validated with experimental data of flow depth/pressure head and discharge inside the structure. The comparison of the pressure and discharge has been presented in an earlier study (Beg et al., 2016a). The comparison was checked with the pressure data at transducer P23, P22 and P18 (see Fig. 1). Fig. 4 shows the comparison between experimental and numerical simulations in a pressure head vs discharge plot. The dot markers of different colours show data from different experimental models while the pentagons show data from numerical models. Fig. 5 shows the comparison between the two data sets in a scatter plot.
It can be seen from Fig. 5 that the numerical model showed very good match with the recorded pressure head of the experimental data. A maximum of 7% error was reported in pressure data; which shows that the numerical model can give a very good estimation of water depth/pressure head.
The numerical model results were also compared with the experimental study data of the velocity profiles at the gully obtained by the Vectrino acoustic velocimetry. Fig. 6 shows contours at the three different planes and Fig. 7 shows a comparison of longitudinal and vertical velocity profiles at different locations of the gully.
It can be seen from Figs. 6 and 7 that the numerical model shows good agreement with experimental data. The vertical vortex size and location created in the numerical result shows similarity to those observed in the experimental model data. Average statistical comparison between the two data is also shown in Fig. 7. It shows that the model can reproduce the longitudinal velocity component (V x ) very well (average r 2 = 0.95 and BAIS = −0.004 m/s). The representation of vertical velocity component (V z ) in the gully is at a satisfactory level (average BIAS 0.011 m/s and r 2 = 0.56). Fig. 8 shows the transverse velocity (V y ) in the gully at three different transects. It can be seen that transverse velocity is very low, in the range of −0.03 m/s to +0.03 m/s. This velocity component was found insignificant to compare with experimental results.

Velocity field and manhole head loss at different surcharge
The velocity field at the manhole was checked for different inflow  M.N.A. Beg et al. Journal of Hydro-environment Research 19 (2018) 224-236 and surcharge conditions. For this part, different numerical simulations were checked to compare the surcharge effect on the manhole. The gully inflow was disconnected. Four different inflows were applied through the Pipe Inlet with combinations of different piezometric pressure heads at the Pipe Outlet for various numerical simulations (see Fig. 1). The applied inflows were 30, 60, 90 and 120 l/s with different piezometric pressure head ranging from 0.4 m to 1.5 m. Different surcharge levels were obtained from the simulations, ranging from 0.20 m to 0.85 m. Fig. 9 shows velocity contours at both horizontal and vertical planes through the main axis of inlet-outlet pipes from four simulations. The four contours are from simulations having 30 l/s and 120 l/s inflows at the Pipe Inlet (lowest and highest among the simulations) with and 0.45 m and 1.2 m piezometric pressure heads at the Pipe Outlet. The corresponding water depths found at the manhole centre are 0.45 m, 1.21 m, 0.54 m and 1.27 m respectively. Fig. 9 shows that for all the cases, highest velocity can be observed at the central axis of the pipe. The main core velocity does not change much during its travel through the manhole but reaches to the outlet pipe. This denotes a 'short circuiting' as reported by other authors (Stovin et al., 2013). The longitudinal velocity shows symmetric view at the horizontal cross sections at high surcharge conditions. At the vertical sections, the longitudinal velocity of the jet core decreases with the increase in vertical distance from the central line of the pipe. However, at the bottom of the manhole near the outlet pipe, the longitudinal velocity is found very low. At lower surcharge conditions (a1 and c1), almost the whole vertical section has high longitudinal velocity towards the outlet. Comparing a2 and c2 conditions, it can be observed at c2 with higher discharge, the jet core does not show a symmetric pattern. On the other hand, at higher surcharge conditions (b1 and d1), the longitudinal velocity has been found towards the opposite at the upper part of the manhole. The horizontal sections show that the longitudinal velocity turns towards the opposite direction near the manhole wall at both higher and lower surcharges. Comparing b1 and d1 of Fig. 9, it can be seen that the longitudinal dispersion of the jet velocity increases with increasing discharge. The horizontal section with higher inflow and surcharge (d2) shows that unlike c2 case, the velocity is symmetric at a higher surcharge.
Head loss coefficient of the manhole was also checked at different manhole surcharges (s). The surcharge levels are calculated as the difference between water levels and outlet pipe soffit levels. From each simulation, the surcharge ratio was calculated as the ratio between surcharge level and inlet pipe diameter (s/Φ p ). The head loss coefficient (k) was calculated as a ratio between pressure head difference and velocity head. The comparison is shown in Fig. 10. Fig. 10 shows that at higher surcharge, the coefficient stays fairly around 0.3 for all discharges. But the head loss coefficient rises very high when the surcharge ratio goes below a certain limit. The change occurs at around s/Φ m < 0.2. The line s/Φ m = 0.2 is shown with a dashed line in Fig. 10; which is also equivalent to s/Φ p = 0.67, as Φ m / Φ p = 3.33. This phenomenon can be explained using submerged jet theory described by Albertson et al. (1950). Several authors (Guymer et al., 2005;Pedersen and Mark, 1990;Stovin et al., 2013) applied the submerged jet theory for characterising manhole head loss coefficient. According to Authors, the incoming flow from the inlet pipe enters the manhole as a submerged jet. The jet core expands after entering the manhole and creates two jet regime; a) core region and b) diffusion region (Fig. 11). Guymer et al. (2005) and Stovin et al. (2013) verified that the diffusive region diverges from the inlet pipe with the travelling lengths inside the manhole at a ratio of 1:5 and the core region keeps the high jet velocity and maintains a conical shape. The core region diminishes as the jet travels through the manhole and remains until the distance of 6.2Φ p , where Φ p is the diameter of the inlet pipe.
For an inlet pipe of 300 mm, the core region of the submerged jet can travel (6.2 × 0.3=) 1.85 m before diminishing completely. As in this work, the manhole diameter is 1.0 m, the core region of the jet should enter directly to the outlet before diminishing completely. This condition can be compared with Fig. 9.
The diffusive region expands from the line of axis of the inlet pipe. The downward facing diffusive region can expand for 0.10 m towards the bottom of the manhole then reaches the manhole floor. Upon reaching the manhole floor, the flow reflects back to upward direction. This phenomenon occurs for all the applied discharge and surcharge combinations. But the upward facing diffusive region may vary according to the available surcharge depth. As the region expands to a ratio of 1:5 towards vertical to horizontal, it may expand to 0.20 m towards the vertical direction until it reaches the manhole wall at the opposite side. But in case, the available manhole surcharge is less than 0.20 m, the diffusive region cannot develop completely. It interacts with the surface violently and creates additional head losses. This phenomenon is reflected in Fig. 10 as the head loss coefficient increases drastically in the cases where the available surcharge is low (s/ Φ m < 0.2). Although, the change of head loss coefficient at different inflow conditions do not follow any particular trend. As for example, the head loss coefficient at 60 l/s at low surcharge was found less than those of at 30 l/s. Similarly, at surcharge ratio within 0.1 to 0.2, flow rate of 90 l/s made higher head loss in compared to those of at 120 l/s. However, in other low surcharge simulations, the higher discharges showed higher head losses. No justification could be drawn for these variations.

Effect of manhole surcharge on the gully flow
The flow distribution from the gully to the manhole at different surcharge heights were checked from the numerical models. A schematic diagram of the gully-manhole system is drawn in Fig. 12. The   11. Velocity distribution and diffusion region in circular free jet (adapted from the works of Guymer et al. (2005) and Pedersen and Mark (1990)). depth of the gully was 0.315 m with the bottom level at 1.3 m. This made the bottom level of the drain 1.615 m (Z). The water depth in the drain at the immediate upstream of the gully was 68 mm (h), which made the water level at the drain 1.683 m (Z + h). The gully outlet pipe has a diameter of 80 mm with a contraction for discharge meter at the end of the pipe, which has a diameter of 60 mm. The soffit level of the outlet pipe is at a level of 0.945 m (Z o ). The manhole water depth was increased from 0.53 m to 2.2 m (H) at different numerical simulations and corresponding changes in the intercepted flow of the gully (Q) were checked. In all the numerical simulations, the pressure flow through manhole inlet (Q 1 ) and surface flow through the drain inlet (Q 3 ) were kept the same; 42.2 l/s and 19.8 l/s respectively.
Examining the discharge from the gully, the manhole can be divided into three zones according to different surcharge levels.
1. Zone 1: the surcharge water level in the gully is less than the soffit level of the gully outlet pipe (H < Z o ). At this zone, the gully flow falls into the manhole as a plunging jet. 2. Zone 2: the surcharge level is more than the gully outlet pipe soffit level but less than water level at the drain (Zo < H < Z + h). At this zone, the gully flow enters the manhole as a submerged jet. 3. Zone 3: The surcharge level of the manhole is more than the water level at the drain. At the zone, the gully is unable to drain itself. The surcharge flow from the manhole pushes the water back to surface.
The discharges from the gully to the manhole at different conditions are shown in Fig. 13. Fig. 14 shows the gully jet at three different surcharge zones mentioned above.
It can be seen from Fig. 13 that the decrease in intercepted flow with the increase in surcharge height is not linear. From the simulation results, it can be observed that out of the 19.8 l/s flow on the drain, a maximum of 7.44 l/s flow was intercepted by the gully. The rest of the flow is overflowed through the drain outlet. The intercepted flow remains the same up to the surcharge height of 0.945 m at the manhole. Until surcharge height of 0.945 m, the flow from gully enters the manhole as a plunging jet. When the surcharge in the manhole increases more than 0.945 m, the gully starts to lose its efficiency. From 0.945 m to 1.683 m the jet from gully acts like a submerged jet. At surcharge level more than 1.683 m, the discharge is negative, i.e. the gully cannot drain to the manhole, and instead, the surcharge flow from the manhole creates reverse flow and floods the surface drain.
The gully flow can be characterised like an orifice and can be described according to the following equation: where, • Q = discharge from the gully, variable at different manhole surcharge.
• h o = Head difference from the surface drain to the gully outlet.
Here, at zone 1, h o is constant, which is equal to (h + Z − Z o =) 0.786 m. At zone 2 and 3, h o is a variable and can be calculated as the difference between (Z + h) and H.
• C d = Coefficient of discharge of the orifice and is different at the  M.N.A. Beg et al. Journal of Hydro-environment Research 19 (2018) 224-236 three different zones. For zone 1, the condition is a free orifice and C d = 0.67, as reported in different literature.
• g = Acceleration due to gravity, which is 9.81 m s −2 .
The gully discharge vs square root of head difference (h o ) has been drawn in Fig. 15 (left) in order to find the discharge coefficient of the gully. The calculated discharge coefficients are used in Fig. 15 (right) to draw the best fit curve of the discharge vs head difference relationship of the gully. The found discharge coefficients are also listed in Table 2.
It has been observed from Fig. 15 and Table 2 that the coefficients of discharge for the discussed type of gully and gully outlet are 0.677, 0.755 and 0.820 at zone 1, 2 and 3 respectively, Different water depths in the manhole also create significant changes in the shear stress at the manhole bottom. The shear stress maps of the manhole bottom at different depths are shown in Fig. 16. Depths of 0.53 m, 0.64 m and 0.93 m, represent simulations of zone 1 surcharges. The depths from 1.07 m to 1.58 m represent the conditions at zone 2. The depths more than 1.62 m represents the conditions of zone 3.
It can be seen from Fig. 16 that shear stress maps at the manhole bottom in case of surcharge heights at zone 1 have significant differences among them. The changes in these simulation results are due to the changes in surcharge heights of the manhole only as the gully discharge at zone 1 is always the same. Moreover, with the change of surcharge heights, the jet impact location in the manhole also changes. This could also be a reason for the change in the shear stress. The maximum shear stress at this zone is around 1 N m −2 . However, in the rest of the stress maps at zone 2 and 3, the changes in shear stresses are probably an effect of both changes in surcharge level as well as changed flow from the gully. At zone 2, the gully discharge decreases with the increase in water depths in the manhole. The shear stresses at this zone maintain a trend keeping higher stresses close to the manhole outlet pipe. The stress map patterns are not symmetrical at manhole water level below 1.58 m. This indicates an asymmetric circulation pattern at the bottom of the manhole resulting from the gully outlet flow. At water level of 1.58 m, the gully flow is so small that it cannot make any effect on the shear stress pattern of the manhole bottom. With high surcharge and low flow from the gully at water level of 1.58 m, the manhole shear stress was found very low and symmetrical to the main axis of manhole inlet-outlet pipe. The increase in surcharge water level reduces the maximum shear stress at zone 2 to approximately 0.3 N m −2 . At zone 3, the gully starts the reverse flow. The simulation results did not indicate very significant changes in the shear stress pattern at this zone. The shear stresses at this zone always remain asymmetrical but very low, having maximum stress in the range of 0.3 N m −2 .
The shear stress diagrams are indicative of possible sedimentation inside the manhole. As the shear stress reduces for higher surcharges, sedimentation is more likely.

Conclusions
The article presents the ability of VOF model to reproduce the flow phenomena of a manhole and gully flow at different surcharge conditions. The numerical model reproduces quasi-real scale model using open source CFD modelling tools OpenFOAM® with interFoam solver. Standard − k ε model and RNG − k ε model was used to replicate the turbulence condition at the gully and manhole respectively. The model validation was examined using velocity data at the gully as well as discharge and water depth data at the manhole.
The study results show that different surcharge level has a significant effect on gully-manhole flow. The head loss coefficient through the manhole was found very high at low surcharge condition and the flow could be explained using submerged jet theory. The flow through the gully can be comparable to an orifice flow and the surcharge level can be divided in three distinctive regions. At low surcharge, the flow acts like a plunging jet and the discharge coefficient was found as 0.677. When the surcharge level stays between the gully water surface level and gully outlet pipe soffit level, the gully discharge acts like a submerged jet to the manhole and a discharge coefficient of 0.755 was found. At a very high surcharge, the gully starts reverse flow with a discharge coefficient of 0.82.
The shear stress variation at the manhole bottom has been compared at different surcharge depths. The stress maps were found to be varying due to surcharge depth as well as the gully flow. The change in shear stress pattern showed distinctive variation in different surcharge zones and indicated that sedimentation is likely to occur at the manhole bottom for a higher surcharge.
For the future work, it will be interesting to see the effect of wall roughness on manhole head loss coefficients. The discharge coefficients for gullies with different outlet pipes may also be investigated. Lastly, the calculated discharge and head loss coefficients may be implemented in urban drainage network models. agreement no 607000.
The author would also like to acknowledge for the contribution of FCT (Portuguese Foundation for Science and Technology) through the Project UID/MAR/04292/2013 financed by MEC (Portuguese Ministry of Education and Science) and FSE (European Social Fund), under the program POCH (Human Capital Operational Programme).