## Abstract

Modelling air movement in sewer networks is needed in order to address the issues related to sewer odour complaints and sewer corrosions due to hydrogen sulphide in sewers. Most of the existing air flow models can only be applied in small sewer networks or the trunk lines of sewer systems. The purpose of this paper is therefore to propose a theoretical approach to formulate a general governing equation set for modelling steady air movement in large sewer systems. This approach decomposes the sewer system of interest into its basic physical components as pipes and nodes, and builds local topology of each pipe and each node based on geographic information system data as the fundamentals of model formulation. It avoids manually identifying each branch of the sewer system, eliminates the effect of physically closed networks in sewer systems on the governing equations, and considers key sewer components and all known driving forces. The proposed approach was applied to a real sewer system with over 500 pipes. The results show that the proposed model is applicable in modelling air movement in a large sewer system and provides a general idea of sewer gases moving through the system and their emission.

## HIGHLIGHTS

The proposed framework is a generalized model for modelling air flow in a large real sewer system.

The proposed model can build up the local topology and the skeleton of the sewer system within the computational domain automatically from the geographic information system data of a sewer system.

The equation matrix for the proposed model solves its unknown variables, air flowrate and air pressure, simultaneously, which increases its computational flexibility.

The proposed model can be applied to a combined sewer system even under the circumstance that the system has physical loop(s).

## INTRODUCTION

Hydrogen sulphide (H_{2}S) can be generated under anaerobic conditions in sanitary sewers and transported by the air in the sewer headspace (Hvitved-Jacobsen *et al.* 2013). H_{2}S has been identified as the leading cause of sewer odour complaints (Edwini-Bonsu & Steffler 2004) and sewer corrosion (Adkins *et al.* 2019). Thus, a general air flow model is urgently needed to model air flow in a large sewer system to estimate the distribution of air flowrate. However, most of the existing models are only able to simulate the headspace air flow along the sewer trunk line or in a small sewer network (Edwini-Bonsu & Steffler 2006a; Brocard *et al.* 2012; Wang *et al.* 2012; Hentz *et al.* 2013; Lowe 2017; Qian *et al.* 2018).

Modelling on air flow in a sewer system was based on the mechanistic understanding of air flow dragged by wastewater in sewer pipes. Pescod & Price (1982), Monteith *et al.* (1997), Edwini-Bonsu & Steffler (2004, 2006b) and Lowe (2017) proposed empirical correlations for estimating the air flow velocity in a single pipe from the known wastewater velocity based on their experimental data or computational fluid dynamics simulation results. Together with a continuity equation, Lowe (2017) proposed a ventilation model by using an empirical relation to model air movement in a small branch sewer system which only includes sewer pipes, junctions between pipes and standard manholes. However, this model cannot simulate the air movements in specific sewer components such as drop structures, ventilators and pump stations, and the correlation relations do not include other driving forces such as pressure differential and friction.

To model air movement in a sewer system with specific sewer components, Edwini-Bonsu & Steffler (2006a) proposed a theoretical framework of a ventilation model. The model considered the driving forces of wastewater drag and pressure differential using a correlation relation, and estimated the pressurization of drop structures and ventilators with assumed second order polynomial functions. It first took the ambient atmosphere as a reference line, then made each branch in the sewer system to be an assumed closed loop with the related junctions, sewer pipes and the reference line. Based on the assumed loops, using zero energy loss around a loop and air mass continuity at each node, equation sets were then derived. However, the model was only applied to a sewer system with 45 pipes, 15 standard manholes and 9 drop manholes. Hentz *et al.* (2013) employed this model to compute air flowrates and pressures along a trunk line in the city of Los Angeles and concluded that this model could be used for planning and decision-making.

The model of Edwini-Bonsu & Steffler (2006a) may not be suitable for a large real sewer system. First, real sewer systems can have over hundreds or thousands of pipes, manholes, junctions, etc. Consequently, it could be quite difficult, even impossible, for users to manually identify each pipe and node and then assemble them into assumed closed networks with the presumed reference line. Secondly, the model did not consider the effect of the existing physically closed networks in a sewer system on the solving of the equation set. Thirdly, the effect of dropshafts on air pressurization in sewers (Ma *et al.* 2016) needs to be incorporated in the model. Additionally, for a large real sewer system, the topology of the sewer system needs to be automatically built based on the geographical information system (GIS) data as GIS has been applied in the design and management of urban sewer systems for decades (Abraham *et al.* 1998; Venigalla & Baik 2007).

To account for known driving forces on headspace air flow, Brocard *et al.* (2012), Wang *et al.* (2012) and Qian *et al.* (2018) developed ventilation models to simulate the air movement along a trunk line based on the balance of forces bearing on headspace air in a sewer pipe. Wastewater drag, pressure differential and pipe wall friction were included in their models. However, the model of Brocard *et al.* (2012) was Excel-based and only available for a simple sewer system without drop structures involved. Wang *et al.* (2012) proposed a full dynamic air flow model for unsteady air flow in gravity sewer systems but was only validated along the trunk line. Qian *et al.* (2018) mainly proposed an air leakage model to estimate the air flowrate via the cut-off laterals and evaluated the leakage effect of laterals on the air flow in the trunk. However, these models have not been applied to or validated in a large sewer system with over hundreds of pipes and various sewer components.

This paper is therefore to propose a general theoretical framework to model air movement in a large sewer system incorporating known driving forces and key sewer components. The framework formulates the governing equation set based on the basic physical components of the sewer system and the local topology at each node, and then solves the equation set with the built-in solver in MATLAB. Particularly, the skeleton of the sewer system of interest is built up for visual checking, the equation set is formulated by applying force balance in each pipe and mass continuity at each node, and the equation set has two variables (air pressure and air flowrate) which are solved simultaneously. Finally, this proposed theoretical framework is applied to a prototype sanitary sewer system in the city of Edmonton, Alberta, Canada and validated. The results show the model is applicable for modelling air flow in a large sewer system.

## FRAMEWORK ON AIR FLOW MODELLING

### Governing equations

In this paper, the air flow in a sewer system is assumed to be a one-dimensional incompressible steady flow. The temperature and the air density in the entire sewer system are assumed to be constant. The governing equations for air flow in a sewer system is the momentum equation in each pipe, and the continuity equation at each node.

*i*and is air flowrate into/out of a node; for a node without opening. The flowrate is assumed to be positive when air is out of a node and negative when air flows into a node; subscript

*i*is the index of pipes connected to this node.

*i*th pipe, of (

*i*+ 1)th pipe and at the bottom of the manhole equal the pressure at node . Based on the continuity equation, the air flowrate in the manhole shaft equals the air flowrate via the manhole cover opening. Thus, can be estimated using the orifice equation (Edwini-Bonsu & Steffler 2006a; Qian

*et al.*2018):where is a correction factor;

*C*is discharge coefficient through a pickhole;

_{o}*N*is number of pickhole opening; is area of each pickhole at a manhole cover; is the air density;

*g*is the gravitational acceleration;

*H*is the manhole depth; is the atmospheric pressure.

*A*is cross-sectional area of air flow; is cross-sectional area of wastewater flow;

*C*is drag coefficient on air by wastewater at the air–water interface in pipe; is wastewater flowrate in pipe;

_{D}*L*is pipe length;

*B*is width of air–water interface in pipe;

*f*is friction factor;

*χ*is air perimeter; is slope of pipe.

Force terms in Equation (3) vary in different situations. In all pipes without wastewater flow, the drag on air does not exist, i.e., , thus . If Equation (3) is applied to some special sewer structures, e.g., drop manholes and ventilators, the determination of in Equation (4) needs special attention. In this paper, in drop structures is estimated using the pressurization formula proposed by Ma *et al.* (2016), while across ventilators is estimated by their performance function as indicated in Edwini-Bonsu & Steffler (2006a).

*m*and

*n*are the nonlinear index of air pressure and air flowrate respectively.

The equation set Equation (7) differs from the one proposed by Edwini-Bonsu & Steffler (2006a). In this paper, it does not need to assemble the assumed closed networks; each governing equation is not necessary to be transferred to unify the unknown variables into *Q*; consequently the unknown variables for Equation (7) are a combination of and .

### Flow chart of air flow model

Based on the equation set Equation (7), the simulation of air flow in a sewer system was programmed following the procedure shown in Figure 2. The model developed by MATLAB reads the topology of sewer system in the form of GIS data, identifies the computational domain, classifies pipes and nodes using the types of facilities located at different GIS layers, and builds local topology at each node and each pipe. Then, the model establishes the equation matrix and initiates the unknown variables using boundary conditions and empirical equations. Finally, the model solves the nonlinear equation matrix by the solver built in MATLAB and outputs the parameters related to the air movement in a sewer system.

Based on the GIS data of the topology of the sewer system provided by the system owner/operator, a module was first programmed to decompose a sewer system into its basic sewer elements and to build local topology. The GIS data have several layers. Each layer includes one kind of facility (pipes, manhole, junction, pump station, ventilator, outfall, etc.) that the system owner/operator has classified. Thus, with the polygon surrounding domain defined by the GIS coordinates, the module can identify each pipe located within the domain based on the GIS coordinates of pipes, then number them in sequence as pipe ID and store the related pipe information withdrawn from GIS data to the structure array of pipes. Then, the module identifies upstream and downstream facilities of each pipe based on GIS data, removes the repeated ID of each facility to make it unique, numbers them in order as node ID and stores the related node information to the structure array of nodes. Finally, the module builds the local topology by identifying the pipe IDs connected to each node ID and the node IDs connected to each pipe ID, and saves them in the related structure arrays. Thus, the topology of the entire sewer system is readable and accessible by the module at any time.

The module can automatically build up the skeleton of the studied sewer network for the users to visually check if the selected computing domain is correct. Figure 3(a) and 3(b) show examples of the skeleton of a sewer system in Edmonton, Alberta, Canada and its trunk line, in which T1 to T10 are manhole IDs along the trunk line. Figure 3(c) is the local zoom-ins near manhole T1 and T5 shown in Figure 3(a) (in which P1 to P3 are pipe IDs), and shows the local topology of each pipe and each node. For example, node T1-2 has two upstream pipes (P1 and P2) and one downstream pipe (P3), while pipe P1 has one upstream end at manhole T1 and one downstream end at node T1-2. The local topology of each pipe and each node is fundamental to formulate Equation (7) for headspace air flow model in sewers in this study.

Equation (7) is a nonlinear algebraic set but has a unique solution as it is square. However, it has two unknown variables, air pressure and air flowrate, and can only be solved simultaneously. In MATLAB, two built-in solvers, *fsolve* and *lsqnonlin*, are capable of solving a large-scale nonlinear equation set. *fsolve* solves a system of equations using a quasi-Newton method and starting at *x _{0}* to find a vector

*x*that makes all function

*f*(

_{i}*x*) = 0 individually, while

*lsqnonlin*expects that function

*f(x)*is vector valued, and implicitly squares the values of

*f(x)*, then minimizes the sum of the squares

*.*For comparison, both solvers with all of their built-in algorithms were employed in this study.

*et al.*2012), or MIKE URBAN model (Edwini-Bonsu & Steffler 2006a), or field measured data saved in a file in a required format. The air pressure at each node can be simply initiated as 0, but the air flowrate in each pipe can be initiated using the explicit power fit of experimental data of air flow velocity of Pescod & Price (1982), which was fitted by WERF (2009), to speed up the iterations

Based on the initial values of air flowrates and air pressures, *fsolve* or *lsqnonline* can start its iterations.

## RESULTS AND DISCUSSION

The proposed air flow model in this paper was applied to a sanitary sewer system in Edmonton, Alberta, Canada for preliminary validation and to examine its capability for simulating air flow in large-scale sewer systems. Two series were designed for comparison. Their skeletons shown in Figure 3 were built up from the pipe attribute and manhole attribute in 2014 with GIS coordinates provided by the system owner. Series A includes all sewer components in the study sewer system. All pipes in the system are concrete pipes. The sewer trunk has a total length of 3.0 km with a diameter of 1,500 mm from T1 to T7 then reducing to 1,200 mm after T7. The laterals have a diameter varying from 200 to 375 mm and a length varying up to 179.8 m. The upstream end of its trunk line is a dead end with manhole T1 sitting above it. Manholes T2 to T5 are connected to the trunk line via short pipes similar to the one in Figure 3(c). Manholes T7 to T10 are sitting above the trunk line with lateral pipes connected to it. Manhole T6 is connected to the trunk line with a short pipe and also has laterals connected to it. The pump station is located at the downstream end of the trunk line. The trunk main between manhole T6 and T8 oscillates between full pipe water flow and partially full pipe flow in responding to pump operation (Guo *et al.* 2018; Qian *et al.* 2018). For the simulated cases in this paper, the trunk main after manhole T8 is considered to be full. Series B is a simplified system of series A and only includes the trunk mains with all laterals cut off. The details of the trunk line, the field measurement and the operation of the pump station can be found in Guo *et al.* (2018).

Table 1 lists a summary of the statistical data of sewer components of two series. The simulated time is at 9 p.m. on 7 July 2016. The pressure at manhole T1 was initiated with the field measured value. For all manholes and junctions with connected lateral pipes cut off, they were set to open to the atmosphere, which permits air to be inducted into the sewer or emitted from the sewer freely. The wastewater flow parameters in the study system were the simulation results of MIKE URBAN. The initial values of air pressures in Equation (7) were assigned to be 0, and the air flowrates were initiated using Equation (8).

Parameters . | Series A . | Series B . |
---|---|---|

Brief description | All trunks and laterals in the domain of interest | Only trunks |

Number of sewer pipes | 518 | 27 |

Number of drop structures | 6 | 4 |

Number of standard manholes | 413 | 6 |

Number of junctions | 72 | 15 |

Number of pump stations | 1 | 1 |

Boundary | Initial pressure at manhole T1 | Initial pressure at manhole T1 |

Parameters . | Series A . | Series B . |
---|---|---|

Brief description | All trunks and laterals in the domain of interest | Only trunks |

Number of sewer pipes | 518 | 27 |

Number of drop structures | 6 | 4 |

Number of standard manholes | 413 | 6 |

Number of junctions | 72 | 15 |

Number of pump stations | 1 | 1 |

Boundary | Initial pressure at manhole T1 | Initial pressure at manhole T1 |

Both series were solved by two solvers with their built-in algorithms in MATLAB: *fsolve* with algorithms of trust region (TR), trust region with dogleg (TRD) and Levenberg–Marquardt (L-M), and *lsqnonlin* with algorithms of trust region reflective (TRR) and L-M. All cases were run by a desktop computer with an Intel Core i7-7700 CPU @3.60 GHz and 16 GB memory. The version of MATLAB is R2018b for academic use. The parallel option was ‘on’. Table 2 lists the statistical data for solving and shows that both series can only be solved by the solvers of *fsolve* or *lsqnonlin* with the L-M algorithm. Two series solved by solvers with other algorithms, i.e. *fsolve* with algorithms of TR and TRD, and *lsqnonlin* with TRR, did not converge. The running processes quitted automatically without a solution once its iteration step size (i.e. StepTolerance) was less than 10^{−8}, where the StepTolerance in MATLAB is defined as a lower bound on the size of a step.

Series No. . | lsqnonlin. | fsolve. | |||
---|---|---|---|---|---|

L-M . | TRR . | L-M . | TR . | TRD . | |

Series A | 14,882 s/441 | N/A | 14,078 s/ 441 | N/A | N/A |

Series B | 161 s/118 | N/A | 165 s/118 | N/A | N/A |

Series No. . | lsqnonlin. | fsolve. | |||
---|---|---|---|---|---|

L-M . | TRR . | L-M . | TR . | TRD . | |

Series A | 14,882 s/441 | N/A | 14,078 s/ 441 | N/A | N/A |

Series B | 161 s/118 | N/A | 165 s/118 | N/A | N/A |

*Note:* N/A for not solved.

For the two converged cases solved with the L-M algorithm in the same series, their computing times are quite close to each other, e.g., 14,078 s for series A with solver *fsolve*, and 14,882 s for series A with solver *lsqnonlin*, and the simulated air flowrates and air pressures at the same locations are identical. Figure 4 shows the variation of the maximum absolute residual against the iteration number. Two curves overlay each other. The maximum absolute residuals of both cases increase from 3,452 to 4,077 at the first iteration, then drop sharply to 250 in 3 iterations, down to 2.47 at the 60th iteration, and the solvers finally converge at the 441st step with a total absolute residual less than 10^{−10}. It indicates that in this study, both built-in solvers, *fsolve* and *lsqnonlin* with the L-M algorithm, can be used to solve a large-scale nonlinear equation set for the steady air flow in a large sewer system with almost the same efficiency.

The model proposed in this paper has been preliminarily validated by comparing the simulated air pressures with the field measured data at the same time (Guo *et al.* 2018) and the field data measured at 3 a.m. presented in Qian *et al.* (2018) for tendency comparison as shown in Figure 5. The simulated pressures at manhole T5 and T6 for series A agree well with the measured data. The air pressures in manhole T2 to T5 in series A have similar magnitudes at around 570 Pa. It agrees well with the field measured result that the measured air pressure distribution along the trunk line between T2 and T5 has a similar magnitude, varying around 255 to 275 Pa at 3 a.m. as presented in Qian *et al.* (2018) and shown in Figure 5. However, a relative large discrepancy occurs at manhole T7 due to the air flow complexity around manhole T7. The trunk diameter suddenly decreases from 1,500 mm upstream of T7 to 1,200 mm downstream, and manhole T7 itself is a drill drop manhole. With the effects of both, the air movement at manhole T7 is complex and energy loss may be large across this node. The current model ignores these effects and thus might cause this large discrepancy.

For series B, the air pressure from manholes T2 to T5 decreases from 493 to 291 Pa quickly, and the pressure at manhole T5 is much smaller than the field measured value. This is because the lateral pipes connected to manhole T6 have been cut off, and manhole T6 is assumed to be open to the atmosphere, which permits air moving through freely without resistance. Therefore, to model the air flow in the trunk line, it is necessary to use a model to estimate the leakage air flowrate via the cut-off laterals as developed in Qian *et al.* (2018).

Figure 6 shows the general air movement in the study sanitary system of series A and B and the locations where air emits from the sewer. In the figure, for sewer pipes, the line width presents the magnitude of air flowrate, the color means the flow direction: green indicates air moving in the direction of water flow, red indicates the opposite, and black indicates no air flow. For nodes, the size of marker indicates the magnitudes of air flowrate, the red marker indicates the air leaking out of manholes from the sewer system, and the green for the air being inducted into manholes.

For both series, the trunk after manhole T8 runs full by water and no air can pass (black line) due to the existence of the pump station at the downstream end of the trunk. Therefore, the air in the trunk ahead of manhole T8 will be diverted from the trunks to laterals or be emitted to the atmosphere via manhole openings or cut-off laterals upon it arriving at manhole T8.

For series A, the air is inducted into the trunk via manhole T1-L1 at 2.11 m^{3}/s. Then the air moves downstream along the trunk line, and is emitted via the openings of manholes T2 to T8 and the laterals connected to manhole T6 to T8 due to the blocked headspace in the trunk main after manhole T8. The larger air emissions happen at manholes T2 to T4, T7 and T8-L1, at about 0.20 m^{3}/s. This gradual decrease of the air flowrate in the trunk from upstream to downstream shows a thinner and thinner green line in Figure 6(a). At manhole T8, the remaining air (0.43 m^{3}/s) in the trunk is forced to be emitted via the laterals connected to manhole T8.

Due to the cut-off of the laterals, the air pressure variation along the trunk line for series B is quite different from the one in series A as shown in Figure 5, as is the air flowrate as shown in Figure 6. For series B, pipe P1 is cut off from junction T1-2, from which air is inducted into the trunk at 5.35 m^{3}/s. This air flowrate is much larger than that of series A (2.11 m^{3}/s). Air is also emitted from the trunk via manholes T2 to T5 but at a lower flowrate than that of series A as its air pressure is lower than that of series A. However, the air flowrate released via manhole T6 for series B is quite large (2.22 m^{3}/s), about 10 times that of series A (0.14 m^{3}/s). This is because lateral pipes connected to manhole T6 in series B have been cut off and manhole T6 is assumed to be open to the atmosphere, while the laterals connected to manhole T6 in series A limit the movement of air moving through it. This difference between the simulated air flowrates of series A and series B indicates that a model to estimate the leakage flowrate of air via the cut-off laterals is necessary if air movement is only modelled in a trunk line, or in a sewer system where lateral pipes are cut off.

## CONCLUSIONS

A theoretic framework on air flow modelling was proposed in this paper to simulate the air flow in a large complex sewer system which can have hundreds or thousands of pipes, junctions and other structures, bypass pipes or connections to a storm system. A model has been developed based on the framework using MATLAB and applied to a sewer system in Edmonton, Alberta, Canada. The simulated results show that the proposed model is capable of simulating the air movement either along a trunk line or in a complete sewer network with over 500 pipes. Both solver *fsolve* and *lsqnonlin* with the L-M algorithm built in MATLAB are suitable for solving the large-scale nonlinear equation set in this study.

Compared to the previous models, the proposed model can automatically build up the local topology and the skeleton of an entire sewer system from GIS data with a user-defined computing domain, and solve the air flowrate and air pressure simultaneously, which makes the model have much more computational flexibility. It means the model can be applied to a combined sewer system even under the circumstance that it has physical loop(s). However, as the equation matrix has two unknown variables, air flowrate and air pressure, the solving of the matrix may encounter slow convergences or convergence difficulty. Also, the result may have larger discrepancy if the computed sewer system has some laterals that are cut off and set as free ends. Thus, further research is needed on the solving convergence, the leakage model for cut-off laterals and the further validation of the model.

## ACKNOWLEDGEMENTS

The authors gratefully acknowledge the financial support from Natural Sciences and Engineering Research Council (NSERC) of Canada, and the City of Edmonton/EPCOR Utilities Inc. Qi Zhang was sponsored by China Scholarship Council.

## DATA AVAILABILITY STATEMENT

Data cannot be made publicly available; readers should contact the corresponding author for details.