For the past 70 years, researchers have dealt with the investigation of odour in sewer systems caused by hydrogen sulphide formations and the development of approaches to describe it. The state-of-the-art models are one-dimensional. At the same time, flow and transport phenomena in sewers can be three-dimensional, for example the air flow velocities in circular pipes or flow velocities of water and air in the reach of drop structures. Within the past years, increasing computational capabilities enabled the development of more complex models. This paper uses a three-dimensional two-phase computational fluid dynamics model to describe mass transfer phenomena between the two phases: water and air. The solver has been extended to be capable of accounting account for temperature dependency, the influence of pH value and a conversion to describe simulated air phase concentrations as partial pressure. Its capabilities are being explored in different application examples and its advantages compared to existing models are demonstrated in a highly complex three-dimensional test case. The resulting interH2SFoam solver is a significant step in the direction of describing and analysing H2S emissions in sewers.
Wastewater in sewers undergoes a lot of physical and biochemical processes. One important factor is the formation of hydrogen sulphide (H2S), which can cause health risks for sewer workers. The tendency towards more complex and longer sewer networks can lead to longer retention times, which enhance the emission of H2S. Climate change at the same time causes higher temperatures in the wastewater, which increases emission rates.
In the past 70 years, extensive research has been performed to increase the knowledge on H2S formations and to develop approaches that describe the development of odour in sewers (e.g. Gilchrist 1953; Thistlethwayte 1972). The state-of-the-art models, which have been developed within the last 20 years, are horizontal one-dimensional. These are the SeweX model from Australia (Rootsey & Yuan 2010; Rootsey et al. 2012) and the WATS (wastewater aerobic/anaerobic transformations in sewers) model from Denmark (Hvitved-Jacobsen et al. 2013). Both are not public domain.
An overview of existing model approaches has been given in Carrera et al. (2017) and the need for further research has been highlighted.
To begin with, the mass transfer approach of the existing models is based on the so-called two-film theory, which uses different assumptions. The WATS model additionally uses different approaches to account for turbulent H2S transfer rates across the water surface in various applications. These different approaches are empirical or theoretical connections between oxygen and H2S transfer on the one hand and empirical models linking H2S emissions to flow properties in the pipe on the other hand (Carrera et al. 2017).
Wang et al. (2018) highlight the shortcomings of the two-film theory. It cannot account for local changes of the flow regime or variations of fluid properties. Furthermore, the theory is based on a constant liquid film that can change in real-life conditions due to flow instabilities. The most limiting factor, however, is assumed to be the one-dimensionality of the approach. More advanced approaches, the penetration theory and the surface renewal theory, can account for the variability of the flux over time but do not account for local variations, the change of fluid properties or flow regimes (Wang et al. 2018). This has already led to wide applications of computational fluid dynamics (CFD) models for mass transfer applications in the chemical industry (Wang et al. 2018).
Carrera et al. (2017) identified the models' shortcomings in describing mass transfer processes across the water surface, the current approaches of which were considered to be simplified, especially when considering hydraulic structures such as gravity sewers, junctions and water falls. Recent research on water falls or drop structures in sewer systems led to improved formulations to account for the effect of local turbulence (Matias et al. 2017a, 2017b), but these approaches are still empirical equations which are fed into the model.
This short overview leads to the question whether a three-dimensional simulation model could help in increasing the process understanding, especially when analysing complex and turbulent flows in a sewer. Another benefit could be the in-depth analysis and design optimization in hotspots of H2S emissions.
In order to address this question, a volume of fluid (VOF) approach as implemented in OpenFOAM's solver interFoam has been chosen to describe the two-phase flow of water and air. This solver has already been used for a number of demanding hydraulic applications (e.g. Thorenz & Strybny 2012; Bayón et al. 2015) and enables a stable, robust and accurate description of complex flow phenomena.
The VOF method is often used to describe mass transfer processes in CFD applications (Wang et al. 2018). Therefore, Haroun et al. (2010a, 2010b) have developed an approach to describe mass transfer processes across the interface between two fluids using the Henry coefficient for the VOF method. This approach has been implemented in OpenFOAM's solver interFoam by Nieves-Remacha et al. (2015), Yang et al. (2017) and Severin (2017), resulting in a solver which will be called interHarounFoam in the following.
A short outline of the driving biochemical processes leading to H2S formation shall be given to highlight important factors. When anaerobic conditions occur in the wastewater, sulphate-reducing bacteria, which reside in the biofilms of sewer walls, can reduce sulphate to sulphide (Sharma et al. 2008). From the biofilm, sulphide is then diffused into the wastewater as H2S. In the water, equilibrium conditions, depending on the pH value and temperature, determine which amounts of sulphide are present as H2S and as bisulphide ion (HS−); together they are described as total dissolved sulphide. The air–water equilibrium, which can be described by the Henry coefficient for a volatile compound such as H2S, can cause emissions of H2S from the water into the air phase. The rate of the transfer process is influenced by factors such as the flow velocities within the different phases, the pH value, temperature and the concentration of oxygen and nitrate. The Henry coefficient describes the relative amount of a volatile compound in the gas phase as a function of its relative occurrence in the water phase under equilibrium conditions and at constant temperature. The temperature dependency of Henry's law can be described by different equations, for example by the van 't Hoff equation. The concentration of H2S in the air phase defines the intensity of odour (Hvitved-Jacobsen et al. 2013).
As this overview of the relevant processes shows, a sole consideration of the Henry coefficient when describing H2S emissions is not sufficient. Therefore, relevant extensions have been made to the solver, resulting in a new, specialized solver, interH2SFoam. This solver is able to account for the temperature dependency of the Henry coefficient. Further extensions enable the user to describe the equilibrium between HS− and pH value in the water phase and to compute the partial pressure of H2Sg in the air phase in ppm in order to gain a better comparability between simulations and measured values. The assessment of turbulent flow effects on mass transfer will be subject to future research.
In the following, after an introduction on the methods used, the capabilities of the interH2SFoam solver are explored in three simple application examples of vertical one-dimensional flow. Then, mass transfer in a rectangular pipe is simulated. In a final example, the new solver is applied to a highly complex sewer geometry.
OpenFOAM version 2.4.0 has been used for the work presented in this paper. Additionally, a supplementary library called swak4Foam has been used to generate customized function objects to calculate the equilibrium conditions between H2S and HS− as well as the partial pressure in the air phase. This approach makes the use of this function optional for the user. Depending on the framework of the model, the user can then decide whether these functions are needed or not. The temperature dependency of the Henry coefficient of H2S has been directly implemented in the solver and makes a definition of the temperature as an input parameter mandatory.
The mass transfer solvers are based on the two-phase flow solver interFoam which is based on a VOF approach that considers both phases as one fluid with changing fluid properties. One set of Navier–Stokes equations is solved. The volume fraction of a phase is stored as an additional variable and the phases are distinguished by an additional transport equation. The equations are defined as follows (Rusche 2003).
The water surface is defined as the area where .
A turbulence model based on the Reynolds-averaged Navier–Stokes (RANS) equations (standard k-ɛ) is applied to consider the turbulent part, and the near-wall turbulence is modelled by so-called wall functions. More advanced turbulence models, such as large eddy simulations or direct numerical simulations, would offer the opportunity of resolving small-scale velocity variations but would come with the price of a highly increased computation time. As we expect that the application of more advanced turbulence effects would not change the insights on the equilibrium conditions of the mass transfer simulations addressed in this publication, a RANS turbulence model has been considered to be sufficient as well as the best way to save computational resources. Even with the standard k-ɛ turbulence model, the computation time of 10 seconds simulation for the complex sewer geometry amounted to 12 hours on 80 parallel processors using the high performance computing clusters of TU Berlin.
The accuracy of the hydrodynamic simulations has been assessed in Teuber et al. (in press).
In general, the transport of a passive tracer with a concentration c is examined with an advection–diffusion equation that can be implemented into the interFoam solver (see Equation (9)). The physical diffusivity as well as the turbulent Schmidt number , which defines the turbulent diffusivity coefficient , then have to be defined by the user (Equation (10)).
The diffusion coefficients for Daq and Dg are defined by the user. Note that these coefficients are temperature dependent, which has to be taken into account when defining the values.
The Henry coefficient, also known as Henry constant, is a temperature dependent variable which is reported in many different forms in the literature and is often expressed in different units.
Temperature dependency of Henry coefficient
The Henry coefficient depends on the overall temperature in the domain. Therefore, the temperature dependency has been added in a way that the solver takes one global temperature value as an input parameter.
Here, C is a temperature coefficient, which depends on the enthalpy of dissolution and is defined as 2,100 K (Sander 2015); is the standard temperature 298.15 K corresponding to 25 °C.
The equilibrium conditions are implemented following Hvitved-Jacobsen et al. (2013). The aim is to describe the water-phase concentration of H2S depending on the amount of total dissolved sulphide and the pH value since those values are usually measured in field investigations.
The equilibrium between H2S in the gas phase (H2Sg) and H2S in the water phase (H2Saq) is described by the Henry coefficient. In the water phase, an equilibrium between hydrogen sulphide H2Saq and bisulphide ion (HS−) exists, where the total amount of both is described as total dissolved sulphide.
Only H2Saq can be transferred across the air–water interface, not the ionized form HS−; however, usually the concentration of total dissolved sulphide and the pH value are measured. Therefore it is useful to derive a way to calculate the concentration of cH2Saq when cS and pH are given.
The value of pKa2 = 14.0 indicates that measurable amounts of the sulphide ion S2− only exist at a value above a pH of about 12. Therefore, only the equilibrium value of Ka1 is important for wastewater and has been implemented for the interHarounFoam solver in OpenFOAM using the utilities funkySetFields (for initial conditions) and funkySetBoundaryFields (as boundary conditions) from swak4Foam.
Note that the equilibrium constants Ka1 and Ka2 are temperature dependent (Yongsiri et al. 2004), which is not considered in the current version of the code. The value of Ka1 is a user-defined variable and the temperature dependency of the equilibrium constant has to be accounted for by defining the corresponding Ka1 value for the temperature analysed.
Calculation of partial pressure of H2Sg in ppm
The partial pressure of H2Sg is being computed using a function object in swak4Foam.
This conversion is only valid for the gas phase concentration; therefore the expression is multiplied by (1−α) in order to keep the conversion constrained to the air phase.
In the following, three different cases will be used to explore the possibilities of the existing solver, to validate the new features added to the solver and to show the importance of the model compared to the existing model approaches. For all simulations, at a temperature of 25 °C, physical diffusivities for H2S in water of 2.2 10−9 m²/s and in air of 1.74 10−5 m²/s are chosen.
The first setup is a quasi-one-dimensional cubic tank with the measurements 1 m × 1 m × 1 m bounded by upper, lower and sidewalls with no-slip conditions. The tank is partially filled with water (water depth = 0.5 m). Both fluids water and air are at rest. As an initial condition, an H2S concentration of cH2Saq = 1 mol/m³ in the water phase is given, the concentration in the air phase is cH2Sg = 0 mol/m³. The domain is discretized with 100 cells in y-direction, which is the vertical dimension of the domain, and 10 cells in x- and z-direction. At the bottom wall, a concentration source is assumed, using a fixed value boundary condition of 1 mol/m3. The top wall as well as the sidewalls is defined with zeroGradient conditions. This setup is used to illustrate the solver's capabilities in a simple setup. In a first example, mass transfer, as described with the existing interHarounFoam solver, is shown in a vertical one-dimensional case. Then, the extensions leading to the interH2SFoam solver are demonstrated in different examples using this first setup.
In a second setup, mass transfer in a rectangular duct is analysed using two well-documented examples of water–air pipe flow as described by Bentzen et al. 2016 (test cases no. 7 and 21). The investigated pipe has a length of 15.0 m, a height of 0.26 m and a width of 0.3 m with two different water depths and slopes. The air phase is only accelerated by the movement of the water surface. Bentzen et al. (2016) measured resulting velocity profiles in detail using laser Doppler anemometry velocity measurements. The flow characteristics of the two test cases analysed are listed in Table 1. The setup is a relatively simple three-dimensional setup of a pipe. It illustrates the applicability of the model to regular pipes. The computational domain consists of 307,970 cells. The inlet has been divided into two parts: one for the water phase and one for the air phase. For the water phase, a fixed discharge has been defined, and the phase fraction value α has been defined to be α = 1. The pressure boundary condition has been defined as null Neumann condition. For the air phase, a fixed pressure has been defined and the phase fraction value has been set to α = 0. The velocity has been defined using a null Neumann condition. At the outlet, a free outflow has been assumed. A fixed pressure has been defined and the remaining boundary conditions were defined as null Neumann conditions. At the walls, no-slip conditions were applied. Hydrodynamic simulations (without mass transfer) were run for 200 s, until quasi-steady-state conditions were reached; afterwards a concentration cH2Saq = 1 mol/m³ has been defined for the water at the inlet, assuming contaminated water flowing into the domain. The upper fluid then has a concentration of cH2Sg = 0 mol/m³ at the inlet; all remaining boundaries were defined with null Neumann conditions.
|Test no. .||Duct slope (%) .||Water depth (cm) .||Uaq (m/s) .||Ug (m/s) .||Reynolds number Reaq .||Reynolds number Reg .|
|Test no. .||Duct slope (%) .||Water depth (cm) .||Uaq (m/s) .||Ug (m/s) .||Reynolds number Reaq .||Reynolds number Reg .|
The third setup describes a complex sewer geometry with an overall length of 93.3 m, a width ranging from 6.0 m to 7.5 m and a sewer height between 4.3 m and 5.3 m. The setup is shown in Figure 6. This geometry has been simulated in OpenFOAM and compared to experimental results from a 1:20 scale model by Bayón et al. (2015) and Teuber et al. (in press). The setup describes a highly three-dimensional pipe diversion, including bends and geometry changes as well as a hydraulic jump. The computational mesh consists of 3,029,223 cells. The setup of boundary conditions is similar to the rectangular pipe of Bentzen et al. (2016). The hydrodynamic model has been simulated for 200 s, until steady-state conditions were reached. Then, a concentration cH2Saq = 1 mol/m³ has been defined for the water phase, and the simulations using the interH2SFoam solver have been carried out for a simulation time of 10 s. A temperature of 25 °C is assumed.
The results of the numerical simulations are presented in the following section.
RESULTS AND DISCUSSION
Saturation of H2S in a tank
Mass transfer modelling
In our first test case, we present the application of the model to a vertical one-dimensional problem. It illustrates the advantage of the new model in describing vertical concentration profiles in contrast to the existing horizontal one-dimensional approaches. The simplicity of the test case enables a first illustration of the model's capabilities. The simulation has been carried out assuming normal temperature (25 °C).
Figure 1 shows the presence of the two phases within the domain (α = 1: water, α = 0: air) and the development of the concentration profile over time. After t = 50 s, a decrease of the overall concentration in the water phase can be observed. This is due to the concentration jump at the interface, which has to be fulfilled by the solver. This concentration jump occurs in the first second due to a direct flux of concentration across the interface. After several seconds, the concentration in the water phase is re-established by the source term at the bottom, and after t = 1,000 s, a steady-state has developed and a constant concentration profile is achieved. The concentration in the water phase is then equal to the source term concentration and the air phase concentration is defined by the Henry coefficient. A detailed validation of the flux under transient conditions has been performed by Haroun et al. (2010a).
In this test case, we will analytically analyse the temperature dependency of the Henry coefficient, which has been implemented. The application example is based on example 4.2 in Hvitved-Jacobsen et al. (2013). The Henry coefficient at a temperature of 15 °C is being calculated.
Equilibrium conditions and unit conversion
In order to validate the solver extensions regarding the equilibrium conditions and the partial pressure in the air phase, example 4.3 by Hvitved-Jacobsen et al. (2013) is simulated. The resulting H2Saq and H2Sg concentrations for a measured concentration of dissolved sulphide cS- and pH value have been simulated. Again, the basic setup of the case is the same as for the first application example. A temperature of 15 °C, pKa1 = 7.0, pH = 7.0 and a dissolved sulphide concentration of 0.001 kg/m³ are given.
In Hvitved-Jacobsen et al. (2013), the Henry coefficient for the given temperature of 15 °C is assumed to be the same as the Henry coefficient from a previous calculation for a temperature of 20 °C, i.e. 433 atm. For comparing the analytical solution with the simulated values, the exact Henry coefficient for 15 °C has been calculated. Using this value and performing the same calculation steps with the corrected Henry coefficient, the analytical solution leads to a water phase H2S concentration of cH2Saq = 0.0075 mol/m³, a gas phase concentration of cH2Sg = 0.0027 mol/m³ and a corresponding partial pressure of pH2Sg = 66 ppm.
Figure 3 shows the results of the numerical simulations. In the water phase, the concentration of H2S is cH2Saq = 0.0075 mol/m³; in the air phase the concentration reaches a value of cH2Sg = 0.0027 mol/m3 and a corresponding partial pressure (in the gas phase) of pH2Sg = 66 ppm and thus agrees well with the analytical solution. The implemented approach predicts the resulting concentrations accurately.
Mass transfer in a rectangular channel
In this test case, we present mass transfer simulations in a rectangular pipe. Figures 4 and 5 present the resulting phase fraction, velocity and concentration profiles along the height of the domain in the middle of the pipe. The simulated velocity profiles indicate a good agreement with the measured values of Bentzen et al. (2016). For a pipe with a length of 15 m and the analysed flow velocities, the concentration profiles show that almost no mass transfer across the water surface into the air phase can be observed. This can be explained by small velocities in directions other than the main flow direction (i.e. in the yz-plane) which cause advective transport to occur mostly in the main flow direction (x-direction). Furthermore, the small diffusion coefficients cause mainly advective transport. This example raises the question as to how simulated mass transfer is influenced under highly turbulent conditions or in cases with higher velocities in the yz-plane, which will be analysed in the next example.
Complex sewer geometry
In a final example, the advantages of the new model are demonstrated by applying the solver to a complex and highly three-dimensional sewer geometry. The existing models are not public domain; therefore a direct comparison to the one-dimensional models is not possible, but the results of the CFD model will be used to highlight the advantages compared to the concepts of the existing approaches.
The results of the simulations at t = 10 s are displayed in Figure 6. Figure 6(a) gives an overview of the computational domain and the water phase behaviour. The location of highest turbulence occurs in the hydraulic jump, which is displayed in Figure 6(b) and 6(c). The velocity vectors in Figure 6(b) indicate the highly three-dimensional flow behaviour in this location and show the complex water surface movement. In Figure 6(c) and 6(d), the isosurfaces of the resulting H2S concentration in the domain are displayed. The concentration range between 0 mol/m³ and 1 mol/m³ has been divided into 10 surfaces. The value range inbetween is not displayed, leaving white spaces for better illustration of the surfaces. The contour plots show that a more diverse and highly three-dimensional concentration profile develops at the reach of the hydraulic jump. This indicates the increased mass transfer (i.e. higher distance of concentration isolines to the water surface) in the location of the hydraulic jump.
Because the existing model approaches are not public domain, a direct comparison to simulation results is not possible; however, the advantages of the new CFD-based mass transfer approach are the following:
In respect to the hydrodynamic behaviour, the new model can describe the three-dimensional flow velocities in the air and water phase. The sewer geometry analysed consists of a bent pipe structure with varying shapes and a hydraulic jump. A hydrodynamic one-dimensional approach would describe this geometry as one connection pipe between beginning and end point. The flow velocity would be calculated as a uniform value without accounting for the complexity of the geometry. The existing model would not account for the highly complex interaction of water and air phase in the hydraulic jump.
Regarding the mass transfer, the model would then account for advection and molecular diffusion and for turbulence in the free-stream flow areas as well as in drop structures in a very simplified way. This would lead to a simplified assumption of the actual mass transfer occurring in the pipe, since the effect of turbulence on the mass transfer is substantial. A validation of the actual mass transfer rate due to turbulence effects is performed in Teuber et al. (submitted).
Most sewer stretches in urban areas are not as complex as the previously shown example and wide networks without high levels of turbulence justify the use of one-dimensional models. However, locations of high turbulence can enhance H2S emissions, and the three-dimensional approach presented in this paper can help analyse the effect of local design aspects on the resulting H2S emissions and improve the sewer network design.
H2S emissions and their consequences are an important topic when considering urban drainage and the design of sewer networks. In the past, different model approaches, from empirical to conceptual, have been developed in order to describe and predict H2S emissions and resulting odour. These models are horizontal one-dimensional, therefore neglecting the occurrence of three-dimensional effects.
In this publication, a model approach has been introduced that can describe H2S emissions across the water surface using a mass transfer approach based on the Henry coefficient, which is implemented in the open source software OpenFOAM. Two-phase flow has been simulated using a VOF method. The solver has been extended by different key features that are crucial when describing H2S emissions. The temperature dependency of the Henry coefficient has been taken into account. Equilibrium conditions between HS− and H2S are described and enable the usage of the measured value for total dissolved sulphide and the pH value as input parameters. The solver also computes the partial pressure of H2S in the gas phase based on the simulated concentration of H2Sg.
The new solver has been applied to different simple test cases and the results have been compared to analytical solutions. Furthermore, it has been applied to a highly complex three-dimensional test case to highlight the advantages of the new model approach. Compared to one-dimensional formulations, it can account for highly complex flow effects in a sewer stretch and describe mass transfer in such environments. The analysis of the results showed an increased mass transfer in the location of highest turbulence, which agrees with existing observations. The exact quantification of local mass transfer rates has been validated in Teuber et al. (submitted) and has led to a good agreement with experimental results.
Overall, the new solver enables an analysis of mass transfer in complex three-dimensional test cases, the description of which has so far only been possible with major simplifications.
Future research will deal with further extensions of the solver to account for temperature effects in the fluids and reactive transport modelling.
The complex sewer geometry has been computed on the supercomputers of Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen in Berlin. The funding provided by the German Research Foundation (DFG) within the Research Training Group ‘Urban Water Interfaces’ (GRK 2032-1) is gratefully acknowledged. We thank Professor Arnau Bayón for providing the experimental data and the mesh for the complex sewer geometry.