Numerical investigation of driving forces in a geyser event using a dynamic multi-phase Navier–Stokes model

ABSTRACT A geyser is an explosive flow of air–water mixture shooting out of a manhole. It has been demonstrated experimentally that the releasing motion of confined air pockets in a dropshaft is the key mechanism to trigger a geyser. Other release events, unassociated with air–water mixtures, can occur, but the intensity is significantly smaller than the air–water geysers. Existing numerical models that simulate vertical air movement in mixed-phase flows typically solve a series of lumped-mass continuity, momentum and energy equations, greatly simplifying the interactions between the water and air phases. Hence, existing models are unsatisfactory in capturing the complex dynamics of a geyser because of the violent interactions between the water and air phases. In this work, a two-phase numerical model solving the Navier–Stokes Equations was applied to investigate the driving forces in an air–water geyser formation in storm sewer system. The simulated dynamics include buoyancy, air compressibility, momentum and pressure. The numerical model revealed the key factor that triggers an air–water geyser, which involves compressed air pockets that are pushed into the dropshaft by pressure surges from the main pipe. The numerical model also captured the two distinctive features of an air–water geyser, which are a violent mixture of water–air and a high-speed jet. This study also revealed how a pressure head in the main pipe, which is much lower than the ground elevation, could lift the water to the ground and push it out of the manhole.


Introduction
Geyser phenomenon has been observed in several stormwater sewer systems (Guo & Song, 1991;Hamam & McCorquodale, 1982). Geysers are a mixture of air and water shooting out of a manhole in a pressurized stormwater sewer system. Stormwater sewer systems are designed to flow under gravity within standard design return periods. However, during extreme storms, the stormwater system can experience a transition from gravity flow to pressurized flow trying to convey the additional flowrate. During the transitions, air may be trapped and compressed in the system by the pressurization front (Vasconcelos, Wright, & Roe, 2006). At ventilation points, such as a manhole or a dropshaft, the trapped air pockets escape from the system by breaking through the water layer above. If the air pockets do not break through the free surface of the water layer before they reach ground level, a water-air mixture shoots out the manhole outlet.
CONTACT Zhiyu S. Shao shaozhiyu@cqu.edu.cn In the past, it was believed that a geyser feature was directly connected to the pressure fluctuations in the main pipe, and that the water was 'lifted' out of the manhole purely by the pressure spikes in the water phase (Guo & Song, 1991;Hamam and McCorquodale, 1982), similar to pressures caused by hydraulic transients in closed conduit flow. Both experimental and numerical studies have shown that severe pressure spikes could be generated due to flow instability and water hammer effects during flow transitions and air escaping processes (Ferreri, Ciraolo, & Lo Re, 2014;Liu, Zhou, Karney, Zhang, & Ou, 2011;. However, recent field work by Wright, Vasconcelos and Roe (2011) showed that at the moment when a geyser occurred, the pressure head reading within the main pipe was much lower than the ground elevation, and hence was not enough to lift the water to the ground. Further experimental work by Lewis (2011), Vasconcelos et al. (2006), Vasconcelos, Wright, & Roe (2011), and Muller and Vasconcelos (2016) indicates that the escaping of compressed air pockets in the dropshaft is a key factor to form a geyser.
However, limited numerical models have been developed to simulate the escaping motion of a compressed air pocket in a dropshaft. Existing numerical work includes models based on one-dimensional mass and momentum equations by Lewis (2011), Vasconcelos et al. (2011) and Guo (1991). Li and McCorquodale (1999) also studied the air pocket release from a horizontal pipe into a vertical shaft. The existing models typically assume that the air pocket within the dropshaft is continuous with a distinct shape, similar to a Taylor bubble. Another typical assumption is that the air pocket occupies the whole dropshaft cross section surrounded by a steady thin layer of water. Hence the theoretical rising speeds of a static and continuous air pocket in a stagnant water column is derived based on Taylor bubble theory (Davies & Taylor, 1950). The existing numerical models generally predicted the air pocket rising motion. However, the dynamics of the air phase as well as the interactions between the water and the air phases are greatly simplified or ignored in the existing models. Hence, they are unsatisfactory in the simulation of geysers because of the violent interactions between the water and air phases. The assumption of a distinct interface between the air and water phases normally does not hold in a geyser, which is characterized by a mixture of water and air. Secondly, the lumped mass governing equations solved in these models are only applicable under the assumption that the two-phase flow is a slug flow or annular flow (Lewis, 2011). However, as indicated by Collier and Thome (1972), the slug flow regime only exists in the early stage of the air pocket evolution process. The flow turns into churn flow or annular flow regime as the air pocket pushes up. The flow structure becomes unstable, especially when the air breaks through the water surface. Thirdly, shear stress between the two phases is typically ignored in these models and it is assumed that there is no restriction for air to move upwards. In fact, as the air pocket pushes up and water flows down, the shear stress could be significant and tends to slow down the air pocket movement. This drag effect slows the air breakthrough time and causes more momentum to be transferred from the air phase to the water phase. As a result, the breakthrough point is closer to the ground, creating a greater chance for a geyser to occur.
In order to capture the complex dynamics between the water and air phases, a numerical model based on solving two-phase Navier-Stokes equations is needed. A few numerical studies based on Navier-Stokes equations have been reported to simulate the water-air two phase flows in the vertical dropshaft. Bugg and Saad (2002) solved the two-dimensional Navier-Stokes equations to study the rising motion of a Taylor bubble in a stagnant liquid in a vertical tube. The interface of the bubble is tracked using a Volume of Fluid (VOF) twophase model. However, this study was performed in a closed dropshaft system that contains a still water column instead of in a geyser where driving forces come from a main pipe. Choi, Leon, and Apte (2014) used CFD software called Star-CCM+ v7.0 and developed a 3D numerical model to study the flow characteristics in a compressed air escaping motion. This study captured the initial shock and pressure fluctuations in the main pipe after the compressed air was introduced in the main pipe. However, this model did not demonstrate the explosive behavior of a geyser. Cataño-Lopera et al. (2014) simulated a geyser reported in Chicago tunnel and reservoir plan (TARP) by coupling InfoSWMM-based hydrologic model, 1D transient model (ITMLab) and 3D FLUENT model. The hydrologic model was used to obtain the flow hydrograph in the main pipe. The transient model was used compute the transient pressure spikes due to rapid filling in the main pipe system. The 3D FLUENT model solving the two-phase Navier-Stokes equation was developed for one specific dropshaft, where a geyser was reported, to provide detailed air-water dynamics. This study revealed important dynamics of a geyser such as the process of trapping air pockets and generation of pressure spikes. This model indicated that pressure surges from transient flows in the main tunnel compress the trapped air pocket and make it act like a spring in the venting process. The compressed air pocket pushes water out of a manhole, forming a geyser. However, this numerical model did not reveal details of pressure distribution characteristics in the air pocket rising motion. It was inadequate in explaining energy transfer as the air pocket is compressed and air-water mixture is pushed toward the ground. The author suggested more detailed two-phase flow studies are needed to further understand this phenomenon. Field data analysis of a geyser in a storm-water tunnel in Minneapolis by Wright, Lewis, and Vasconcelos (2011) shows that the pressure head within the main pipe is about 20 meters below the ground, not high enough to push water out of the vertical shaft when a geyser was observed. These findings indicate there might be other driving forces in a geyser.
This work is a further expansion of the existing numerical models simulating geysers in a stormwater sewer system. Specifically, this work explores four possible geyser-triggering forces by isolating each individual force in the numerical simulations using a twodimensional numerical model based on Navier-Stokes equations. The Navier-Stokes equations solver is coupled with a VOF two-phase model to simulate the complicated air-water interactions. Whereas FLUENT treats a trapped air pocket as incompressible in a typical twophase flow simulation, this work incorporated the Ideal Gas law to reflect air compression effects inside of an air pocket. Results of this numerical study provide insight to help guide future stormwater sewer system design to alleviate possible geyser occurrences.

Methodology
The model development and verification is detailed in previous work of Shao (2013) and Shao and Yost (2013), but is highlighted here for convenience. The twodimensional Navier-Stokes equations (Equations 1 and 2) with variable density and viscosity are the basis of the model. Hence both the air and water phases are solved simultaneously. A Volume of Fluid (VOF) two-phase model is implemented to track the interface between the air and water. In the VOF method, density and viscosity are represented by a scalar variable called volume fraction, f ij , representing the volume fraction of the heavier fluid phase (i.e. water) in each computational cell. Volume fraction has the following properties, namely 0 ≤ f ij ≤ 1, with 1 representing a water cell and 0 representing an air cell. Composite density and viscosity for each computational cell are calculated from the volume fraction, as shown in Equations (3) and (4). Equation 5 is the governing equation of VOF model. This is a volume fraction transport equation and is solved at every time step to track the air-water interface. The overall governing equations of the model are where U is a velocity vector, p is pressure, ρ is density (for respective phases), B is body force and μ is viscosity (for respective phases). A turbulence model was not included in this current study. In this work, an air pocket will go through compression/expansion process as the pressure within the air pocket varies with time. It is assumed that the change of the air pocket volume is a pseudo-adiabatic process. The Ideal Gas Law (Equation 6) is used to describe the relationship of pressure and volume Here p is absolute pressure, n is amount in moles, R is the specific gas constant, T is absolute temperature, and V is volume. The temperature is assumed to remain constant, simplifying Equation 6 to Here p n+1 and V n+1 represent the absolute pressure and air pocket volume at time step n + 1, respectively, while p n and V n represent the absolute pressure and volume at time step n. The resulting pressure within the air pocket is used as a reference pressure in the solution procedure of Pressure Poisson Equation, where the location of the reference point is assumed to be at the center of the air pocket.
The Navier-Stokes equations are solved using a 1storder projection method introduced by Chorin (1968) in a staggered grid system. The solution procedure for the Navier-Stokes equations contains three parts: 1. solve the momentum equation in the absence of pressure to obtain an intermittent velocity field; 2. solve a Pressure Poisson Equation (PPE) using the intermittent velocity; 3. project the intermittent velocity field to a divergencefree velocity field. The discretization is performed using a finite volume method. A quadratic upstream interpolation for convective kinematics (QUICK) scheme is used to discretize the advective terms in the momentum equations. In cases where a reference pressure is needed for the air pocket, it is calculated before solving the Pressure Poisson Equation.
A Shuman Filter is applied to the auxiliary velocity field after the momentum equations are solved to handle the 'Cell-Re' and aliasing problems. This filter convolutes an irregular solution to a smooth one by providing additional dissipation needed to dampen aliasing. The solution procedure for the Navier-Stokes equations, VOF equation and Shuman Filter are detailed in other works by Shao (2013) and Shao and Yost (2013). The model presented in this work was extensively tested and validated (Shao, 2013). These validations included fundamental fluids applications (i.e. liddriven flow, Poiseuille and Couette flows, etc.), multiphase flow applications (i.e. rigid body rotation, shearing flow, Rayleigh-Taylor instability, sloshing tank, dam break, etc.), transient flow applications (i.e. free surface hydraulic bores, pressure surge waves, etc.). Associated numerical stability and grid convergence tests have been performed as part of the model development effort. In all cases the results compared extremely well with the analytical/experimental results.

Dam break problem
Dam-break simulation is a benchmark problem for twophase flow models due to the simple initial and boundary conditions. Dam break problems involve significant interface deformation such as overturning, breaking up and air entrapment. The interface dynamics are a challenge to some two-phase models. The experimental work of Martin and Moyce (1952) is commonly used to check the numerical results. In the dam break problem, initially a rectangular column of still water is contained between a vertical wall and a gate. At time t = 0+, the gate is suddenly removed and the water column starts collapsing under gravity. This collapsed water column forms an advancing water wave, propagating to the right. The numerical computational domain size is the same as Martin and Moyce's experiments (1952), which is a square of 4a × 4a, where a is the initial water column width. The initial water column is 0.05175 m (a) wide and 0.1035 m (2a) tall. The densities of water and air used in the simulation are 1000 kg/m 3 and 1.23 kg/m 3 , respectively. Viscosities for water and air are 1.0 × 10 −3 kg/m·s and 1.8 × 10 −5 kg/m·s, respectively. The gravitational acceleration is 9.8 m/s 2 . The grid size is 64 by 64 with a time step size of 5 × 10 −5 sec. A non-slip boundary condition is used for the two vertical walls and the bottom boundary. The top boundary is open.
To compare with experimental results, the wave front location (x), remaining water column height (h) and time (t) are converted to the non-dimensional numbers published in Martin's work. This is done by normalizing wave front location (x) using x * = x/a, and the corresponding time (t) using T * = t ((2g)/(a)). Figure 1(a) shows the surge front locations and compares the numerical predictions to Martin's experimental results. Similarly, the water column height (h) is normalized by h * = h/a and the corresponding time (t) by t * = t (g/a). Figure  1(b) shows the water column height with time compared with the experimental measurements. It can be seen from Figure 1 that the numerical prediction of wave front location and remaining water column heights agree very well with the experimental data.
Snapshots of the interface from the experiments by Koshizuka, Tamako and Oka (1995) and the numerical prediction at selected times are displayed in Figure 2. It can be seen that numerical predicts the whole complicated interface deformation process very well. This case demonstrated that the numerical model is capable of simulating complex two phase flow with air entrainment, interface breaking up, merging and overturning.

Uniform open channel film flow
The second validation simulation is a two-dimensional open channel film flow with a known analytical solution. In this case, a constant thickness film of viscous liquid flows down a plane that is inclined at angle θ , as shown in Figure 3. The flow, driven by gravity, is a uniform flow, frequently found in a stormwater systems. In the uniform channel flow, gravity acceleration is balanced by frictions from channel bottom. As a result, flow within the channel maintains the same profile along the flow direction.
The analytical solution of this type of flow could be derived from two-dimensional Navier-Stokes equations by assuming a fully developed flow along the flow direction (all variable gradients are zero). Additionally, the presence of a free surface indicates a zero flux across the interface. The shear stress at the top of the flow (water   surface) is zero. The analytical velocity profile can be derived theoretically as (Chaudhry, 1993).
where H is the total flow depth in channel and y is the flow depth. The computation domain includes a main pipe and a manhole as shown in Figure 3. The main pipe is 8 m long and 0.4 m tall. The manhole domain is 0.2 m wide and 0.1 m tall. The manhole center is located at x = 5.05 m. The heavier fluid has a viscosity value of 0.05 kg/m·s and a density value of 1000.0 kg/m 3 . The lighter fluid has a viscosity of 0.001 kg/m·s and a density of 1.0 kg/m 3 . Body force along flow direction is B x = g*sin(θ ), and normal to the flow direction is B y = −g*cos(θ ) with g being gravity acceleration. Initially, the flow field is uniform with an analytical profile. The initial flowrate for the heavier fluid is 0.01 m 3 /s and the initial flow depth is 0.05 m.
A non-slip wall boundary is used for channel bottom boundary. A free-slip wall boundary is used for the top of the channel to simulate an atmosphere open field. A specified velocity profile inflow boundary is used for the inlet. The inflow velocity profile is a uniform velocity obtained by dividing the flowrate by depth. The flowrates are 0.01 m 3 /s for the heavier fluid and 0.007 m 3 /s for the lighter fluid. At the outlet, an outflow boundary condition is implemented. Since gravity is included in the simulation, pressure at the outlet has a hydrostatic distribution. It can be verified that Froude number for this problem is Fr = 0.7 . The downstream boundary depth is the normal depth corresponding to the channel slope and flowrate, which is 0.05 m. The time step size is set to be δt = 0.001 s.
To observe the grid convergence trend, three sets of grid size are used in the simulation; namely, 80 × 20, 160 × 40 and 320 × 80. Figure 3 plots the velocity profile at the outlet predicted by the model using three grid sizes, and compares them with the analytical profile. The numerical solution matches the analytical result very well, and as the grid size is refined, the numerical solution converges to the analytical solution.

Numerical simulations
This work explores the four air-water geyser driving forces by strategic application of the bottom boundary conditions, representing conditions in the main pipe and associated impact to the air pocket rising motion in a vertical riser. The four possible driving forces include: buoyancy, air compressibility, momentum, and pressure. The riser is 0.1 m wide and 4.0 m tall. The computational domain is shown in Figure 4, as well as showing the initial set up of each simulation, strategic boundary conditions and the important observations (i.e. free surface displacement) during the simulations Table 1 shows a summary of the results. The densities of water and air in the simulation are 1000 kg/m 3 and 1.23 kg/m 3 , respectively. Viscosities of water and air are 1.0 × 10 −3 kg/m·s and 1.8 × 10 −5 kg/m·s, respectively. Surface tension was not included; Based on Weber Number values, surface tension has an extremely minor impact on the dynamics during the whole simulation period due to the fact that the simulation stops as soon as the surface breaks.
The numerical simulations are carried out as follows. Initially, an air pocket at rest is contained below a stagnant water column and is separated by a hypothetical film, as shown in Figure 4. The air phase has a height of H a and the water phase has a height of H w. The shaded area closed by bold solid lines in Figure 4 represents the initial water column. The same initial condition is used for all simulation cases. At t = 0 + , the hypothetical film is removed and the air pocket starts rising in the vertical riser due to buoyancy and other driving forces defined at the bottom boundary. The bottom boundary is a key parameter in all simulations as it represents the individual driving force that is fed from the conditions in the main pipe system. The bottom boundary is defined in a way so that different driven forces could be isolated and studied. Depending on the simulation case, the bottom boundary is defined as a non-slip wall if there is no inflow coming in the dropshaft or an inflow boundary if there are external driving forces coming from the main pipe. The top boundary is defined as an open boundary with a constant prescribed pressure (i.e. atmosphere pressure). The pressure values shown in all figures are gauge pressures (relative pressure). The walls of the riser are non-slip.
The dashed line in Figure 4 represents the 'final' free surface at the moment when the air pocket breaks through the free surface. Since the final water column is a mixture of water and air with an irregular shape, the final free surface height refers to the maximum y-coordinate value at the moment when the air pocket breaks through the free surface and escapes. The distance between the initial and final free surface is defined as 'free surface displacement', which represents the height of water lifted by the corresponding driving force.
It should be pointed out that there are two distinctive interfaces between the water and air before the air pocket breaks through. The first one is at the top of the water column that contacts with the atmosphere; this interface is referred to as the free surface. The second interface is between the water column and the air pocket within. This interface is referred to as the air pocket front or air pocket nose. The upward speed of the free surface or water-air mixture can be calculated based on the displacement of the free surface (y 1 and y 2 ) in any time interval (t 1 and t 2 ) as v fs = y 2 − y 1 t 2 − t 1 The final free surface height is usually called geyser strength (Vasconcelos et al., 2006). It is normalized by the initial pressure head, P o , to characterize the effects from each individual driving force: y * fs = y fs P 0 (10)

Case 1: Buoyancy effect
The first simulation investigates an air pocket rising motion induced by buoyancy effects. In this case, the initial pressure (P 0 ) within the air pocket is set equal to the hydrostatic pressure of the water column, namely, P 0 = ρgH w , where, H w is the height of the water column. The bottom boundary is a non-slip solid wall. Hence the air pocket is in a close system without external forces, isolating this from other driven forces. This would represent a unique situation in the main pipe where air just happens to be trapped below water and it trying to escape due to density differences (buoyancy). In this simulation case, the initial height of the air pocket is 0.5 m. The initial height of the water column is 1.0 m, which results in an initial pressure head of 1.0 m within the air pocket. The computational grid size is 20 by 160 and the time step size is 5 × 10 −5 s. Table 1 lists the numerical parameters used in this simulation.
The numerical results in Table 1 shows that the speed of the free surface reached 0.4 m/s at the time when the air pocket broke through the surface. This speed is consistent with the results of a Taylor bubble. According to Davies and Taylor (1950), the speed of a well-rounded nose Taylor bubble in a vertical cylindrical tube has a velocity of where, g is gravity and D is the diameter of the tube. In this configuration, the speed of a Taylor bubble is calculated to be 0.35 m/s. The value predicted by the numerical model is slightly higher due to the acceleration of the air pocket near the free surface. This free surface speed is fairly low comparing to a typical water air mixture jet speed in a geyser, which is generally in the order of 10 ∼ 20 m/s (Cataño-Lopera et al., 2014;Wright et al., 2011). However, as the simulations are based on a much smaller dropshaft diameter (0.1 m versus 2.5 m in the cited field observations), scaling will have some effect. Given that there are two different fluids, and hence gravity plays an important role, Froude number similarity governs the comparison between scales. Here the simple Froude number is used (versus the densimetric Froude number) because the density difference between air and water is nearly equivalent to the density of water. Thus, in this scale comparison, the velocity in this modeled dropshaft would scale to the velocity in the field by 1/5 ( √ 0.1 m/2.5 m). Based on this similitude, the velocity results would be compared to 2 ∼ 4 m/s, 1/5 of the literature reported 10 ∼ 20 m/s. The gap between the simulated value (0.4 m/s) and reported free surface velocity (scaled to 2 ∼ 4 m/s) indicates that the momentum generated purely from buoyancy is minor. Buoyancy effect, which could enhance geyser strength, is clearly not a dominant force in geyser formation.

Case 2: Air compressibility effect
The second study case involves a compressed air pocket. This case is designed to investigate a rising motion driven by air compressibility after isolating this from other forces. The air pocket is initially compressed by a pressure that is higher than the hydrostatic pressure. In this simulation, a pressure that is twice the hydrostatic pressure is initially set within the air pocket. This case is not a typical situation found in a real system, as this would be consistent with an air pocket of very limited size. But this case is presented to specifically isolate what compressibility does to geyser formation. The pressurization of a real system is limited so this case presents a reasonable value for a pressure surcharge. The bottom boundary is a non-slip solid wall. Hence the air pocket rises and expands in a closed system. Initial conditions and other numerical parameters are listed in Table 1.
It is observed from Table 1 that when the air pocket breaks through the surface, the free surface displacement is insignificant. The free surface is lifted slightly due to air pocket expansion. The speed of the air pocket front reaches 1.1 m/s when the air pocket breaks through the surface. It is slightly higher than the speed of a Taylor bubble in Case 1. This is expected because a compressed air pocket contains more energy that could be transferred to the water column. However, this speed is still significantly lower as compared to reported geyser jet speed, even after scaling. This indicates that the momentum generated by the air compressibility effect is minor and not enough to trigger a geyser.

Case 3: Momentum effect
In this simulation, the air rising motion is driven by momentum (i.e. velocity) that is introduced by inflow into the vertical riser at the bottom boundary. This simulation is designed to mimic the momentum gained from the main stormwater sewer pipe. In a stormwater sewer system, momentum could be transported from the horizontal main pipe to the vertical dropshaft when a water column or an air pocket is pushed into the vertical shaft. The velocity in the main pipe could originate from a propagation surge front during a pressurization from rapid filling, or pressure fluctuations due to water hammer effects. Early researchers believed that the momentum gained from the main pipe 'lifts' the water and air in the dropshaft and forms a geyser. Therefore, this driven force is isolated and studied in detail to investigate the strength and impacts.
In this simulation, the air pocket is initially contained beneath a water column under hydrostatic pressure. At time t = 0 + , a velocity inflow boundary with a uniform profile is applied at the bottom boundary. This setup is similar to a plug flow where momentum (or velocity) is continuously fed by a large air pocket that is pushed and propagating toward a manhole. In this simulation, it was assumed that the air pocket is simply pushed into the dropshaft with a typical velocity, but it is not compressed. Table 1 lists the initial conditions and numerical results. Figure 5 plots the interface evolution at various times. It can be seen from Table 1 that when the air pocket breaks through, the free surface is pushed up to an elevation (2.7 m) that is much higher than the previous two driven forces (Cases 1 and 2). In this case, the free surface rises significantly high with a moderate speed. It is shown in Table 1 that the speed of the interface reaches 2.1 m/s at the time of breakthrough (1.7 m displacement), which is higher comparing with the surface lifted by buoyancy (Case 1) and compressibility (Case 2). Under Froude number scaling conditions, this velocity starts to reach the low end of field observed velocity ranges. Figure 6 plots the interface displacement of the air pocket rise at different times. It is observed that the interface rises with an almost same speed as the inflow. This result is expected because, based on plug-flow theory, the air pushes the water column and all move together with the same speed. The results show that the momentum gained from the main pipe could push water to a fairly high elevation with a moderate speed. Hence momentum gained from the main pipe has an important influence on the geyser strength, but may not be the key force.

Case 4: Pressure effects
In this case, the characteristic of air rising motion driven by pressure effects is investigated by prescribing a pressure inflow boundary at the bottom boundary. This configuration represents a scenario where a manhole receives a compressed air inflow from the main pipe. This scenario typically occurs at times when a portion of a large air pocket is compressed and pushed into the manhole while the rest of the air pocket, namely, the air pocket tail, is still contained in the main pipe. While the pressurevolume relationship can be captured by the numerical model, the implementation of the bottom boundary condition creates a very large (infinite sized) air pocket in main channel. The pressure in the main pipe could originate from a pressurization surge front or pressure spikes due to water hammer effects. The distinction between this case and Case 3 is that the air pocket is now pushed and compressed whereas before the air pocket is simply pushed, but not compressed. The difference between Case 2 and Case 4 is the extent of the pressurized air pocket. In Case 2, the air pocket is small and hence is fully contained in the vertical shaft (as given by the closed boundary at the bottom). Case 2 simulation details how a small pressurized pocket (from the conditions in the main channel) impacts the geyser when it moves into the shaft. In this case the boundary condition applied at the bottom of the shaft creates an infinite sized air pocket. The model uses the ideal gas law to determine the pressure conditions in the air pocket, hence with the pressure being fixed at the boundary in this case, it is consistent with a very large air pocket.
In this simulation, a continuous pressure inflow with a fixed value is specified at the bottom boundary. In a real stormwater sewer system, the pressure value at the bottom of the dropshaft varies with the pressure head in the main pipe. In this numerical setup, the main pipe is pressurized and a pressure head higher than the hydrostatic pressure is specified at the inlet boundary. This is similar to a manhole that is surcharged while a compressed air pocket is pushed into the dropshaft by a pressure surge. To investigate the impacts caused by different levels of compression, two initial pressure head values are simulated in this case. In Case 4.a, an initial pressure head of 1.0 m is used while in Case 4.b a value of 2.0 m is used.
The instantaneous interfaces at different times are plotted in Figures 7 and 8. From the interface evolution history show in these two figures, a thin layer of film-flow around the rising air pocket can be observed. At the time when the free surface breaks up, the air flow mixes with the falling water column and forms a violent water-air mixture. Table 1 lists the numerical results of both Case 4.a and 4.b. One feature illustrated in Table 1 is the extremely high velocity near the air pocket front. The speed of the interface reaches 10 m/s at the time of breakthrough in Case 4.a and 18.5 m/s in Case 4.b. Comparing to Case 3, this simulation result indicates that the force generated by an air pocket that is compressed and pushed into the dropshaft is more significant than an air pocket that is simply pushed into the dropshaft. The numerical result is more consistent with the unscaled velocities extrapolated in the field as reported in the literature  although it exceeds the velocities observed in previous experiments (Lewis, 2011). However, the numerical result demonstrates an important factor that the compressed air pockets released from a surcharged main pipe are a much more dominant factor in a geyser formation. It also revealed that the velocity is highly dependent on the degree of pressurization because Case 4.b was double that of Case 4.a. This can also be observed from the comparison between Case 4 and the Case 2, as both dealt with compressed air. However, Case 2 did not have nearly the impact on the velocity as Case 4. In Case 2, the simulated pressurized air pocket is relatively small and the pressure is transferred rather quickly (using ideal gas relationship), but very limited (size of the pocket), in driving the geyser. In a real air-water geyser system, as the air expands, the pressure drops. Hence a smaller volume of compressed air cannot sustain  a geyser event. The setup of Case 4 is consistent with a very large amount of air in the main tunnel moving into the dropshaft. Hence the source and movement of compressed air appears to the key feature in geyser strength. Another feature demonstrated by Figure 8 is the violent mixture of air-water phases and flow chaos that distinguish a geyser from a typical stormwater overflow. Hence, it is reasonable to conclude that the two distinct features of a geyser event have been captured by this numerical model.
Comparison of Case 4.a and Case 4.b also indicates that higher pressure within the pocket leads to a jet with higher velocity. Higher pressure means the air pocket is further compressed by the pressure spikes in the main pipe. This result demonstrates a linkage between the surcharged state in the main pipe and the geyser strength. The free surface displacement and free surface velocity of Cases 4.a and 4.b at different time are plotted in Figure 9. From Figure 9, a noticeable acceleration could be observed as the air pocket moves closer to the free surface. These results agree with previous experimental results (Muller and Vasconcelos, 2016;Vasconcelos, Wright and Roe 2011) in which they showed a similar acceleration process as the air pockets approach the free surface.
A surprising result is observed in Table 1 when comparing the final free surface height with the pressure head in the main pipe, which is represented by the pressure head at the bottom boundary. In both Cases 4.a and 4.b, water is lifted to an elevation much higher than the actual pressure head prescribed at the boundary. In Case 4.a, the pressure head at the inlet boundary is 1.0 m while the final free surface reaches a height of 1.8 m, which is nearly twice of the pressure head in the main pipe. Similar pattern is observed in Case 4.b. In Case 4.b, the pressure head is 2.0 m while the water surface is lifted to an elevation of 3.1 m. In fact, this result parallels what was observed in the field study by Wright et al. (2011), in which the pressure head in the main pipe was only 6.0 m but with the compressibility of an air pocket it could push a water-air mixture to an elevation of at least 28.6 m, which is the distance between the tunnel invert and the ground surface.
This phenomenon could be explained by the significant pressure gradient near the air pocket front, as shown in Figure 10. When the air pocket rises and approaches the free surface, the pressure near the air pocket front rapidly drops to atmosphere pressure. This generates a significant pressure gradient between the nose and the tail of the air pocket. The pressure gradient results in high  flow acceleration near the air pocket front. If the breakthrough point is close to the ground level, the water-air mixture will shoot out of the manhole and form a geyser.

Conclusions
A two-phase numerical model solving the Navier-Stokes equations was applied to investigate the driving forces in air-water geyser formation in stormwater sewer systems. The numerical model was applied to simulate the violent two-phase water-air mixture flow during an air pocket rising and releasing process in a vertical dropshaft. The numerical simulations explored the primary geyser driving forces including buoyancy, air compressibility, momentum and pressure.
The numerical model reveals the key factor that triggers a geyser feature, namely that when a compressed air pocket is pushed into a vertical dropshaft and escapes from the system, a violent water-air mixture with a high velocity forms, and hence could trigger a geyser. The results also indicate that just having captured air in the system is not enough to produce the air-water geyser. Simple buoyancy effects, or small pockets of compressed air have minimal impact on the geyser formation. Furthermore, momentum gained from the movement of air or water in main pipe could promote the strength of a geyser but is not the main triggering force. The results show clearly that the occurrence and strength of an air-water geyser depends on a source and movement of a compressed air pocket. Hence if a system that is prone to strong air-water geysering, removing the large source, movement and/or the pressurization of the air pocket seems to be the key to reducing/eliminating the air-water geyser event. The two-phase numerical model captures the two distinctive features of an air-water geyser, (violent mixture of water-air and high-speed jet) reasonably well.
Another important contribution of this study is that it addresses the existing knowledge gap, which is explaining how a pressure head much lower than the ground level in the main pipe could potentially lift the water to the ground and push it out of the manhole.

Notation:
The following symbols are used in this paper:

D
Diameter (

Disclosure statement
No potential conflict of interest was reported by the authors.