Multiphase Rapid Filling Conditions of Tunnel System in Columbus, Ohio
Abstract
The City of Columbus, Ohio is implementing a tunnel system to reduce the number of episodes of combined sewer overflows into the Scioto River. The tunnel systems provide relief to the existing Olentangy Scioto Interceptor Sewer. Two new tunnels being implemented are the OSIS Augmentation and Relief Sewer (OARS), in service since July 2017, and the Lower Olentangy Tunnel (LOT) that is planned to be in service in 2025. The performance of these tunnels in respect to high inflow conditions was investigated with the use of the HAST mixed flow model and the OpenFOAM CFD model to determine the magnitude of surges, the possibility of air pocket entrapment, air–water surging, and the consequences of uncontrolled air pocket releases through shafts. Inflows into the systems were obtained from a calibrated collection system SWMM model. Modeling results quantified surging in the tunnel dropshafts and their mitigation from built-in surge control chambers. HAST simulations also pointed to locations where air pockets could form. These results were used in OpenFOAM to determine the effects of uncontrolled air release through the shaft that links the two tunnels. It was shown that proper ventilation at the shaft will mitigate the growth of air phase pressure to damaging levels.
1 Introduction and objectives
Stormwater management is a very complex task in highly urbanized areas and involves a combination of draining excess runoff efficiently while minimizing environmental impacts to the receiving water bodies. Such a task is more complex when storm sewers and wastewater sewers are combined in the same system. Around 860 municipalities in the United States have combined systems (EPA 2020), and during intense rain events episodes of discharge of combined sewage in waterbodies occur that increase the severity of the environmental impacts.
The City of Columbus, Ohio is one municipality that operates combined sewer systems, and it is implementing a long term control plan to reduce the frequency of CSO discharges into the Scioto River. A key component of Columbus infrastructure is the Olentangy sewer interceptor system (OSIS). Given the limited availability of storage at grade level, OSIS storage is augmented using the OSIS augmentation and relief system (OARS) and the Lower Olentangy tunnel (LOT) to reduce the frequency of CSO episodes. With these tunnels systems, the City of Columbus aims to prevent two billion gallons (910 ha m) of combined sewage from reaching the Scioto River without treatment (Pratt 2019).
While deep tunnel systems are a technically sound approach to manage large flows in highly urbanized areas, there are some design concerns over their implementation. Since the implementation of these tunnels in the late 1970s and 1980s, researchers have pointed to potential problems linked to surging in these systems (Song et al. 1983; Guo and Song 1990). More recent research has also recognized the important role of the entrapped air phase within conduits, such as in air–water surging (Zhou et al. 2002; Vasconcelos and Leite 2012).
Much is still unknown with regards to air–water interactions in stormwater systems. There are various uncertainties with regards to the mechanisms for air pocket formation in closed conduits (Vasconcelos and Wright 2006; Schulz et al. 2020). Once entrapped, the fate of air pockets is still largely unknown. One hypothesis has been to assume that entrapped air pockets move in closed conduits as discrete gravity currents (Chosie et al. 2014; Hatcher et al. 2014). Upon reaching vertical shafts and other ventilation points, air will be released in an uncontrolled fashion, potentially leading to structural damage (Zhou et al. 2002), geysering episodes (Wright et al. 2011; Vasconcelos and Wright 2011; Muller et al. 2017) and slab or access cover displacement (Wang and Vasconcelos 2020).
Few numerical models able to simulate rapid filling in stormwater tunnel systems have explicitly incorporated air–water interactions, and this has significantly affected overall model accuracy (Vasconcelos et al. 2015). Considering all these aspects, the simulation of the rapid filling conditions in the OSIS–OARS–LOT system, which is summarized in this paper, used state-of-the-art methodology. First, a 1D model able to track air pocket motion was used to predict the potential regions of air pockets formations within the system. Then, facing the insufficient capability of these models to simulate multiphase air pockets discharge through vertical structures, a 3D model was used for predicting the pressure variations within these structures. The reason for a limited usage of 3D model is that, to date, it is computationally infeasible to use a 3D model in system-wide applications. Therefore, this sequential application of 1D and 3D models combines a reasonable approach to identifying pocket formation and anticipating potential issues associated with the uncontrolled air release in complex structures.
2 Method
2.1 Geometry of OSIS–OARS–LOT tunnel system
The connections between OSIS and OARS and LOT are through several dropshaft structures as the tunnels are deep underground. In 2020, the connection between OARS and OSIS was made through three dropshafts at locations referred to as shafts 4, 5 and 6. A future drop shaft connection (shaft 3) is scheduled to be in service by the year 2025. When LOT is in service, by 2025, OSIS will connect with LOT at three dropshafts known as the Tuttle, Gowdy and Vine shafts. Vine shaft is also connected with OARS at shaft 6 through a 3.66 m diameter connector. The schematics of these tunnels are presented in Figure 1.
Figure 1 Schematic of OARS and LOT tunnel systems and their points of connection with OSIS.
Most of OSIS is of circular cross section that varies from 3.1 m to 3.6 m in diameter, transiting to a 3.2 m × 4.2 m rectangular cross section toward the downstream end of the combined flow portion. The cross section of OARS is also circular, with a much larger diameter of 6.10 m through the entire extension of 7.05 km. The LOT cross section is circular, with a diameter of 3.66 m and length of 5.12 km. Detailed geometry of OARS shaft 6 and the LOT Vine shaft are presented here as it is relevant to the CFD modeling component of this study. Geometry from remaining junctions are not included here for brevity.
Shaft 6, shown in Figure 2, is a large structure with 18.4 m diameter located at the upstream end of OARS. It serves four different purposes: inflow admission through an approaching channel and a 4.9 m diameter vortex structure; air ventilation through a 6.10 m diameter vertical shaft; surge relief through a built-in overflow chamber with area ~270 m2; and a near-future connection with the LOT system.
Figure 2 Schematic of OARS shaft 6, including the dropshaft, surge chamber, ventilation shaft and future connection to LOT.
Vine shaft, presented in Figure 3, has a 33.5 m deep lower shaft with a diameter of 7.6 m, which is divided between a baffled shaft to dissipate inflow energy and a ventilation shaft for air discharge. The upper portion of Vine shaft is 21.3 m high and has a 15.2 m diameter. It is fitted with control gates to regulate inflows from LOT into OARS, and a flap gate to prevent backflow from OARS into LOT.
Figure 3 Schematic of Vine shaft in the downstream end of LOT and its downstream connection to OARS.
2.2 Numerical Simulation
OSIS–OARS–LOT performance was evaluated using 1D system-wide modeling and 3D computational fluid dynamics (CFD) modeling of Vine shaft. The former used the Hydraulic Analysis of Stormwater Tunnels (HAST) model, which is based on the solution of the Saint-Venant equations. The main features of the numerical solution used in HAST are:
- The use of a nonlinear numerical scheme based on the Roe scheme, as presented in Macchione and Morelli (2003), renders a robust and accurate solution even with low local Courant numbers, which is typical for flow regime transition conditions.
- Enabling the Saint-Venant equations to represent both free surface and pressurized regimes, through application of the two-component pressure approach (Vasconcelos et al. 2006).
- Explicit treatment of regions where air pockets appear in closed conduits, tracking compression and expansion processes, spreading, motion and release through junctions.
- Representation of the unique geometric characteristics of junctions in the tunnels through construction of custom boundary conditions in the model as required.
To date, most mixed flow models do not explicitly track the formation and interactions of the air phase in closed conduits. As we show in this work, this feature was very important in determining potential adverse air–water interactions within the tunnels. The mathematical implementation of the model can be found in Vasconcelos et al. (2015). For OSIS–OARS–LOT modeling, 59 junctions were represented in the model, with 63 reaches, totaling 26.5 km of conduits. The initial conditions of the simulation considered the tunnel initially empty, and inflow hydrographs were admitted at selected locations according to the results of surface hydrological modeling. Each junction within the system was represented through specific boundary conditions. These boundary conditions apply equations that include mass and energy balance, characteristic equations, and empirically-derived rating curves. An in-depth discussion of these boundary conditions is outside of the scope of this work. After discretization, 4127 computational cells were used in the model, with time steps in the order of few milliseconds once the pressurization conditions developed.
A CFD tool was used to describe the effects of a potential air pocket release through Vine shaft from OARS into LOT. Since the work of Choi et al. (2014), there have been many subsequent studies using relatively simple geometries, such as Muller et al. (2017) and Qian et al. (2020). However, results in these more fundamental studies are difficult to generalize for physical shaft geometries such as Vine shaft. This motivated the development of a CFD model to assess the effects of air release through that structure.
The CFD model was developed using OpenFOAM, an open source C++ object-oriented library that can perform multiphase flow simulation. In this application, Navier–Stokes equations were solved, tracking the free surface in using the volume of fluid (VOF) method (Hirt and Nichols 1981). There are various multiphase flow solvers within OpenFOAM, and here the continuity, momentum, and energy equations were solved using the compressibleInterFoam solver (Svenungsson 2016). Equations 1, 2 and 3 present the two-phase continuity, momentum, and energy expressions solved in OpenFOAM. Equation 4 is used to track the free surface and is modified to reduce issues associated with the convection of a step function.
(1) | |
(2) | |
(3) | |
(4) |
where:
ρ | = | fluid density, |
t | = | time, |
U | = | 3D velocity vector, |
p | = | pressure, |
µ | = | dynamic viscosity, |
SU | = | momentum source term, |
Cp | = | specific heat, |
T | = | temperature, |
ST | = | energy source term, |
α | = | volume fraction (0≤α≤1), and |
Ur | = | velocity field to compress the interface. |
In the CFD formulation, the values of α represent the liquid fraction within a cell, with unity representing pure water and zero pure air. The velocity field Ur is included to counter a disadvantage of VOF in solving for the free surface that is related to interface smearing near the free surface, as described by Rusche (2003). Turbulence was represented initially in this research using a k–ε model but subsequently a k–Ω SST model was used.
Mesh generation was completed using the snappyHexMesh utility supplied by OpenFOAM. This tool generates 3D meshes that are in the stereolithography (STL) or wavefront object (OBJ) format from triangular surface geometries (Greenshields 2015). Since different geometries were evaluated, mesh details will be further discussed in section 4, OpenFOAM modeling results. CFD modeling was performed in different stages, which varied according to the volume of air released in the tunnel–shaft system, initial water volume, or ventilation. These values were selected with the intention of representing different scenarios ranging from less critical (e.g. 12.5 m of water level and 209 m³ of air volume) to highly critical (e.g. 35 m of water level and 576 m³ of air volume) in regions where the HAST model predicted a tendency for formation of air pockets. A total of 14 simulation conditions were tested, and these are presented in Table 1.
Table 1 CFD study variables and ranges of variation.
Variable | Range of variation considered in the modeling |
Initial water level (m) | 12.5, 25.0 and 35.0 |
Ventilation conditions | Ideal ventilation, fully blocked, 2 ports 2.4 m × 2.4 m |
Initial air volume in system | Various, ranging from 209 m3 to 576 m3 |
Turbulence model | k–ε model but subsequently the k–Ω SST model |
Number of computational cells in CFD model | Up to 5.09 million |
For various system-wide 1D rapid filling simulations performed with HAST, we noticed a tendency for air pocket formation in the tunnel reaches between shaft 5 and shaft 6. Given the slope of the tunnel and that there is no route for air escape in shaft 6, entrapped air pockets in this region are likely to be discharged through the 6.1 m diameter ventilation shaft of shaft 6. Muller et al. (2017) demonstrated using CFD that it is possible that, even with such large ventilation, there was a residual volume of air that could not be ventilated in the shaft and remained in the system. In the present case, this air could be routed through the connection between OARS and LOT and through Vine shaft, as illustrated in Figure 4.
Figure 4 Potential air pocket entrapment in OARS and release of a fraction of the air pocket in Vine shaft.
3 HAST modeling results
Inflows used in the HAST modeling were obtained through a calibrated SWMM model with a time series spanning 20 y. The City of Columbus model was recently updated using the innovative Model at the Source approach (Gheith 2014; Herr and Gheith 2015), which leads to a high quality flow prediction compared to observed data. Inflows from 10 rain events were selected based on the flows that were most likely to generate the strongest surges. According to Vasconcelos and Wright (2017), such conditions are expected when inflows are highest while horizontal storage in tunnels is depleted (i.e. full pressurization of horizontal reaches). Table 2 lists the events in a ranked order based on higher to lower inflows that were modeled in HAST based on the summation of inflows at the time when the tunnel reaches were pressurized. The normalized value for the tunnel system diameter of 3.85 m was computed using the diameters DR and lengths LR of each tunnel reach through Equation 5:
(5) |
Table 2 Rank of inflow events from OSIS to OARS and LOT when pressurization of the tunnel reaches occurred.
Rank | Date | Q inflow (Qi) at pressurization (m3/s) | Normalized Qi = Qi/(gD5)0.5 |
1 | 2003-08-30 | 59.8 | 0.66 |
2 | 2009-08-28 | 56.1 | 0.62 |
3 | 2004-06-11 | 37.3 | 0.41 |
4 | 2005-01-11 | 30.1 | 0.33 |
5 | 1997-05-13 | 25.7 | 0.28 |
6 | 2000-01-03 | 24.5 | 0.27 |
7 | 1998-12-21 | 22.6 | 0.25 |
8 | 2004-01-03 | 18.5 | 0.20 |
9 | 2011-12-05 | 14.6 | 0.16 |
10 | 1995-01-15 | 13.6 | 0.15 |
In the early stages of a storm event, OSIS begins to fill before flows start to divert into OARS and LOT over relief weirs located upstream of the dropshafts. In normal conditions, as OARS is filling, the reaches between the downstream end and shaft 6 become pressurized. Typically, the upstream end of OARS at shaft 6 is the last portion to undergo pressurization, as illustrated in Figure 5.
Figure 5 HGL in OARS, LOT and a portion of OSIS immediately prior to the full pressurization of OARS (2005-01-11).
Following the rapid filling of the horizontal reaches, deep tunnels begin to experience a rapid rise of water in the vertical structures. If this filling is too fast, surging processes can develop with a risk of water reaching grade if adequate design does not control the formation of surge conditions. The dropshaft structures in OARS shafts 4, 5 and 6 were designed with surge chambers to control the surge propagation. The performance of the surge chamber at shaft 6 during a large storm event is illustrated in Figure 6. HAST computed that up to 6.1 m3/s is diverted from the dropshaft into the surge chamber, causing a significant attenuation in the surge as it propagated up to the surge chamber weir crest.
Figure 6 Water level at shaft 6 and surge chamber together with flow rate diverted to the surge chamber (2003-08-30 event).
One aspect of LOT design that was evaluated using HAST was the benefits of including a flap gate at Vine shaft to prevent OARS from backing up into LOT. Simulation results from HAST indicated that without flap gates the filling of LOT would occur at least 1.5 h sooner, with potential propagation of surge in OARS into LOT. Alternatively, with flap gates, the two systems are disconnected and filling in LOT would occur due to inflows diverted directly into it. The differences between these two outcomes, computed with HAST, are presented in Figure 7.
As pointed out earlier, in different simulation conditions we noticed a tendency of air pockets to form at particular locations, including near shaft 6. Despite the ventilation in shaft 6, we hypothesized that a large pocket could form in this location and not be fully ventilated. In such an event, the connecting tunnel between shaft 6 and Vine shaft in LOT would be a possible location for air to escape to Vine shaft. The CFD study was performed to assess the effects of an uncontrolled release of a large air pocket through Vine shaft.
Figure 7 Water levels upstream (US) and downstream (DS) of Vine shaft for scenarios with or without flap gates (2005-01-11 event).
4 OpenFOAM modeling results
The first step in the OpenFOAM CFD modeling of the uncontrolled release of air pockets through Vine shaft involved the development of a mesh using the tool SnappyHexMesh. Mesh development was incremental with respect to the ventilation conditions in Vine shaft. At this point it is important to highlight that a systematic mesh convergence analysis was not performed in this work due to time constraints. However, several meshes were tested in order to guarantee the grid independence. The initial mesh development iterations included perfect ventilation, in which the top slab did not exist. Such simulations were intended to assess whether the release of air could potentially raise the water level within the shaft above grade. As it is shown in Figures 8 and 9, the release of large air pockets through a water filled Vine shaft did raise the free surface within Vine shaft. However, through all the simulated conditions, the change in the water level was not large enough to reach grade elevation and trigger geysering. This is attributed to the much larger plan area of Vine shaft at elevations that are close to grade.
The second group of CFD simulations assumed zero ventilation conditions in Vine shaft, so vertical water motion created by the air pocket release caused compression of the air in the headspace under the slab. Figure 10 shows air pressure at the headspace for different air pocket sizes. Large air pockets were released in separate segments, creating separate air phase growth spurts. Of most importance is that the magnitude of the air phase pressure head, depending on the scenario, can increase to up to 6 m of water in the absence of ventilation, which in turn could create potentially significant forces in the slab of Vine shaft.
Figure 8 Air pocket release impact on water level in Vine shaft for a 25 m initial water level assuming no top slab.
Figure 9 Impact of air pocket release on water level in Vine shaft for a 35 m initial water level assuming no top slab.
Figure 10 Air pressure head under Vine shaft slab in different scenarios of uncontrolled air pocket release if no ventilation is available.
The last CFD simulation considered the most adverse conditions for headspace air pressurization when a proposed ventilation in Vine shaft was added. As shown in Figure 11, two ports with dimensions 2.4 m × 2.4 m near the top of Vine shaft were proposed in the LOT design. This modeling is the most representative condition for the uncontrolled release of air in Vine shaft.
Figure 11 Geometric details of the ventilation underneath the top slab at Vine shaft.
Figure 12 shows the CFD simulations of the air phase pressure head, expressed in meters of water, underneath the slab and at the two ventilation ports. As can be seen in the figure, the release of the air phase at the bottom of Vine shaft was initiated at time T = 15 s. Within 3 s, the free surface at the top of Vine shaft began to experience additional pressure because of the water displacement created by the air pocket. Pressures increased very rapidly, achieving a peak at T = 18.4 s, after which a sudden drop occurred as the water interface collapsed after the pocket opening. The values of the peak pressures at the three locations varied slightly, with the maximum of 2.3 m under the Vine shaft slab.
Figure 12 Air phase pressure results in Vine shaft with proposed ventilation.
Simulation results of the pressures on the baffles indicated in general an immediate drop as the air pocket started rising in Vine shaft, as shown in Figure 13, which is consistent with previous studies (Muller et al. 2017). As the air pocket reached the surface, pressures were minimum, and a series of pressure oscillations caused by mass oscillations in Vine shaft are seen.
Figure 13 Water pressure on different baffle structures in the Vine shaft during air release.
5 Final remarks and conclusion
From the development of these studies that involved the numerical simulation of different rapid filling scenarios of the OSIS–OARS–LOT systems, the following conclusion are drawn:
- Rapid filling conditions underlie the complex interactions between junctions and reaches in these tunnel systems. Adequate modeling using tools such as HAST enable the accurate quantification of surges in the tunnel junctions and dropshafts. The new OARS and LOT systems did not show excessive surging as the structures were adequately designed to mitigate such events.
- Despite the adequate sizing of junctions and tunnels, numerical simulation of the filling processes indicated possible air pocket formation at specific locations in OARS, LOT and OSIS. In particular, the formation of air pockets near shaft 6 can potentially lead to air release through shaft 6 in OARS and Vine shaft in LOT.
- The behavior of Vine shaft in scenarios of uncontrolled air pocket release was evaluated using a state-of-the-art OpenFOAM CFD modeling approach. Our work indicated that the occurrence of geysering is unlikely at that location. However, this study also demonstrated the need for an adequate ventilation structure at that point to avoid potential damage created by air compression underneath the Vine shaft top slab.
We hypothesize that the types of flow conditions and air–water interactions during intense rain events are more widespread than is realized by designers and city officials. However, as monitoring tools and design strategies improve, such adverse air–water interactions will become better recognized. The development of modeling studies such as this offers one possible path for designing future systems or evaluating existing systems. It is expected that with the improvement of this method, stormwater systems will achieve better operational conditions and greater resiliency.
Acknowledgements
This work was supported in part by the U.S, National Science Foundation under Grant No. 2048607.
References
- Choi, Y.J., A.S. Leon, and S.V. Apte. 2014. "Three-dimensional numerical modeling of air–water geyser flows". In World Environmental and Water Resources Congress 2014, 1535–48.
- Chosie, C.D., T.M. Hatcher, and J.G. Vasconcelos. 2014. “Experimental and numerical investigation on the motion of discrete air pockets in pressurized water flows.” Journal of Hydraulic Engineering 140 (8): 04014038
- EPA 2020. Combined Sewer Overflows (CSOs). Washington, DC: U.S. Environmental Protection Agency. https://www.epa.gov/npdes/combined-sewer-overflows-csos
- Gheith, H.M. 2014. “An Innovative Approach to Quantify RDII Sources Using SWMM Groundwater Module", Presentation at the International Conference of Water Management Modeling, Toronto, ON, February 26-27, 2014 (abstract only). https://www.icwmm.org/Archive/2014-C023-33/innovative-approach-to-quantify-rdii-sources-using-swmm-groundwater-module
- Greenshields, C.J. 2015. OpenFOAM User Guide. OpenFOAM Foundation Ltd, 230.
- Guo, Q., and C.C. Song. 1990. "Surging in urban storm drainage systems". Journal of Hydraulic Engineering 116 (12): 1523–37.
- Hatcher, T.M., and J.G. Vasconcelos. 2014. "Alternatives for flow solution at the leading edge of gravity currents using the shallow water equations". Journal of Hydraulic Research 52 (2): 228–40.
- Herr, R., and H. Gheith. 2015. “Blueprint Columbus collection system model.” OWEA Buckeye Bulletin 2: 60–2.
- Hirt, C.W. and B.D. Nichols. 1981. “Volume of fluid (VOF) method for the dynamics of free boundaries.” Journal of Computational Physics 39 (1): 201–25.
- Macchione, F., and M.A. Morelli. 2003. “Practical aspects in comparing the shock-capturing schemes for dam break problems.” Journal of Hydraulic Engineering 129 (3): 187–95.
- Muller, K.Z., J. Wang, and J.G. Vasconcelos. 2017. “Water displacement in shafts and geysering created by uncontrolled air pocket releases.” Journal of Hydraulic Engineering 143 (10): 1–13.
- Pratt, N. 2019. “OARS–OSIS augmentation & relief sewer tabbed as OCEA finalist.” ASCE News. https://news.asce.org/oars-osis-augmentation-relief-sewer-tabbed-as-ocea-finalist/
- Qian, Y., D. Zhu, L. Liu, W. Shao, S. Edwini-Bonsu, and F. Zhou. 2020. “Numerical and experimental study on mitigation of storm geysers in Edmonton, Alberta, Canada.” Journal of Hydraulic Engineering 146 (3): 04019069.
- Rusche, H. 2003. Computational fluid dynamics of dispersed two-phase flows at high phase fractions. London: Imperial College, University of London. PhD Dissertation.
- Schulz, H.E., J.G. Vasconcelos, and A.C. Patrick. 2020. “Air entrainment in pipe-filling bores and pressurization interfaces.” Journal of Hydraulic Engineering 146 (2): 04019053.
- Song, C.C., J.A. Cardie, and K.D. Leung. 1983. “Transient mixed-flow models for storm sewers.” Journal of Hydraulic Engineering 109 (11): 1487–504.
- Svenungsson, J. 2016. “Solving electric field using Maxwell’s equations and CompressibleInterFoam Solver.” CFD with OpenSource Software.
- Vasconcelos, J.G., and S.J. Wright. 2006. “Mechanisms for air pocket entrapment in stormwater storage tunnels.” In World Environmental and Water Resources Congress 2006, 40856 (1999), 9.
- Vasconcelos, J.G., and S.J. Wright. 2011. “Geysering generated by large air pockets released through water-filled ventilation shafts.” Journal of Hydraulic Engineering 137 (5): 543–55.
- Vasconcelos, J.G., and G.M. Leite. 2012. “Pressure surges following sudden air pocket entrapment in stormwater tunnels.” Journal of Hydraulic Engineering 138 (12): 1081–9.
- Vasconcelos, J.G., and S.J. Wright. 2017. “Anticipating transient problems during the rapid filling of deep stormwater storage tunnel systems.” Journal of Hydraulic Engineering 143 (3): 06016025
- Vasconcelos, J.G., P.K. Klaver, and D.J. Lautenbach. 2015. “Flow regime transition simulation incorporating entrapped air pocket effects.” Urban Water Journal.
- Wang, J., and J.G. Vasconcelos. 2020. “Investigation of manhole cover displacement during rapid filling of stormwater systems.” Journal of Hydraulic Engineering 146 (4): 04020022.
- Wright, S.J., J.W. Lewis, and J.G. Vasconcelos. 2011. "Geysering in rapidly filling stormwater tunnels". Journal of Hydraulic Engineering 137 (1): 112–5.
- Zhou, F., F.E. Hicks, and P.M. Steffler. 2002. “Transient flow in a rapidly filling horizontal pipe containing trapped air.” Journal of Hydraulic Engineering 128 (6): 625–34.