Air entrapment modelling during pipe filling based on SWMM

Abstract

The paper proposes a novel methodology to locate and quantify entrapped air pockets created during pipe-filling events often found in intermittent water supply systems. Different filling conditions were tested in an experimental pipe with a high point. Measurements were taken and video recordings were carried out to assess air pocket volumes for different air release conditions at the downstream end of the pipe. The stochastic nature of air pocket creation resulted in varying air volumes. A new numerical model capable of simulating the air pocket creation, dragging and entrainment has been proposed. The new model, AirSWMM, was implemented as an extension of the Stormwater Management Model (SWMM) with stochasticity of air pocket formation reproduced by simulations with different air entrainment rates. The obtained numerical results show that the proposed model, even though based on a single-phase one-dimensional flow, can accurately locate and approximately quantify the entrapped air pocket volumes.

KEYWORDS:

• Air–water interface interaction
• air pocket entrapment
• flow visualization and imaging
• intermittent water supply
• SWMM
• water pipelines

1 Introduction

The presence of air severely affects water supply, urban drainage and stormwater systems (Fuertes-Miquel et al., 2019). Intermittent water supply (IWS) systems are particularly affected by air since they are not always pressurized, i.e. they continuously go through a cycle of filling, supplying and emptying stages. The filling stage is characterized by water abruptly entering into the system, forcing the air inside to be released. Conversely, the air release rate influences the flow rate going into the system. The air release delays the pipe filling and affects the water supply time to end consumers. During the filling stage (and also due to possible inadequate air valve design and maintenance), air pockets can get entrapped in the pipes. These pockets can cause pipe malfunctioning and local head losses (Ramezani et al., 2016) and can induce high pressure variations during transient events (Erickson et al., 2022; Ferreira et al., 2021; Martins et al., 2015; Martins et al., 2017). This behaviour has also been observed in water supply systems after the implementation of IWS: increased pipe failures and higher leakage levels are observed after a water supply system starts operating intermittently (Charalambous & Laspidou, 2017; Christodoulou & Agathokleous, 2012).

Despite all the associated problems with air pockets, the present state of engineering practice does not offer a numerical methodology to detect, locate and quantify air pockets using more traditional one-dimensional models. Work has been developed to locate and quantify entrapped air pocket volumes in a bridging pipe using three-dimensional (3D) computational fluid dynamics (CFD) model (Kaur et al., 2023). However, using such models for network assessment and engineering practice is not a viable solution due to their computational effort and specialist knowledge required to set up and use these models.

Stormwater Management Model (SWMM) is a one-dimensional (1D) model developed for urban drainage and stormwater design and management that has been used as a computationally inexpensive tool to simulate IWS systems (Cabrera-Bejar & Tzatchkov, 2009). Campisano et al. (2019) investigated pipe filling using SWMM and obtained a satisfactory agreement between field data and numerical results. However, the model does not incorporate the air phase during the filling stage. Gullotta et al. (2021) expanded the use of SWMM to determine the optimal location of pressure-reducing valves to equitably distribute water under IWS conditions. Still, this work did not consider the air phase. A method has been developed to identify possible entrapped air pockets locations using SWMM, but overestimated the entrapped air pocket volumes (Ferreira et al., 2022). An air accumulator model has been incorporated in SWMM to simulate the air pressure and density with a piston equation to track the waterfront’s position (Ferreira et al., 2023); this model allows simulation of the air pressurization in pipes but is not able to detect or quantify air pockets due to the assumption of a perpendicular waterfront to the pipe axis.

This paper presents and proves a novel methodology for improved one-dimensional modelling of air pockets during pipe filling events in IWS systems. Unlike in Ferreira et al. (2023), this methodology incorporates a more robust waterfront tracking method and two-phase flow dynamics mechanisms, i.e. pocket creation, dragging and entrainment. These, in turn, enable the new methodology presented here to locate air pockets and quantify their respective volumes, which was not possible in the previously published method. Novel experimental tests are conducted in an existing pipe rig to understand the process of air pocket creation and volume variation at a high point in the system. Pressure and flow rate measurements as well as video recordings are carried out. Image processing is used to estimate the actual entrapped air pocket volumes and air travelling in the pipe. Collected data are used to calibrate and validate the new air entrapment model.

The outline of the paper is as follows. Section 2 describes the experimental data collection and corresponding analysis. Section 3 gives a summary of the original SWMM model and proposes the methodology to locate and quantify entrapped air pockets. Section 4 presents the model assumptions, describes the calibration process and shows the validation results by using the proposed methodology. A summary of the results and conclusions are presented in Section 5.

2 Experimental data collection and analysis

2.1 Experimental set-up and instrumentation

The experimental rig is composed of a straight acrylic pipe, an elevated water tank and a pneumatically actuated valve (Figure 1). The pipe has an inner diameter of $D_{\text{p}}$ = 21 mm with wall thickness of 2 mm, is 12.5 m long and it features a high point with a maximum elevation of 0.1 m, with the rising pipe 4.5 m distant from the actuated valve. The pipe segments leading to the high point are sloped with a 45° angle. The pipe is longitudinally supported and fixed to minimize its movement. The system is gravity driven, supplied by the 50 l tank located at the upstream end. A pneumatically actuated quarter-turn ball valve, with an inner diameter of 20 mm, is located at 0.2 m downstream of the water tank. Acrylic plates with centrally drilled orifices with diameters of d = 2.2, 3.0 and 4.5 mm are installed at the downstream end of the pipe. They are used to simulate a contraction in the flow cross-section just before the water is discharged into the atmosphere. Additionally, the system is tested without the presence of an orifice at the downstream end of the pipe, simulating a fully free discharge into the atmosphere. Orifice sizes larger than d = 4.5 mm are not tested as they behave as if the pipe was almost fully open, not showing relevant air pressure variations (Ferreira et al., 2023).

Figure 1. (a) Experimental rig and (b) downstream section schematic. Photographs of (c) the downstream boundary condition and (d) high point.

Three types of observations are carried out: pressure and flow rate measurements and video recordings. Pressure measurements are carried out using Siemens (Nürnberg, Germany) SITRANS P series Z pressure transducers with a range of 0–2.5 m, with a 0.5% full-scale error and a response time lower than 0.1 s. Four of these transducers are installed in the pipe (Figure 1): PT1 is located immediately upstream the actuated valve to monitor the initial tank head; PT2 and PT3 are located before and after the high point of the system, at 4.15 and 6.2 m from the upstream valve, respectively; and PT4 is located at 10.35 m from the upstream valve. Flow rate measurements are carried out using a Dynasonic (Neuffen, Germany) ultrasonic flowmeter that has a full-scale error of 1% and a response frequency of 50 Hz, installed at 5 m from the actuated valve. All pressure and flow rate measurements are acquired at a 1 kHz frequency. Videos are recorded using a GoPro 7 Black (San Francisco, [California], United States of America) with a resolution of 2074 × 1520 pixels and a frame rate of 24 frames per second.

2.2 Experimental tests

Twenty pipe filling tests are carried out for each downstream orifice size to ensure replicability and observe the stochastic nature of the pipe filling and the respective air pocket sizes. The pipe is empty at the beginning of the experiments and all tests start with a constant water head $H_{\text{0}}$ = 0.5 m in the upstream tank. The valve is opened to start the test with an effective opening time of 0.0021 s (Ferreira et al., 2018), allowing water to flow into the pipe. After the pipe filling is completed and a steady state flow is reached, the test is considered to have finished. After each test the pipe is emptied with compressed air and the upstream tank level is reset to 0.5 m, so that the system is ready for the next test. Water temperature varied less than 0.5°C during all conducted tests.

2.3 Experimental data analyses

2.3.1 Air pocket formation and development

Three types of air pockets are observed during the pipe filling process: travelling air pockets, dynamic entrapped air pockets and steady-state entrapped air pockets. Travelling air pockets are formed during the filling process and are carried with the flow due to the velocity profile’s shape, pipe layout and diameter. These often collide and merge with other entrapped air pockets of the other two categories. Dynamic entrapped air pockets are created or dragged into high points during the pipe filling process but are unable to be dragged by the flow. These pockets can have their air–water interface disrupted from the collision with other pockets or bubbles. Air entrainment is also observed in these pockets’ tail during the pipe filling. Steady state entrapped air pockets are formed by the coalescence of all dynamic and travelling air pockets converging to higher elevation points due to the lack of flow momentum to drag them or to compensate for the air buoyancy after reaching steady-state flow. Steady state entrapped air pockets show clearer and more static boundaries than dynamic entrapped air pockets.

An example of collected pressure-head signal and air pocket images are shown in Figure 2 for each air pocket type for a test with a downstream orifice diameter d = 3.0 mm. Figure 2a shows the pressure-head time series for the sampled pipe filling test at each transducer and Figure 2b shows the corresponding images of the high point.

Figure 2. Test sample d = 3.0 mm and one high point pipe layout: (a) images at different pipe filling moments and (b) pressure-head signal at each pressure transducer.

The images shown in Figure 2a come from the video recordings for the same test. The timestamps of each image are shown in Figure 2b. The high point of the rising pipe is empty until t = 10 s, i.e. the moment when a waterfront reaches this location. Air pockets are observed on the waterfront due to the turbulence of the filling behaviour. At t = 13 s, an air pocket is created after the high point and does not get dragged due to the high pipe’s downward slope, the low flow velocity and the large air pocket size. The filling continues and an unexpected behaviour is observed at t = 17.83 s. The highlighted air bubble coming from upstream collides with the entrapped air pocket and splits into two smaller-sized ones as observed at t = 18.16 s. The air entrainment process keeps occurring at the downstream bottom side of the air pocket in the form of air bubbles until the waterfront reaches the downstream end at t = 19.25 s. The air pockets rise to higher points due to the lack of drag force from the steady state flow and, consequently, take a more stable shape. Air pockets close to the slope change (from downstream sloped pipe to the horizontal pipe) rise to the high point. The air pocket that stays at the high point (dynamic entrapped air pocket) and the one rising from the downwards sloped pipe coalesce and reach an equilibrium when a steady state flow is achieved, t = 25 s. A small local head loss is observed in the pressure-head signal from the difference between PT2 and PT3 due to the air pocket in the high point.

Figure 2b shows pressure-head time variations in four transducers. A pressure-head drop from 0.50 to 0.30 m is observed at t = 0 s in PT1, corresponding to the valve opening. As the water enters the pipe, air becomes pressurized, as observed in the pressure head rise at PT2 – PT4. The latter transducers measure the same pressure-head until t = 7.5 s since these values correspond to the air pressure ahead of the waterfront. The waterfront arrival to each transducer is identified by the pressure-head stabilization or increase. Pressure-head signals diverge when the waterfront reaches each transducer: PT2 is reached around t = 7.5 s, PT3 at t = 13 s and PT4 at t = 16 s, since air pressure decreases as is released and water pressure depends on the upstream water tank head and head losses. Two more pressure-head features are identified. The pressure-head rise at PT2 at t = 10 s corresponds to the waterfront rise towards the high point between PT2 and PT3 by building up pressure to the upstream section of the pipe. The sharp pressure variations following t = 19.5 s correspond to the moment the waterfront hits the downstream acrylic plate and generates a pressure wave (water hammer), thus concluding the filling process. The pressure wave goes back and forth along the pipe until the flow energy is dissipated by friction, heat transfer by air compression and expansion and pipe wall viscoelastic behaviour, reaching a new steady state. From this moment onwards, the flow becomes steady and all entrapped air pockets are formed and attain their final volume.

In this paper, all the comparisons between experimental and numerical entrapped air pocket volumes will refer to steady state air pockets. Initial and evolving air volumes were not quantified because of the lack of a second camera while running the tests. The air pocket splitting and eventual further coalescent phenomena (seen at t = 18.16 and t = 22.17 s in Figure 2b example) have been disregarded for the following reasons. Firstly, air pocket splitting does not occur in all tests, as it only occurs in ca. 35% of the tests, very much depending on the tested orifice size. When this behaviour occurs, most air pockets with a disrupted interface are dragged downstream. Thus, no substantial difference is expected between split and non-split air volumes at the final steady state air pocket volume. Secondly, dynamic air pockets seem to depend on the air bubble size that causes the split, the split air pocket size, its shape and the angle of incidence when both collide. The numerical modelling of this behaviour is not possible using a one-dimensional (1D) model and, therefore, it is out of the proposed work scope. More complex modelling, such as 3D computational fluid dynamics models, is needed to simulate these phenomena.

2.3.2 Flow rate variation

Maximum and steady-state flow rates are higher for larger downstream orifice sizes (Figure 3a and b). Large orifice cross sections allow higher initial air flows and, consequently, higher initial (maximum) water flows. Also, larger orifice sizes lead to lower the local head losses are and, thus, higher (final) steady-state flow rates.

Figure 3. Distribution of the (a) maximum and (b) final steady-state flow rates during the pipe filling process for each downstream orifice size considering all tests; and (c) examples of flow rate time series and (d) corresponding head signal for each downstream orifice size.

Figure 3c shows the flow rate time series for each of the four orifice size tests. As the valve is opened (at t = 0 s), the waterfront takes some time to reach the ultrasonic flowmeter which is located 5 m from the upstream valve, which is why flow rate rise occurs at different times. The flow rate is high during the pipe filling, and sharply decreases to the smallest three orifice sizes (2.2, 3.0 and 4.5 mm). There is a short period when the flow rate fluctuates at the end of the filling due to the pressure wave going back and forth until steady state is reached. No sharp flow rate variation is observed for d = 21 mm because no orifice is at the downstream end (i.e. the downstream end is fully open into the atmosphere). Figure 3d shows the corresponding pressure signal for each orifice size presented in Figure 3c. The pressure in the air increases for smaller orifice sizes. A higher pressure at the initial stage of opening also indicates a slower filling, hence a delayed arrival to the downstream end of the system where the waterfront hits the orifice and creates a pressure peak.

2.3.3 Entrapped air pocket volume

The air pocket volumes at the high point are determined based on the processing of the video recorded images. Air pockets are assumed to be axisymmetric (i.e. of cylindrical shape) around the air pocket axis. The air pocket height, $l_{\mathrm{AP}}$, and diameter, $D_{\mathrm{AP}}$, are obtained once the pixels corresponding to the air pockets are identified by cropping (done using Gaussian filters and binarizing the images). Each air pocket pixel dimension is obtained and the air pocket pixels are converted into lateral air pocket area. This area allows the final entrapped air pocket volume, $V_{\mathrm{AP}}$, to be obtained when the air pocket resembles a cylinder by: (1) $V_{\mathrm{AP}}=\frac{\pi {D_{\mathrm{AP}}}^2}{4}{l}_{\mathrm{AP}}=\frac{\pi D_{\mathrm{AP}}}{2}{S}_{\mathrm{AP}}$(1) where $S_{\mathrm{AP}}$ is the air pocket area in the image given by $S_{\mathrm{AP}}=D_{\mathrm{AP}}{l}_{\mathrm{AP}}$. Air pockets that more resemble a cone (Figure 4a) have their volume calculated by: (2) $V_{\mathrm{AP}}=\left(\pi r_{\mathrm{AP}}{S}_{\mathrm{AP}}\right)/3$(2) where $r_{\mathrm{AP}}$ becomes the radius of the cone base.

Figure 4. Example of image treatment: (a) original image with ROI and (b) detail of processed and binarized image.

Figure 4a presents an example of the original image with the analysed region of interest (ROI) and Figure 4b the final cropped, filtered, binarized and filled air pocket area. The area of each air pocket is calculated by counting the number of white pixels in Figure 4b and converting such value into the air pocket cross section area and volume according to Equations (1) and (2).

Steady state entrapped air pocket volumes show a considerable variation between tests, as depicted in Figure 5a. Maximum and median volumes of final entrapped air pocket decrease with the increase of downstream orifice size (Figure 5a). As the orifice size increases, the less steep the waterfront becomes, originating lower air pocket volumes (Guizani et al., 2006). This is reinforced by the progressively higher maximum filling flow rate, $Q_{\text{max}}$, in the flow rate time series for a sampled tested for each downstream orifice diameter. The maximum filling flow rate decreases as the downstream orifice size decreases, meaning that it takes longer for the waterfront to reach the downstream orifice (note the sharp maximum flow rate drop from d = 2.2 to 4.5 mm in Figure 3a). For the case of the fully open downstream end (d = 21 mm), the average is 280 mm³ and median air pocket volumes is 0. The maximum and minimum entrapped air pockets for d = 2.2 mm are shown in Figure 5b and c.

Figure 5. (a) Air pocket volume boxplots for each orifice size and (b) maximum and (c) minimum air pockets for d = 2.2 mm, across the repeated experiments.

2.3.4 Air pocket dragging

An additional category of air pockets is identified during the experimental tests. Small air pockets are created in the horizontal pipe section during the pipe filling. Their formation is not directly observed in the video recordings but their dragging in the flow over time allows for their identification. Image processing like that used for steady state air pockets is carried out to quantify their volume. Each test’s recordings are analysed to understand the influence of the air release conditions on the small air pocket creation in the horizontal pipe section.

The distribution of the volumes of travelling air pockets for each downstream orifice is presented in Figure 6 in non-dimensional terms as a ratio of the dragged air volume and the upstream pipe volume. Maximum and average dragged air volumes decrease with the downstream orifice diameter increase, whilst the minimum is within the same order of magnitude. This decrease of the volumes can be explained by the waterfront wave becoming gradually steeper and, consequently, less air pockets being created on the pipe crown to be later dragged. Maximum and average values of dragged air volumes increase after d = 4.5 mm, caused by the considerably higher flow rate during the filling process and when the steady state is reached. This means all the air pockets created in the upstream pipe section are effectively dragged, pass by the ROI and are released at the downstream end, whilst that is not possible for the d = 4.5 mm and below. Staggered air pockets are visually observed at the upstream section of the pipe at the end of each test for d = 2.2–4.5 mm but not for d = 21 mm. No videos were recorded other than from the high point.

Figure 6. Normalized travelling air pocket volume for each downstream orifice.

Travelling air pockets can, partially or fully, coalesce with air pockets created at the high points or further down. Ultimately, the travelling air volume is dragged until the high point where it can coalesce with the air pocket there. Image processing was carried out to quantify the coalescence ratio between the dragged volume and the volume that overcomes the high point air volume and goes downstream. However, no conclusive results could be obtained. The air–water mixing behaviour at the air pocket tail and the relatively low contrast image did not allow a clear delimitation to be established between the air and the water in that region. The same reasoning applies to the air entrainment observed at the entrapped air pocket tail, where small air bubbles are observed to detach from the original air pocket and are dragged downstream (Figure 2b, t = 19.25 s). Other researchers managed to quantify the air entrainment in pressurized flows, but for steady state flows; and the air volume was artificially injected in the pipe rather than as a consequence of a pipe filling event (Escarameia, 2007; Kalinske & Bliss, 1943; Wisner et al., 1975), leading to lower air–water mixtures and more observable pockets.

3 Numerical model

3.1 Original SWMM

In this section the original SWMM model is briefly summarized. This model is based on a 1D simplification of the Saint-Venant equations corresponding to the mass and momentum continuity equations of free-surface flows: (3) $\frac{\mathrm{\partial }A}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{Q}_w}{\mathrm{\partial }x}=0$(3) (4) $\frac{\mathrm{\partial }{Q}_w}{\mathrm{\partial }t}+\frac{\mathrm{\partial }\left({Q_w}^2/A\right)}{\mathrm{\partial }x}+\mathrm{gA}\frac{\mathrm{\partial }H}{\mathrm{\partial }x}+\mathrm{gA}{S}_f=0$(4) where $A$ is the flow cross-section area, $Q_w$ is the water flow rate, $g$ is the gravitational acceleration, $x$ is the length and $S_f$ is the friction slope.

When the pipe gets pressurized, the model’s assumption of a free surface flow is no longer valid, and a surcharge method is required to continue the iterative process. Thus, the user can resort to one of two surcharge methods: Extended Transport (EXTRAN) or SLOT. Only the EXTRAN method is used herein, since the SLOT method did not show good results when applying an air accumulator in SWMM (Ferreira et al., 2023). Further information on the general SWMM engine and its numerical implementation can be found in Rossman (2017) and more details on the EXTRAN surcharge method can be found in Roesner et al. (1988). SWMM software version v5.1.015 is used herein.

3.2 AirSWMM

The above original SWMM model is modified by adding post-processing calculations that are performed at each time step, to account for different aspects of air pockets formation and fate during the simulation. Figure 7 presents a flowchart of these additional calculations composed of four main steps. Step 1 corresponds to the tracking and quantification of the initial air pocket volume at the moment the pipe filling is completed (i.e. waterfront reaches the downstream pipe end). In Step 2, air release and accumulator calculations are carried out to correct the air pocket volume according to the surrounding pressure. This step changes the hydraulic conditions (pressure and water flow rate) that will influence the air pocket creation and its initial volume. Step 3 consists of the verification of whether the air pocket is being dragged and if any air entrainment is occurring. In Step 4, it is identified if any air volume has arrived at the air pocket section and, if it has, both volumes are merged.

Figure 7. Flowchart of the proposed air pocket creation, transport and entrainment methodology.

This methodology builds upon the one presented in Ferreira et al. (2023) by modifying Step 1 to a more robust method that simulates waterfronts not perpendicular to the pipe axis and by adding Steps 3 and 4.

The final entrapped air pocket volume at time t is estimated based on the balance between the initial air pocket volume calculated in Steps 1 and 2, the air volume that is lost by the entrained air volume in Step 3 and the incoming bubbly flow to the air pocket from Step 4, as follows:

(5) $V_{\mathrm{AP}\text{,}f}=V_{\mathrm{AP}\text{,0}}-\sum _{t=0}^{t}a(F-1{)}^b\ast Q_{w,t}+\sum _{t=0}^t\frac{V_{\mathrm{AP}\text{,}\mathrm{drag}}}{V_p}{Q}_{w\text{,}t}$(5)

where $V_{\mathrm{AP}\text{,}f}$ and $V_{\mathrm{AP}\text{,0}}$ are the final and initial entrapped air pocket volumes, respectively, $a$ and $b$ are entrainment function parameters, $F$ is the Froude number, $V_{\mathrm{AP}\text{,}\mathrm{drag}}/V_p$ is the air-pipe volume ratio from the small air pocket creation from Figure 6. Details of the calculations in each of the four steps are detailed in the following sections.

3.2.1 Step 1: air pocket volume tracking and quantification method

Step 1 aims at identifying which pipes may contain entrapped air pockets and quantifying their volume. The method is initialized by going through each pipe to check if the water depths of its upstream or downstream nodes are lower than the pipe diameter (i.e. if water has with a free surface flow). Pipes with free surface flow are added to a “pool” of non-pressurized pipes. Pressurized pipes are assumed to continue to be pressurized in the subsequent time steps (pipe filling up) and, thus, are no longer required to be checked. The method proceeds to iteratively connect identified non-pressurized pipes and to add them to a “pool” of connected pipes. Each “pool” corresponds to an air pocket. When no pipes remain to be added to the “pool”, the air pocket tracking process is interrupted, and air release conditions are assessed for each air pocket. If any pipe pool is connected to an orifice, that air pocket features air release. Alternatively, an air pocket is considered an entrapped air pocket, if the pipe-set is between two pressurized nodes and not in contact with a non-pressurized orifice. The above process is repeated until no pipes are left in the “pool”, from where the tracking finishes and Step 2 starts.

The formation of entrapped air pockets occurs as follows (Figure 8): (a) the pipe fills with a free surface flow; (b) when a hydraulic jump is identified due to a slope change or an obstacle, an air pocket is created, and the jump creates a void upstream of the slope change; (c) the “void” pipe section between those nodes (marked using a different colour) corresponds to an entrapped air pocket. Once an air pocket is formed, several variables are initialized: the air pocket’s centre of mass, the water depth in the pipes where the air pocket is contained (pipe diameter minus the average depth the air pocket occupies in the pipe), the air pocket pressure and density, and the entrapped air pocket volume. The entrapped air pocket volume is computed by running the average water depth in the air pocket’s pipe using linear interpolation between nodes and discounting it to the total volume of where the air pocket is contained. This is schematized in Figure 8d.

Figure 8. Air pocket creation conceptual representation: (a) pipe filling with free surface flow, (b) sudden pressurization of the pipe with an empty volume in the sloped pipe, (c) filled pipe with entrapped air pocket and (d) numerical implementation of entrapped air pockets.

Three relevant assumptions are: (i) entrapped air pockets can increase and reduce volume over time; (ii) an entrapped air pocket has the volume equally distributed between the pipes where the air pocket is identified; (iii) an air pocket can move between pipes but keep its initial shape. An important remark is that air pocket location and volume are obtained simply by using flow rate and water depths that the original SWMM already calculates.

3.2.2 Step 2: air pressurization, release and coupling

3.2.2.1 Step 2.1: air release model

This step calculates the air release from each orifice from the system. Air release depends on the downstream boundary conditions and the air pressure inside the pipes. The air inside the pipes is initialized at atmospheric pressure $p_{\mathrm{atm}}$. Once the valve is opened, the pipe filling starts and the air inside the pipe compresses depending on the system boundary conditions. When the air pressure of the air pocket inside the pipe $p_{\mathrm{AP}}$ is such that $p_{\mathrm{AP}}/p_{\mathrm{atm}}<1.89$, the air release from an orifice occurs under subsonic conditions (Binder, 1955). The airflow rate is then described as follows (Zhou et al., 2002):

(6) $Q_{\mathrm{AP}\text{,}t}=\sum _{j=0}^{A{P}_{\mathrm{exits}}}{C}_{d,j}{A}_{0,j}{Y}_{j,t}\sqrt{2g\frac{{\rho }_{w,t-\mathrm{\Delta t}}}{{\rho }_{\mathrm{AP},t-\mathrm{\Delta t}}}\frac{(p_{\mathrm{AP},t-\mathrm{\Delta t}}-p_{\mathrm{atm}})}{{\gamma }_{w,t-\mathrm{\Delta t}}}}$(6) (7) $Y_{j,t}={\left[\frac{k}{k-1}{\left(\frac{p_{\mathrm{atm}}}{p_{\mathrm{AP},t-\mathrm{\Delta t}}}\right)}^{2/k}\frac{1-{\left(p_{\mathrm{atm}}/{\rho }_{\mathrm{AP},t-\mathrm{\Delta t}}\right)}^{(k-1/k)}}{1-\left(p_{\mathrm{atm}}/{\rho }_{\mathrm{AP},t-\mathrm{\Delta t}}\right)}\right]}^{1/2}$(7) where $Y$ is the air expansion factor after the orifice (Martin, 1976), $C_d$ is the discharge coefficient, $A_0$ is the cross-section area of the orifice, $g$ is the gravitational acceleration, ${\rho }_{\mathrm{AP}}$ is the air density, ${\rho }_w$ is the water density, $j$ is the orifice index from which the air is released and $k$ is the polytropic coefficient.

Conversely, whenever $p_{\mathrm{AP}}/p_{\mathrm{atm}}\ge 1.89$, the flow through the orifice becomes supersonic and the flow becomes choked with a maximum airflow rate being released as follows (Binder, 1955): (8) $\begin{array}{rl}{Q}_{\mathrm{AP}\text{,}t}& =\sum _{j=0}^{A{P}_{\mathrm{exits}}}{C}_{d,j}{A}_{0,j}\sqrt{g\frac{{\rho }_{w,t-\mathrm{\Delta t}}}{{\rho }_{\mathrm{AP},t-\mathrm{\Delta t}}}{p}_{\mathrm{AP},t-\mathrm{\Delta t}}}\\ & \times \sqrt{k{\left(\frac{2}{k+1}\right)}^{(k+1)/(k-1)}}\end{array}$(8) The air pressure $p_{\mathrm{AP}}$ is obtained from the air accumulator model in the next subsection.

3.2.2.2 Step 2.2: air accumulator model

This subsection complements the air release model by calculating the air pressure and density of the air mass downstream of the waterfront. Air masses in between a waterfront and any kind of boundary or another waterfront are described by the following equations assuming the air behaves as an ideal gas (Vasconcelos & Leite, 2012): (9) $\begin{array}{rl}\frac{d{p}_{\mathrm{AP},t}}{\mathrm{dt}}& =\frac{p_{\mathrm{AP},t\text{- }\mathrm{\Delta t}}k}{V_{\mathrm{AP},t}}\left(\sum Q_{\text{w},t-\mathrm{\Delta t}}-\sum Q_{\mathrm{AP},t}\right)\end{array}$(9) (10) $\begin{array}{rl}\frac{d{\rho }_{\mathrm{AP},t}}{\mathrm{dt}}& =\frac{{\rho }_{\mathrm{AP},t\text{- }\mathrm{\Delta t}}k}{V_{\mathrm{AP},t}}\left(\sum Q_{\text{w},t-\mathrm{\Delta t}}-\sum Q_{\mathrm{AP},t}\right)\end{array}$(10) (11) $\begin{array}{rl}{p}_{\mathrm{AP},t}& =p_{\mathrm{AP},t-\mathrm{\Delta t}}+d{p}_{\mathrm{AP},t}\mathrm{\Delta t}\end{array}$(11) (12) $\begin{array}{rl}{\rho }_{\mathrm{AP},t}& ={\rho }_{\mathrm{AP},t-\mathrm{\Delta t}}+d{\rho }_{\mathrm{AP},t}\mathrm{\Delta t}\end{array}$(12) where $V_{\mathrm{AP}}$ is the air pocket volume, $Q_w$ is the water flow rate filling and compression the air pocket, $Q_{\mathrm{AP}}$ is the released air flow rate obtained according to the equations in the previous subsection.

Air pockets at dead ends are described by Equations (6–12) with $Q_{\mathrm{AP}}$ being zero. Thus, the air pocket volume and density of dead ends and entrapped air pockets will vary over time but follow the ideal gas law for assumed near constant ambient temperature: $p_t(V_t{)}^k=p_{t+\mathrm{\Delta t}}(V_{t+\mathrm{\Delta t}}{)}^k$. A polytropic coefficient of $k$ = 1.2 is used here since the polytropic coefficient does not considerably influence pipe filling when there is air release (Ferreira et al., 2023).

3.2.2.3 Step 2.3: air–water interaction in flow rate calculations

This step presents the coupling between the air release and accumulator models from the previous subsections and the original SWMM. The original SWMM model calculates the water flow rate in each pipe as follows: (13) $Q_{w,t+\mathrm{\Delta t}}=\frac{Q_{w,t}+\mathrm{\Delta }{Q}_{w,t}^{\mathrm{inertia}}+\mathrm{\Delta }{Q}_{w,t}^{\text{p}\mathrm{ressure}}}{1+\mathrm{\Delta }{Q}_{w,t}^{\mathrm{friction}}}$(13) where $Q_{w,t}$ is the water flow rate in a specific pipe p at a given time step, $Q_{w,t}^{\mathrm{inertia}}$ is the flow rate change during the analysed time step corresponding to inertial forces, $Q_{w,t}^{\mathrm{pressure}}$ is the corresponding water flow rate change based on pressure forces and $Q_{w,t}^{\mathrm{friction}}$ is the water flow rate change based on the friction forces. More information on each of these terms can be found in Rossman (2017).

Two changes to above flow variation terms are required to account for the air compression during pipe filling events and entrapped air pockets. Firstly, the flow rate variation associated with the pressure component needs to be adjusted as proposed by Ferreira et al. (2023): (14) $\begin{array}{rl}\mathrm{\Delta }{Q}_{w,t}^{\mathrm{pressure}}& =-g\overline{A}\frac{H_2-H_1}{L}\mathrm{\Delta t}\text{to }\\ \mathrm{\Delta }{Q}_{w,t}^{\mathrm{pressure}}& =-g\overline{A}\frac{(H_2+H_{\mathrm{AP}})-(H_1+H_{\mathrm{AP}})}{L}\mathrm{\Delta t}\end{array}$(14) where $\overline{A}$ is the change in average flow area between over the analysed time step and $H_{\mathrm{AP}}$ is the pressure-head of the air pocket in gauge pressures. $H_{\mathrm{AP}}$ is null when $H$ is higher than the pipe diameter. Secondly, the original SWMM model calculates the flow rate considering the total water flow cross-section area when using the EXTRAN surcharging method. By introducing the entrapped air pockets and bubbly flow in the model, the pipe flow cross section area must be reduced by the cross-section area occupied by the air volume in the pipe where the flow rate is being calculated. Thus, the flow cross-section area becomes $A_p=A_{p\text{,0}}-V_{\mathrm{AP}}/L$. Flow rate inertial and friction terms in Equation (14) remain unchanged (other than the influence of reduced flow cross section area).

3.2.3 Step 3: air pocket dynamics

3.2.3.1 Step 3.1: air pocket drag model

A minimum critical flow velocity is required to overcome the surface tension near the pipe wall to induce movement to a static air pocket in the pipe. Previous literature contributions propose different values or ranges of critical flow velocities depending on the pipe’s slope and diameter. Theoretical approaches determined the critical velocity for horizontal pipes (Benjamin, 1968; Davies & Taylor, 1950; Dumitrescu, 1943) and experimental approaches of the critical velocity determined the ranges if critical velocity for varying pipe slopes (Goldring, 1979; Thomas M. Liou & Hunt, 1996; Walski et al., 1994). Gandenberger (1957) also proposed different critical flow velocities that are based on the air pocket volume. Most studies present formulations to obtain the air pocket critical flow velocity $U_{\mathrm{AP}\text{,}c}$, as in Escarameia (2004), the formulation that is used herein: (15) $U_{\mathrm{AP}\text{,}c}=(0.55\sqrt{\sin \theta }+0.53)\sqrt{g{D}_P}$(15) where $\theta$ is the pipe slope. Once the flow velocity at the upstream pipe of the air pocket is higher than the critical flow velocity, the corresponding air pocket starts moving at a given velocity.

Conversely, little information is available in the literature about air pocket velocity. Experimental data of Escarameia (2004) show the air pocket velocity increases with the critical flow velocity, but this is very limited for pipe angles higher than 6° and it does not consider the water flow velocity as a variable. Given this lack of experimental data, a model based on air pockets’ drag coefficients, the water flow velocity, the pipe diameter and the pipe slope are used here instead. Once critical flow velocity shown in Equation (15) is reached, the relative velocity of the air pocket, $U_{\mathrm{AP}\text{,}r}$ is calculated by Archimedes’ law accounting for buoyancy as follows: (16) $U_{\mathrm{AP}\text{,}r}=\sqrt{\frac{4}{3}\frac{g{D}_P\sin \theta }{C_{\mathrm{drag}}}}$(16) where $C_{\mathrm{drag}}$ is the drag coefficient of the air pocket. A spherical shape is assumed for co-current air pockets when they are contained in two pipes, thus originating a $C_{\mathrm{drag}}$ = 0.47 (Idel'čik & Steinberg, 2005). Once the air pocket relative velocity is calculated, the final air pocket velocity, $U_{\mathrm{AP}}$, is obtained as: (17) $U_{\mathrm{AP}}=U_w-U_{\mathrm{AP}\text{,}r}$(17) where $U_w$ is the mean flow velocity.

When using the above equations, it is assumed that air pocket velocity is null until the critical flow rate is reached from Equation (15). In addition, it is assumed that when the air pocket centre of mass exceeds the pipe boundaries, the air pocket moves to the next pipe or set of pipes and the air pocket centre of mass is reset to zero.

3.2.3.2 Step 3.2: air entrainment model

Air entrainment from entrapped air pockets is a complex two-phase flow behaviour. The water flow has enough momentum to emulsify part of the air pocket in the form of air bubbles at its tail but not enough to fully drag the air pocket. Several literature contributions were made in this direction based on different experimental set-ups. In all these approaches, no air entrainment is observed for Froude numbers below 1 (i.e. for sub-critical flow) and the higher the Froude number in supercritical flow the higher the air/water flow ratio entrained by the water flow. Previous studies on free-surface flow set-ups (Kent, 1952; Rajaratnam, 1967; USACE, 1980) and pressurized flows (Escarameia, 2007; Kalinske & Bliss, 1943; Mortensen et al., 2011; Rabben et al., 1983; Wisner et al., 1975) agree that the air entrainment ratio, $Q_{\mathrm{AP}\text{,}\mathrm{ent}}$, can be obtained as follows: (18) $\frac{Q_{\mathrm{AP}\text{,}\mathrm{ent}}}{Q_w}=a(F-1{)}^b$(18) However, authors disagree on the $a$ and $b$ values, originating a considerable uncertainty range (Schulz et al., 2020). For this reason, these parameters will be calibrated in the next section to estimate the entrapped air pocket volumes in this work.

Further considerations are required for this model’s implementation. The air entrainment starts as soon as an entrapped air pocket is created, and no entrainment is assumed if the Froude number is below the unit. Air volume due to air entrainment is not included in the model as an entrapped air pocket but as a bubbly flow. This is done because they do not have a clear boundary at the pipe crown but rather travel within the flow as a mixture. These two types of two-phase flow are tracked and analysed separately in our model because they move at different velocities (although these two are still subjected to volume changes according to the ideal gas law $\mathrm{pV}=\mathrm{kRT}$). The assumption presented in the previous section on entrapped air pocket drag is applied to the bubbly flow as well: a centre of mass is assigned to this bubbly flow that travels inside the pipe and moves from pipe to pipe according to its centre of mass.

3.3 Simulating entrapped air pocket stochastic nature

The previous section presented the deterministic model to simulate the air pocket entrapment. This section presents how the stochastic nature of the air pocket entrapment observed in the experimental tests is introduced.

A single final air pocket is obtained by running a single simulation with Steps 1–4 because the model is deterministic. The stochastic nature of the phenomenon is introduced by varying the air entrainment coefficients from Steps 3.1 and 3.2, obtaining a range of entrapped air pocket volume after a predetermined amount of model runs rather than a single volume. Thus, the user should specify the number of simulations, run the model each time with a different combination of $a$ or $b$ values (according to with a predetermined distribution and interval) and obtain the range of volumes.

It is recommended to only modify only one of these two parameters in simulations within the defined interval, since different pairs of $a$ and $b$ might lead to the same final entrapped air pocket volumes. Varying both parameters at the same time would only increases the required computational time and result in the same air pocket volume range.

4 Results and discussion

4.1 AirSWMM calibration

The calibration of the newly proposed model is done in two stages. First, a spatial discretization analysis is carried out to assess the influence of pipe length on the entrapped air pocket volume throughout the pipe. Hydraulic head results are then compared between the proposed methodology with the predetermined spatial discretization and results from the previous air model incorporation in SWMM (Ferreira et al., 2023) to determine if any major changes are observed. Secondly, entrapped air pocket volumes are calibrated adjusting the air entrainment parameters in Equation (18). Hydraulic head and air pocket volume from tests with orifices d = 2.2, 3.0 and 21 mm are used as calibration data for all the above stages. Data from d = 4.5 mm are used only for validation.

4.1.1 Spatial discretization analysis

The analysed pipeline system implemented in the experimental rig is modelled using the modified SWMM model proposed in this work. The “surcharge depth” of every node is set to 100 m for the nodes not to pond. The required spatial discretization is obtained by a mesh analysis. The pipe length is progressively decreasing and the corresponding time step ratio is obtained for pipe filling events according to $\mathrm{\Delta t}=0.1L/\sqrt{g{D}_P}$ using the EXTRAN surcharge method (Vasconcelos et al., 2018). Simulations are run with no air dragging or air entrainment (thus only following steps 1 and 2 of the methodology) to obtain the initial air pocket volumes in the pipe from Equation (4). The initial air pocket volumes are obtained for d = 2.2, 3.0 and 21 mm and for different pipe lengths (Figure 9).

Figure 9. Initial air pocket volume for different normalized pipe lengths.

This analysis has shown that only pipe lengths $L/D_P$ between 1.7 and 2.3 (see grey rectangle in Figure 5a) return the air pocket volumes within the same order of magnitude of those experimentally observed. The average value of $L/D_P$ = 2 is considered for further simulations. This length-diameter ratio corresponds to 300 pipes in the numerical model, each with $L$ = 0.042 m. The time step is obtained using the following equation $\mathrm{\Delta t}=0.1L/\sqrt{g{D}_P}$ = 0.00924 s, proposed by Vasconcelos et al. (2018).

Experimental hydraulic head for a sampled test from d = 3.0 mm is compared with the results obtained by the new proposed model and by the previous model from Ferreira et al. (2023) (Figure 10). The new model provides a better fit with experimental data than that from Ferreira et al. (2023), since the model no longer relies on the assumption of perpendicular waterfront to the pipe which was inherent in the piston equation previously used Ferreira et al. (2023). A lower air pressure is obtained during the filling process because the air volume is better quantified in the proposed version of the model (Equations 9 and 10). This allows for a better estimate of the arrival of the waterfront to the downstream end time because the downstream node of the system can pressurize sooner. The presented model still identifies the existence of a pressure peak but is not capable of describing it due to the limitations presented in the Discussion section. The previously developed model by Ferreira et al. (2023) used a piston equation to simulate the pipe-filling behaviour with air pressurization. However, the assumption that the waterfront is perpendicular to the pipe axis is very restricted, since most pipe systems are undulated. The model presented herein does not require such an assumption, making the model more robust and applicable to rising, horizontal and descending pipes. Thus, this new model is more accurate than the previously proposed one with the exception of continuously rising pipes with a considerable slope for which results from both models are equivalent. This is because, in such conditions, the waterfront's tail length (i.e. the length of water further than the corresponding pipe axis location) is negligible and can be considered perpendicular to the pipe axis. This is numerically demonstrated by the obtained root mean square error (RMSE) values which are 0.0395 m for Ferreira et al. (2023) and 0.0158 m for the new model. The main drawback of the new model is requiring a more detailed spatial discretization resulting in increased computational time (i.e. 10 times slower).

Figure 10. Experimental data and numerical results hydraulic head time series for d = 3.0 mm from the proposed methodology and Ferreira et al. (2023) model at pressure transducers PT1, PT2 and PT3.

4.1.2 Air pocket volumes

Entrainment factors, $a$, are calibrated to obtain maximum, average and minimum entrapped air pocket volumes by running the deterministic model 200 times. Tested entrainment factor values varied between 0 and 30 (i.e. parameter $a$ of Equation 18). Thus, a predicted entrapped air pocket volume range is obtained and should be compared to that experimentally obtained. Entrainment rate parameter $b$ is considered constant, equal to 1.3 (i.e. the average value between different contributions in literature), in all simulations. A total of 200 model runs is carried out because no air pocket volume change over 100 mm³ is observed in each quartile of predicted air volumes. The same number of runs is used in the validation process.

The values of $a$ that lead to absolute errors of air volumes higher than 100 mm³ are discarded. By determining the range of the entrainment factors $a$ for each calibration diameter, regression laws are calculated to estimate values of $a$ for validation. The obtained calibrated parameter values are shown in Figure 11. The minimum value of parameter $a$ obtained for d = 2.2 mm is not used for calibration purposes since the model is not able to reproduce the initial entrapped volume and, thus, this is considered as an outlier. This is a limitation of the model and will be discussed in the Discussion section. The calibration datasets are expressed in the cross-section area of the orifice with the pipe, $s/S$. Figure 11 shows the higher the orifice size ratio, the lower is the entrainment factor range. Also note that values higher than observed in literature parameter $a$ values are required to attain the observed volumes (0.03 in literature and 0.05 here).

Figure 11. Calibration curves for maximum, average and minimum entrainment factors for d = 2.2, 3.0 and 21 mm cross section area ratios s/S.

Orifices d = 2.2 mm and d = 3.0 mm were used to obtain the sharper flow rate variations and their influence on the air entrainment. Orifice d = 21 mm was used as the opposite extreme to obtain the minimum $a$ air entrainment coefficients. These orifices allow covering the full range of flow rates possible to be analysed in this system. Orifice d = 4.5 mm, which is an intermediate one in terms of size given the range used here, is used for validation using an entrainment coefficient according to the fitted laws. More diameters could have been tested and used to calibrate and validate the model. However, as seen in Figure 5a, the air pressurization effect on the entrapped air pocket decreases for larger orifice diameters and, ultimately (for orifices > 10 mm), tends to the air entrainment behaviour for d = 21 mm. The reader can also observe that the air pressure variation during the filling is progressively decreasing with the orifice size increase in Figure 3d. There are three main reasons for the observed discrepancies in the values of parameters $a$ and $b$. Firstly, as shown in literature, $a$ and $b$ values strongly depend on the experimental set-up size, configuration and being pressurized or free surface. Even though the highest $a$ value found is 0.03, the authors mention that downstream boundary conditions are relevant. Higher entrainment factor values are necessary for d = 2.2 mm, where the air release is severely constrained, and the air cushioning effect actively delays the pipe filling. Secondly, the air pressure might influence the entrainment rate, a phenomenon not accounted for in previous studies (Pothof & Clemens, 2010; Pozos et al., 2010). Literature experiments were carried out under steady state flows and with air being injected artificially and not entrapped by hydraulic means. More experimental research is required to validate this hypothesis. Thirdly, incorporating a turbulent and complex 3D phenomenon into a 1D model carries uncertainties and might not be able to fully reproduce behaviours observed experimentally.

4.2 AirSWMM validation

The entrapped air pocket volumes are validated using a different dataset for the orifice size of d = 4.5 mm. The fixed value of $b$ = 1.3 is used in all 200 model runs with uncertain parameter a value represented using a triangular probability density function with the same minimum to the maximum $a$ values and with the mean value of the distribution being the entrainment factor that corresponds to the average air volume.

The comparison of experimental air pocket distributions (based on 20 repeated experiments) with the corresponding distributions of model predicted air pocket volumes (based on the 200 simulation runs) are shown in Figure 12. The respective maximum, average and minimum air pocket volumes are presented in Table 1 together with absolute and relative errors.

Figure 12. Comparison between experimental and predicted entrapped air pocket volumes for calibration and validation.

Table 1. Maximum, average and minimum experimental and numerical entrapped air pocket volumes and their respective errors for all orifice sizes.

		Calibration			Validation
d (mm)		2.2	3.0	21	4.5
Experimental air pocket volume (mm^3)	Min	2835	1832	0	422
	Average	24,270	6773	280	6939
	Max	49,868	15,455	732	19,207
Predicted air pocket volume (mm^3)	Min	2893	1820	0	0
	Average	24,278	6775	350	7739
	Max	36,810	15,495	700	13,673
Relative error (%)	Min	−2.0	0.7	0.0	–
	Average	0.0	0.0	−25.0	−11.5
	Max	26.2	−0.3	4.4	28.8
Absolute error (mm^3)	Min	−58	12	0	422
	Average	−8	−2	−70	−800
	Max	13,058	−40	32	5534

Small differences are observed between the experimental and predicted air pocket volumes for three calibration cases. The only exception is the case of orifice size of d = 2.2 mm. The original SWMM engine is not able to reproduce the maximum volumes for d = 2.2 mm. Nevertheless, these limitations are only observed for the most extreme scenario and constrained air release, which is not expected to be observed in water networks. Additionally, high relative errors are observed for minimum air volumes, though these are not as relevant as average and maximum because of safety purposes. The remaining air pocket volumes for calibration orifice sizes show good agreements between predicted and observed, i.e. small errors are obtained for maximum, average and minimum volumes.

The results obtained for the validation dataset of d = 4.5 mm present relatively small errors for minimum and average air pocket volumes (11–29%). Some larger errors are obtained for the maximum air pocket volume which can be caused by higher uncertainties from: (i) the influence of the flow cross section reduction for the flow rate calculation in the air pocket cross section, which can be much more complex than a simple and linear interface as needed to consider in a 1D model; and (ii) the air entrainment being considered Froude dominated which is a simplification from the observed by Schulz et al. (2020), who observed the pressure at the air pocket location slightly influences the air entrainment.

Thus, the modified SWMM model incorporating air detection, location, dragging and entrainment allows making reasonably accurate predictions of entrapped air pocket volumes in a system analysed here, despite the fact that complex related phenomena are modelled in simplified way based on a one-dimensional modelling approach.

5 Discussion

5.1 AirSWMM limitations

The proposed AirSWMM model has limitations originating from two different sources: the SWMM engine and the new AirSWMM model. There are two limitations associated with the original SWMM model (i.e. SWMM numerical engine). The first is that SWMM cannot simulate sub-atmospheric pressures. SWMM assumes a free-surface flow and uses the Saint-Venant equations to solve the flow under such a regime. When the flow is pressurized, it either uses the EXTRAN surcharge method (assuming fully pressurized flow, solving the flow through another set of equations) or uses the Preissmann SLOT method (that uses an artificial slot to solve the Saint-Venant equations). However, neither of these can simulate sub-atmospheric pressures which other pressurized flow models can. Secondly, intermittent water supply systems are also susceptible to hydraulic transients (Erickson et al., 2022) which cannot be correctly reproduced in the SWMM. The EXTRAN surcharge method considers a wave celerity equals to $\sqrt{g{D}_P}=0.45$ m s^–1 and the SLOT method considers a celerity varying with the slot width, $B$, equal to $\sqrt{\mathrm{gS}/B}=12.59$ m s^–1. None of these formulations reproduce a realistic pipe wave celerity (around 300 m s⁻¹ for acrylic pipes) obtained by $c=\sqrt{\left(K/{\rho }_w\right)/c_1(1+[\left(K/E\right)\left(D_P/e\right)\left]\right)}$, the pipe wave celerity, in which $K$ is the water bulk modulus, $c_1$ is a constant dependent on the pipe support conditions and the pipe material, $E$ is the Young modulus of elasticity of the pipe and $e$ is the pipe wall thickness. Water compressibility would need to be several times lower for the celerities in the model to be representative of reality To the authors’ knowledge, only the SLOT surcharge method can describe the elastic column behaviour by changing the slot width as proposed by Pachaly et al. (2021), making slot width equal to $B=\mathrm{gS}/c^2$. However, this modification has only been verified in conceptual conditions and is still to be compared with measurements.

Additional limitations come from the surcharge method and air pocket shape used in the implementation of AirSWMM. First, this model only provides good results when using the EXTRAN surcharge method, as observed in Ferreira et al. (2023). This limits the AirSWMM model to this surcharge method even though the SLOT method provides better results for fast unsteady events as demonstrated by Pachaly et al. (2019). Secondly, the air pocket shape is imprecise because the 1D nature of SWMM does not allow the introduction of the round shape of air pockets, caused by the surface tension equilibrium laterally and longitudinally, and therefore inhibits the correct air-pocket length estimation in horizontal pipes. As a consequence, no additional local head losses are taken into account where air pockets are estimated due to the shape’s imprecision. The only modification that reduces the flow rate around the air pocket is the flow cross-section area reduction corresponding to the existing air pocket depth.

5.2 AirSWMM applications and developments

The proposed model can be extended to incorporate urban drainage systems features, such as manholes or pipe shafts. That will require further developments, namely: (i) including the air volume in storage components to the overall air being pressurized; (ii) including air release devices possibly connected to such components; and (iii) implementing a local air movement model so that air in the pipes would rise in the manhole or shaft rather than being dragged in the pipes. Another element that could be included is the air valve which can also be found in urban water networks.

6 Conclusions

A methodology for simulating the creation, location and entrainment of air pockets during filling conditions in a pipeline system is proposed. Novel experimental pipe filling tests are carried out in a laboratory set-up to understand the expected pressure heads the entrapped air pocket formation at high points and to quantify their volumes. The new model, labelled AirSWMM, is implemented as an extension of the existing SWMM model. It is an upgrade of the previously published model (Ferreira et al., 2023) allowing the waterfront tracking based on the hydraulics (SWMM) as well as locating and replicating the drag and entrainment of air pockets. The new model AirSWMM was calibrated and validated using the collected experimental data. The main conclusions are as follows:

• The new AirSWMM model captures reasonably well different aspects of air pocket creation and its fate during the pipe filling conditions. It is able to predict the location and final volume of an entrapped air pocket with an average relative error of 20%, which is deemed good given the complex nature of the analysed phenomena and the use of a 1D model. Still, some dynamic behaviours such as air-pocket interface disruption could not be simulated due to the complexities involved and limitations of a 1D model.

• The obtained experimental observations provided evidence of a stochastic nature of air pocket creation and its dependency to air release conditions. Since the proposed AirSWMM model is deterministic in nature, a set of simulations should be run with different entrainment rates to obtain a realistic range of entrapped air volumes. In this work 200 samples were used for this, resulting in the aforementioned prediction accuracy.

• The obtained experimental observations also provide evidence that lower filling flow rates (which are commonly recommended in practice) tend to create larger entrapped air pocket volumes whilst delaying the pipe filling time. Therefore, this advice should be revisited in future work.

The AirSWMM methodology presents considerable scientific advances as no previous 1D model allowed for location and correct quantification of entrapped air pockets. Still, the proposed model has some limitations. The air accumulator model assumes the air inside the pipe pressurizes all at once and does not take into consideration the compressibility rate of the air. This methodology no longer requires the piston equation to track the waterfront, thus allowing each air pocket to be analysed separately.

Additional pipe systems should be tested to further validate this methodology in larger set-ups to account for scale effects. A see-through and unburied set-up would be required to validate this methodology in a real system. This methodology has been implemented in SWMM but could be implemented in other free-surface flow models, provided a stable algorithm coupling between the air and the water phase is obtained. Further research should focus on testing this methodology in other linear and branched pipe layouts.

Notations

$a$, $b$	=	air entrainment parameters (–)
$A$	=	flow cross section area (m2)
$\overline{A}$	=	average flow area variation (m2)
$B$	=	Preissmann slot width (m)
$c$	=	pipe wave celerity (m s−1)
$c_1$	=	pipe constraint coefficient (–)
$C_d$	=	discharge coefficient (–)
$C_{\mathrm{drag}}$	=	drag coefficient (–)
$d$	=	orifice diameter (m)
$D$	=	pipe diameter (m)
$e$	=	pipe wall thickness (m)
$E$	=	pipe Young modulus of elasticity (Pa)
$F$	=	Froude number (–)
$g$	=	gravitational acceleration (m s−2)
$H$	=	pressure-head (m)
$k$	=	polytropic coefficient (–)
$K$	=	water bulk modulus (Pa)
$l$	=	longer axis air pocket length (m)
$L$	=	pipe length (m)
$p$	=	pressure (Pa)
$Q$	=	flow rate (m3 s−1)
$\mathrm{\Delta }{Q}_{t,p}^{\mathrm{inertia}}$,	=	SWMM water flow
$\mathrm{\Delta }{Q}_{t,p}^{\mathrm{pressure}}$	=	rate components
$\mathrm{\Delta }{Q}_{t,p}^{\mathrm{friction}}$	=
$r$	=	air pocket radius (m)
$R$	=	gas constant (J K−1 mol−1)
$R_e$	=	Reynolds number (–)
$s$	=	orifice cross section area (m2)
$S$	=	pipe cross section area (m2)
$S_f$	=	friction slope (–)
$t$	=	time (s)
$T$	=	temperature (K)
$T_s$	=	Preissmann slot width (m)
$U$	=	average velocity (m s−1)
$V$	=	volume (m3)
$Y$	=	orifice expansion factor (–)
$\mathrm{\Delta t}$	=	time step (s)
$\mathrm{\Delta x}$	=	space step (m)
$\rho$	=	specific density (kg m−3)
$\theta$	=	pipe slope (–)

Subscripts

0	=	initial
AP	=	air pocket
atm	=	atmospheric
c	=	critical
drag	=	dragged air component
j	=	orifice index
f	=	final
max, min	=	maximum, minimum
p	=	pipe
r	=	relative
w	=	water

Additional information

Funding

The authors acknowledge the Portuguese national funding agency, Fundação da Ciência e Tecnologia (FCT), for funding the reported work within grant number SFRH/BD/146709/2019 and for the support through funding UIDB/04625/2020 for the research unit CERIS; Civil Engineering Research and Innovation for Sustainability.