## Abstract

Storm geysers have received significant attention lately due to its more frequent occurrences and the induced severe local flooding and infrastructure damages. Previous studies suggested that the air pocket pressure oscillated during geyser events especially in rapid filling process, but only the peak values were studied and the oscillation period was not discussed in detail. In this paper, a theoretical model was developed focusing on the period of the pressure oscillation induced by the expansion/compression of the air pocket below a water column in a vertical riser with film flow. It was found that the oscillation period was a function of the initial air pocket volume, initial air pocket pressure head, the riser diameter, and the initial water column length. The oscillation period increased with the air pocket pressure head and the air pocket volume, but decreased with the riser diameter and the polytropic coefficient. The oscillation period increased then decreased with an increasing water column length. Further, when considering the film flow along the riser, the oscillation period decreased slightly from the analytical solution. It was also found that the inflow rate change did not significantly influence the oscillation period.

## HIGHLIGHTS

Theoretical model was developed on the period of the pressure oscillation induced by the air pocket below a water column in a vertical riser with film flow.

The oscillation period increased with air pocket pressure head and the air pocket volume, but decreased with the riser diameter.

When considering the film flow along the riser, the oscillation period decreased slightly from the analytical solution.

## INTRODUCTION

Commonly seen as air/water mixture erupting out of manholes, storm geysers appear more and more frequently due to a variety of reasons such as climate change, urban growth, system aging, among others. Storm geysers can negatively impact on municipal infrastructures and public safety such as the loss of manhole cover, local flood, or infrastructure damage, etc. (Li & McCorquodale 1999; Shao 2013; Wright *et al.* 2011). Some recent studies reported that the pressure in sewer systems showed an oscillation pattern during geyser events (Liu *et al.* 2020). But these studies mainly focused on the magnitude of the pressure without much discussion on the oscillation period. The pressure oscillation pattern is important for further analysis such as the forces on municipal infrastructures, vibrations, system stability, geyser prevention, system protection, among others. Therefore, establishing an explicit relationship between the oscillation period and the flow conditions is of importance for developing a feasible storm geyser prevention and mitigation plan.

Researches have been conducted on the mechanics of storm geyser events. In general, there are two main streams that lead to geyser events. One is the rapid filling of sewer pipes, which was studied in Zhou *et al.* (2002) where it was found that with orifice plates installed at the pipe end, the impinging pressure can reach up to 15 times of the driving pressure. This phenomenon was also observed in Zhou *et al.* (2011), Li & Zhu (2018), and Qian & Zhu (2020). The other mechanic is the release of air pocket trapped in the storm system, which was investigated in Vasconcelos & Wright (2011), and Cong *et al.* (2017) etc. Liu *et al.* (2020) experimentally studied the generation of geyser events using a physical model. The study explored the factors that triggered the geyser events and proposed a preliminary model for predicting the period and magnitude of the pressure oscillation for the free oscillation scenario. However, the model did not consider the compression/expansion of air pocket during the geyser process. Qian *et al.* (2020) numerically simulated storm geyser events and proposed potential geyser mitigation methods. For the above studies, no explicit relationship was proposed between the period of pressure oscillation and other factors.

For theoretical models, a combined rigid column method (RCM) and method of characteristic (MOC) model was proposed by Zhou *et al.* (2002, 2004). The model was then further developed by Huang *et al.* (2018a, 2018b) to account for the vertical motion of water. Vasconcelos & Wright (2011) developed a theoretical model on predicting the movement of a water column in a riser driven by an air pocket. Qian & Zhu (2020) proposed a theoretical model predicting the pressure surge when using orifice plate on mitigating the geyser events. The models mentioned above were all solved numerically and lack of analysis on pressure oscillations. Huang & Zhu (2020) proposed an analytical solution on the pressurization of rapid filling in a horizontal pipe with entrapped air. For the movement of a water column in a riser driven by an air pocket beneath it, there is no analytical solution yet.

This study was conducted to develop an analytical solution for explicitly show the relationship between the pressure oscillation period and other parameters such as the dimension of the structure, flow conditions, etc. After a comprehensive evaluation, the model explicitly showed the pressure oscillation period. The behavior of the proposed model was then assessed and discussed in detail. The cases considering the mass loss of water column due to film flow as well as the flow rate change were solved numerically and the impact of the film flow and flow rate change on the pressure oscillation period was discussed.

## METHODOLOGIES

### Governing equations

The control volume for the model is defined in Figure 1. The model consists of a riser with a diameter of *D*_{r}. The air volume is confined with the flow system with volume of ∀_{a}, and absolute pressure head of *H**. The tank in Figure 1 is for illustration purpose. The water column in the riser has a length of *h*, and the position of the top and bottom of the water column is *Y*_{u}, *Y*_{d}, with velocity of *v*_{u}, *v*_{d}, respectively. The length of the water column decreases due to the film flow along the riser wall at a thickness of *δ* and flowrate of *Q*_{f}. In prototype systems, the difference in inflow and outflow rates may cause extra compression or expansion on the air pocket. To include this flow rate change scenario, the bottom of the chamber is movable and the air pocket volume change rate is *Q*_{b} which equals to the difference between the water flow into and out of the system. During the process, assuming the pressure head of the air pocket is higher than the weight of the water column, the water column is firstly driven by the air pocket and moves upward along the riser. As the water column moves up, the air pocket expands and the pressure head decreases until the point where the weight of the water column balances with the pressure head in the air pocket. The water column then moves downward and compresses the air pocket resulting in a pressure rise in it and thus forms a cycle.

*k*). The effect of the polytropic coefficient is going to be discussed later in this paper. It is also assumed that the height of the riser is sufficiently large to contain the water column during the entire oscillation process. Additionally, the film flow is the water flow from the water column along the riser wall causing mass loss of the water column, and the length of the water column decreases due to the film flow. It is then assumed that the initial film flow is modelled as Taylor Bubble and the expression is:where: and

*δ*can be written as: with

*R*= 0.99, and

^{2}*g*is the gravity acceleration (Qian & Zhu 2020). The water column length can be expressed as:

The governing equations above are fairly straight forward to follow where Equation (3) shows that the expansion or compression of air pocket induces the upward or downward movement of the water column. Equation (4) suggests that the acceleration of the water column is a result of the pressure difference between the driving pressure and the gravity of the water column. Equation (5) is the classic ideal gas equation which correlates the pressure and volume change of the air pocket. The polytropic coefficient *k* in Equation (5) represents different process for ideal gas, where *k* = 0 is isobaric process, and *k* = 1 is isothermal process.

### Geyser process without film flow

If *H*_{0}^{*} > *h*_{0} + *H*_{atm}, the initial pressure in the air pocket is higher than the weight of the water column. It is to push the water column upward, and oscillates above the initial location. If *H*_{0}^{*} < *h*_{0} + *H*_{atm}, the initial pressure in the air pocket is less than the weight of the water column and it needs to be compressed to increase the pressure, and the equilibrium point is below the initial location. It should be notated that if the difference between *H*_{0}^{*} and *h* + *H*_{atm} is high enough, the equilibrium point can be far from the initial position and the linearized model may not reflect the actual movement. The detailed limitation of the proposed model is discussed later in this paper.

The extreme values are and . To avoid complex values, should be larger than . This means that for , *H*_{0}^{*} < *h* + *H*_{atm}, and the volume of air pocket at the equilibrium state is less than the initial air pocket volume, and for , the initial air pocket pressure is larger than the water column length and the water column is pushed upward in the beginning. After some algebra, the validated region of the dimensionless head can be expressed as . For the term in the right-hand side, it is in the order of 0.01 which is significantly less than 0.5 and therefore can be dropped. Previous study suggested that *k* = 1.4 (Qian & Zhu 2020) and, therefore, the mode limit for *k* = 1.4 is .

### Geyser process with film flow

The additional governing equations along with the previous ones cannot be explicitly solved as analytical solution, and numerical approach is applied for the model. The governing equations were solved numerically using a fourth order Runge-Kutta algorithm with Matlab.

Table 1 shows the cases tested in this study. Runs A were the cases solved by the analytical solution and mainly for evaluating the model with Vasconcelos & Wright (2011) and Liu *et al.* (2020). Run B1 was compared with Vasconcelos & Wright (2011) for model evaluation. Run B2 was for testing the effect of initial water column length under a variety of boundary conditions (, *D _{r}*, and

*H*). Seven values were tested for each variable which results in a total of 2,401 runs for Run B2. Run C1 was for model evaluation with Liu

_{0}*et al.*(2020). Run C2 was mainly for testing the effect of the inflow rate change and a total of 22 cases were tested where the bottom air volume change rate ranged from

*Q*

_{b}= −1,000 L/s to 1,000 L/s.

Run . | Solution type . | h (m)
. _{0} | (L) . | D (m)
. _{r} | H (m)
. _{0} | Q (L/s)
. _{b} |
---|---|---|---|---|---|---|

A1* | Analytical | 0.254 | 3.79 | 0.057 | 0.305, 0.61, 0.915 | N/A |

A2* | 0.356 | |||||

A3* | 0.457 | |||||

A4* | 0.300 | 9.61 | 0.06 | 2.93 | ||

B1* | Numerical | 0.25 | 3.79 | 0.057 | 0.264 | |

B2 | 0.1–10.0 | 10–5,000 | 0.02–2.0 | 0.1–20.0 | ||

C1* | 0.450 | 9.61 | 0.06 | 0 | 62.5 | |

C2 | 5.0 | 2,000 | 1.2 | 5.0 | −1,000–1,000 |

Run . | Solution type . | h (m)
. _{0} | (L) . | D (m)
. _{r} | H (m)
. _{0} | Q (L/s)
. _{b} |
---|---|---|---|---|---|---|

A1* | Analytical | 0.254 | 3.79 | 0.057 | 0.305, 0.61, 0.915 | N/A |

A2* | 0.356 | |||||

A3* | 0.457 | |||||

A4* | 0.300 | 9.61 | 0.06 | 2.93 | ||

B1* | Numerical | 0.25 | 3.79 | 0.057 | 0.264 | |

B2 | 0.1–10.0 | 10–5,000 | 0.02–2.0 | 0.1–20.0 | ||

C1* | 0.450 | 9.61 | 0.06 | 0 | 62.5 | |

C2 | 5.0 | 2,000 | 1.2 | 5.0 | −1,000–1,000 |

Asterisk signs represent model evaluation cases.

### Model evaluation

The analytical solution of the model was compared with Vasconcelos & Wright (2011) for Runs A1–A3, and Liu *et al.* (2020b) for Run A4 in terms of the period of the pressure oscillation. The results for Runs A1–A3 are shown in Figure 2. It is notable that in Vasconcelos & Wright (2011), the pressure oscillation after the air pocket reached the vertical riser was used. The figure suggested a good agreement between the present model and the experimental data where the differences were within 10%. The discrepancies mainly came from the mass loss due to the film flow in the experiment. For Liu *et al.* (2020), *D _{r}* = 0.06 m, = 9.61 × 10

^{−3}m

^{3},

*h*≈ 0.30 m (approximated from the figures when the pressure reaches the maximum), and

*H*

_{0}= 2.93 m (maximum pressure in the air pocket). The calculated period of the pressure oscillation was 0.59 s while the measured period in Liu

*et al.*(2020) was 0.53 s. With a difference of about 10%, the comparison is reasonable.

For the model developed by Huang & Zhu (2020), assuming the driving pressure head equals the water column length in the present model, and the ratio of the present analytical solution and the result of Huang & Zhu (2020) yields where *x _{e}* is the displacement of air/water interface in the horizontal pipe at the equilibrium point relevant to its initial location. The ratio equals 1 when the initial condition of Huang & Zhu (2020) is at its equilibrium point, which means the two models give identical results. For

*x*> 0 (i.e. driving pressure in water tank is higher than the air pocket pressure), the ratio is less than 1 which means the model of Huang & Zhu (2020) gives a larger oscillation period than the present study and vice versa. In Huang & Zhu (2020), the air pocket was trapped in a closed pipe and the water column was driven by a constant pressure head, which may not fully represent the prototype manhole connections. The current study is a reproduction of the physical models that generated storm geyser events and it is believed that it better represents the real storm geyser case.

_{e}For Run B1, the numerical solution considered the loss of the water mass and its solution was compared with Vasconcelos & Wright (2011) in Figure 3(a). The initial water column length and air pocket pressure were determined based on the flow condition immediately when the air pocket entered the riser in the experiment instead of the global initial condition of the experiment. The numerically solved model compared well with the measurement data. It was noticed that the right-hand side of Equation (2) should be multiplied by 1.32 to match the physical measurement. This was mainly because the pressure in the air pocket drove the water column upward and enhanced the film flow causing an increased film flow. The top of the measured and simulated water column varied in an oscillation manner with a slight upward trend. The bottom of the water column moved upward over time and made the water column shorter. With respect to the pressure head in the air pocket, the measured and calculated pressure generally decreased over time due to the expansion of the air pocket pushing the water column upward. The oscillation period and its magnitude could also be simulated using the numerical solution with less discrepancies. Therefore, the numerical solution can be used for further analysis.

For numerically solved model considering the flow rate change, the results of Run C1 is compared with Liu *et al.* (2020) in Figure 3(b). The figure suggests that the pressure in the air pocket can be simulated generally well, especially for the first pressure oscillation cycle. The calculated peak pressure was 3.22 m and the measured maximum *H* was 2.93 m. The difference was about 10% and the comparison was reasonably good. For the numerical solution after the first peak pressure, discrepancy occurred. This was mainly due to the assumption of the current model that the riser is infinitely long to contain the water column. In Liu *et al.* (2020), the riser was 1.22 m long and after the first pressure peak, the water column exits the riser and caused the discrepancy.

The Taylor bubble theory was applied only for estimating the flow rate of the film flow (Equation (1)) which was assumed as a constant for a given riser, and the length of the water column in the riser decreased due to the mass loss induced by the film flow (Vasconcelos & Wright 2011; Qian & Zhu 2020). For Runs A, when deriving the solution, the film flow rate was contained in the equations. Nevertheless, after the linearization, the film flow term in the analytical solution for the pressure oscillation period (Equation (10)) was cancelled out. Therefore, the film flow theory did not affect the analytical solution. For Runs B and C, the film flow was considered. With a generally good comparison between the models and experimentally measured water column change, the Taylor bubble assumption is believed to be reasonable for the study.

## RESULTS AND DISCUSSION

### Parametric analysis for analytical solutions

Equation (10) suggests that the pressure oscillation period increases with the initial air pocket volume () in the power of 0.5, and with the initial air pocket pressure head (*H _{0}**) in the power of , and decreases with the cross-sectional area of the riser (

*A*) in the power of −0.5. There is no explicit trend between

_{r}*T*and

*h*from Equation (10). The derivative with respect to

*h*of Equation (10) suggests that

*T*increases with

*h*before

*h*=

*kH*

_{atm}. Then,

*T*decreases with the increase of

*h*. When

*k*= 1.4, the

*h*corresponding to the maximum

*T*is when

*h*= 14.42 m.

The analytical model was developed with several assumptions and simplifications. Therefore, uncertainties and limitations existed. In the current study, because the analytical solution was only valid near the equilibrium point, the variation of the water column length was ignored. In the experiment of Vasconcelos & Wright (2011), and Liu *et al.* (2020) and the analysis of Qian & Zhu (2020), the length of the water column decreased due to the film flow along the riser wall. However, considering the decrease of the water column length, it may be impossible to have the analytical solution on the oscillation period on the pressure head. Therefore, the pressure oscillation period may be less due to the decrease of the water column if the water column is less than *kH _{atm}*.

Additionally, the model may have potential limitations in terms of the initial conditions. The model only works when the normalized pressure head *H ^{’}* is larger than 0.38 (when

*k*= 1.4). Nevertheless, this is fairly wide range considering the atmospheric pressure head of 10.3 m. For example, to ensure that the water column length (

*h*) is larger than 0 m, the minimum driving pressure head (

*H*

_{0}) can reach as low as −6 m. This is already a vacuum condition where the water column is initially forced to move downward along the riser. For a driving head of atmospheric pressure, the maximum

*h*can reach up to 16.8 m. In this case, the water column initially moves downward due to gravity and then reaches the equilibrium point. Therefore, this limitation may not be a significant restriction when applying the proposed method to prototype systems. Finally, the film flow theory was developed in lab scale experiments and may cause errors when dealing with large scale issues.

### Pressure oscillations with film flow

It has been discussed that the water column length decreases over time during the process and the pressure oscillation varies. When the water column length is less than *kH*_{atm}, the period decreases when the water column losses mass. However, the analytical solution cannot provide an explicit relationship between the changing water column length and the oscillation period. Therefore, numerical method was used to assess the relationship.

Examples of numerically solved pressure variations for Run B2 are shown in Figure 5. Due to a variety of boundary conditions, the pressure variation appears as different patterns. Figure 5(a) shows a typical pressure oscillation pattern where the pressure oscillates with a decreasing period over time due to the mass loss via the film flow. The period *T*_{0} can be clearly found as the second point where the derivative of pressure over time equals to the initial one with the first one as *T*_{1}. The time duration for the first full pressure cycle (i.e. *T*_{0}) is defined as the period of pressure oscillation for further analysis. For Figure 5(b), the pressure keeps decreasing with oscillations. In this case, the first pressure plateau at time *T*_{0} is treated as the oscillation period. Figure 5(c) shows the pressure pattern where the pressure oscillates but with limited cycles. In this case, the water column loses mass due to film flow and the pressure oscillation period decreases so fast that the pressure cannot oscillate more than one cycle before the length of the water column becomes zero. For the pressure oscillation pattern shown in Figure 5(d), the pressure varies less than one full cycle and, therefore, there is no period for these cases. For Figure 5(e), the pressure decreases due to the high-pressure head pushing the water column moving up until the length of the water column reaches zero. In these cases, the pressure does not oscillate and there is no oscillation period.

For the pressure cycle pattern in Figure 5(c), it is defined that if the difference between 0.5 *T*_{0} and *T*_{1} is larger than 0.5 (i.e. ), the period calculated is skewed by the mass loss of the water column due to the film flow and the oscillation period is ignored for further analysis. Out of the 2,401 runs in Run B2, there were 1,687 runs that resulted oscillation periods. The other 714 runs resulted in ‘no oscillation period’ due to either the fact that the length of the water column reached zero before the pressure completed a full oscillation cycle or the oscillation was skewed by the loss of water due to the film flow.

*R*

^{2}= 0.99. Combining Equation (10) and the film flow rate equation, letting the time for one oscillation cycle equals to the time for the water column to lose its mass, the critical water column length (

*h*) for the no oscillation period cases can be written as:where

_{cr}*T*

_{a}is the oscillation period calculated using Equation (10).

*h*

_{0}<

*h*

_{cr}means that the length of the water column decreases to zero before the pressure oscillates for one full period (Figure 5(d) and 5(e)). Out of the 714 ‘no oscillation period’ cases, there were 556 cases fell in this category. The rest 158 cases resulted in ‘no oscillation period’ due to the skewed oscillation period (Figure 6(c)). For these cases, it was noticed that the relative change of the water column after first half oscillation period was higher than 0.3. Therefore, for this scenario, the criteria can be written as:

The critical water column length for generating pressure oscillation is the higher *h*_{cr} value calculated by Equations (16) and (17). When applying the proposed method to engineering applications for the pressure oscillation period, one may firstly use Equations (16) and (17) to determine if the oscillation occurs and then use Equation (10) to calculate the analytical oscillation period. After, the oscillation period is to be multiplied by 0.96 to consider the mass loss due to the film flow.

The analysis above shows a promising approach on predicting the pressure oscillation period in risers when a water column was pushed up by an air pocket. In a prototype system, with the existence of the upstream and downstream pipes and the net flow rate into the system which can be either positive or negative, the dynamics of the process may change and the result from Equation (10) may need to be adjusted accordingly. Figure 6(b) shows the statistic of the calculated 22 periods for Run C2 where *Q*_{b} ranged from −1,000 L/s to 1,000 L/s. It was found that the changing inflow rate did not significantly contribute on the pressure oscillation period. The simulated pressure oscillation for the given condition ranged from 1.15 s to 1.25 s averaging at 1.19 s with a standard deviation of 0.026 s. Comparing with the oscillation period calculated using Equation (10) at 1.29 s, the difference was 8%. The correction factor from the analytical solution to the numerical solution was 0.92 which was close to 0.96 obtained above in earlier analysis. Therefore, it can be concluded that the changing flow rate did not significantly affect the pressure oscillation. The oscillation period was dominated by the initial air pocket size, initial driving pressure, diameter of the riser, and the water column length in the riser.

## CONCLUSIONS

In the present study, a theoretical model was proposed for predicting the pressure oscillation period. The model was solved analytically without considering the film flow and numerically considering the film flow and inflow rate change. After analyzing the results in detail, the following conclusion can be composed.

Five parameters have been identified to dominate the oscillation period (*T*). These are the initial air pocket volume (), initial air pocket pressure head (*H*_{0}), water column length (*h*), the diameter of the riser (*D _{r}*), and the polytropic coefficient (

*k*). The oscillation period explicitly changes with in the power of 0.5,

*H*

_{0}

^{*}in the power of , and

*A*

_{r}in the power of −0.5. The oscillation period reaches the maximum when the water column length (

*h*) equals to

*kH*

_{atm}. For an increased

*k*, the maximum

*T*at

*h*=

*kH*

_{atm}decreases. The relationship presented above was from the analytical solution which was only evaluated with limited cases in lab scale and the overall parameters tested ranged wider than the evaluated physical conditions, which may bring some uncertainties when dealing with prototype scale issues.

The numerical model considering the mass loss of the water column due to the film flow along the riser was further developed. It was found that the pressure oscillation period decreased slightly to 0.96 times of the analytical solution. Cases not completing an oscillation cycle were due to the water column length was less than the critical value, and the criteria were proposed. With respect to the inflow rate changing cases, it was noticed that it did not affect significantly on the pressure oscillation period. The obtained oscillation period mainly depended on the factors that concluded above. The film flow theory was developed in lab scale experiments and may cause errors when dealing with large scale issues. Nevertheless, this study provides a general method on predicting the oscillation period induced by the expansion/compression of an air pocket below a water column in a vertical riser.

## ACKNOWLEDGEMENT

The writers gratefully acknowledge the financial support from the Natural Sciences and Engineering Research Council (NSERC) of Canada.

## DATA AVAILABILITY STATEMENT

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