https://orcid.org/0000-0002-8959-8808
ABSTRACT
This study explores the computational fluid dynamics (CFD) simulation of oxygen (O2) and hydrogen sulfide (H2S) mass transfer in a highly turbulent stirring tank. Using the open-source software OpenFOAM, we extended three-dimensional two-phase flow solvers with a rotating mesh feature to model the mass transfer processes between the water and air phases. The accuracy of these simulations was validated against experimental data, demonstrating a strong agreement in the mass transfer rates of H2S and O2. The investigation highlights the impact of turbulence on mass transfer coefficients, confirming the reliability of the solvers for predicting mass transfer in turbulent conditions. The results suggest that these CFD models can serve as effective tools for understanding and optimizing sewer system designs. Additionally, the study highlights the potential of numerical simulations to reduce the need for extensive and potentially hazardous laboratory experiments.
HIGHLIGHTS
Extended OpenFOAM solvers effectively simulate turbulent H2S and O2 and the first integration of rotating mesh in computational fluid dynamics (CFD) for mass transfer.
A high correlation between CFD predictions and experimental data confirms model accuracy.
Turbulence increases mass transfer rates, crucial for accurate sewer system modeling.
CFD results show consistent volumetric mass transfer rates, aligning well with established experimental data.
Future work should focus on the influence of gas phase resistance.
INTRODUCTION
Modern cities face increasing challenges due to the growing centralization of sewer systems. One of these problems is the occurrence of odor and corrosion due to hydrogen sulfide (H2S) emissions across the wastewater–air interface, which cause a health risk for sewer workers as well as high costs for sewer maintenance. The maintenance costs for sewer networks, including renovation and renewal, were estimated to be around 2 billion Euros in Germany for the period 2019–2023 (Berger et al. 2020). Depending on various factors, including the pH level, temperature, oxygen (O2), and sulfate content at the (waste)water– (sewer)air interface, H2S can be released from the water phase into the air phase. The transfer of O2 in between this interface, called reaeration, should also be considered while modeling. The level of aerobic and anaerobic processes in wastewater depends on the degree of reaeration during transport (Yongsiri et al. 2005; Hvitved-Jacobsen et al. 2013). Due to O2 transfer across the air–water interface, aerobic conditions may predominate in the wastewater phase. In pressure pipes, O2 can be depleted, resulting in anaerobic conditions and sulfate reduction. When wastewater is pumped from the wet well, an intense gas–liquid transfer may occur downstream, at the discharge point into a gravity sewer. When transferred to sewer air, H2S causes sewer corrosion by oxidation to sulfuric acid on the inner surface of concrete pipes as well as odor nuisance and health hazards when emitted to the surrounding environments (Zhang et al. 2023). Therefore, studying the mass transfer and transport of both H2S and O2 becomes important for understanding the odor and corrosion-causing processes in sewers.
In the past years, significant progress has been made in the field of understanding H2S emissions in sewers. On the one hand, conceptual model approaches have been developed to describe the occurrence of odor and corrosion (Hvitved-Jacobsen et al. 2013). On the other hand, numerous experiments have been conducted in an attempt to understand contributory factors such as the influence of turbulence on H2S mass transfer across the water surface (Carrera et al. 2017; Matias et al. 2017; Sun et al. 2023, 2024).
Being able to model and directly quantify these emissions will help to improve existing models. Moreover, investigating a wider range of cases and factors will be of great assistance in the design of sewers, especially drop structures or aeration systems.
State-of-the-art modeling approaches to describe these in-sewer processes have majorly been one-dimensional (Hvitved-Jacobsen et al. 2013). However, previous research has shown that the flow processes occurring are three-dimensional, for example, air phase velocities in the headspace of a circular sewer (Edwini-Bonsu & Steffler 2006), making three-dimensional computational fluid dynamics (CFD) models essential for predicting small-scale effects despite the significantly higher computational effort. One benefit of these models is the detailed analysis of turbulent hotspots in sewers they enable.
High levels of turbulence in sewers can be found at drop structures. Transport and mass transfer phenomena are highly influenced by the level of turbulence; therefore, correct quantification of these effects on the H2S and O2 mass transfer is crucial for the development of reliable models and can only be done using a three-dimensional approach. This is why Teuber (2020) and Dixit (2024) have extended a CFD solver by Haroun et al. (2010a) within the open-source software tool OpenFOAM to account for H2S and O2 mass transfer in sewer systems. The solver has been modified to account for turbulent effects when using Reynolds-averaged Navier-Stokes (RANS) equations such as k–ε-turbulence models.
Turbulent mass transfer has been investigated and quantified on the laboratory scale in different publications (e.g., Matias et al. 2017). Laboratory-scale investigations on mass transfer rates have been performed multiple times using stirring tanks (Wu 1995; Carrera et al. 2017). In the following, the case setup of a stirring tank for which experimental results have been generated by Pacheco Fernández et al. (2020) and Tang (2019) is used to evaluate our solvers' performance. With this setup, the current applicability to real sewer systems is limited. However, the setup enables us to minimize the number of influencing factors and to investigate the model's capability to describe mass transfer under different turbulent conditions.
Geometry and mesh used for mass transfer simulations.
This paper starts with an overview of the methods and materials. Then, the results are evaluated and discussed. In a quantitative analysis, the simulated results are then compared to experimental results by Pacheco Fernández et al. (2020). Finally, the findings are compared to investigations by Carrera et al. (2017) and Springer et al. (2020).
METHODS AND MATERIALS
Numerical model
OpenFOAM V6 was used for the simulations presented in this publication. The official version was extended with two custom solvers for two-phase simulations and mass transfer, which also account for the temperature dependency of the Henry coefficient. The first solver describes H2S (interH2SFoam) (Teuber et al. 2019b); the second, O2 mass transfer (interO2Foam) (Dixit 2024).
Hydrodynamic simulations
Hydrodynamic simulations are based on the two-phase flow solver interFoam, which uses a volume of fluid (VOF) approach (Hirt & Nichols 1981) that considers both phases as one fluid with changing fluid properties. The implementation of the VOF approach in OpenFOAM has been extensively documented in several publications, e.g. Rusche (2003).
is a volume fraction or indicator function [–] and the water surface is defined as the area where α = 0.5. Variations of the free water surface (Lopes et al. 2018) were considered when analyzing the mass transfer rates but were of negligible effect.The accuracy of such hydrodynamic simulations for applications in sewer networks has been assessed by Teuber et al. (2019a).
The k–ε turbulence model was chosen for all simulations, as it has been found most accurate for similar simulations in Teuber (2020).
Mass transport and transfer simulations
Transport is considered by defining the H2S and O2 concentration as a passive tracer with an advection-diffusion equation (Equation (1)).
is the velocity field (m/s); t is time (s); Dphys is the physical diffusivity (m²/s); Dturb is the turbulent diffusivity (m²/s); He is the Henry coefficient (–).
The diffusion coefficients for Dphys,L and Dphys,G are defined by the user. Note that the temperature dependency of these variables has to be accounted for.
The Henry coefficient depends on the temperature in the whole domain. Therefore, the temperature dependency has been added in a way that the solver takes one global temperature value as an input parameter.
To combine the rotating mesh with mass transfer simulations, the respective functionalities for mass transfer had to be transferred to the interDyMFoam solver for the initial version of the solver, resulting in a new specialized solver interDyMH2SFoam in OpenFOAM version 2.4.0 (Teuber 2020). Most dynamic solvers, denoted by the prefix ‘Dy(namic)’, have been since merged or integrated into their corresponding non-dynamic counterparts in OpenFOAM V6. Therefore, the updated solver utilized for simulations was the extension of the in-built solver interFoam (interH2SFoam for H2S and interO2Foam for O2) (Dixit 2024).
is the initial concentration at t=t0. In this publication, a has been extracted from CFD simulations and 
was fitted by automated curve fitting using the programming language Python. For O2, 
has been calculated accordingly using the following equation:
Geometry
In our model, the geometry consists of a circular stirring tank in which a Rushton turbine is placed with a defined stirring rate (Figure 1). For the construction of the geometry for the simulation setup, all the dimensions were manually measured from the lab-scale setup. The OpenFOAM solver interDyMFoam, utilized in a modified version here, is capable of handling mesh motion using both static and rotating mesh domains. Therefore, the mesh motion has to be defined separately by specifying the rotation speed.
Experimental investigations
Experiments were designed as in the published works of Carrera et al. (2017) and were carried out by Pacheco Fernández et al. (2020). The experiment was divided into two methodologies, namely, an absorption method for measuring the mass transfer coefficient of O2 and a desorption method for measuring the mass transfer coefficient of H2S. The principle of the absorption method was to obtain the relationship between the time and the concentration of dissolved O2 in the water phase after an initial concentration was set to zero. The desorption method determines the mass transfer by the decrease of the previously dissolved H2S concentration in the water phase with time.
Experimental setup for the H2S mass transfer from the liquid to the gas phase (top) and O2 reaeration (bottom) experiment.
Experimental setup for the H2S mass transfer from the liquid to the gas phase (top) and O2 reaeration (bottom) experiment.
Case study
The tank had a maximum volume of 52.8 L. Different probes were installed at a depth of 8 cm below the water surface. The stirrer had a range from 50 rpm (revolutions per minute) to 2,200 rpm. 300, 400, and 500 rpm were chosen for the experiments of O2 and H2S mass transfer. The stirring speed results in Reynolds numbers between 50,000 and 83,333, indicating turbulent flows. Table 1 gives the details of the dimensions of the experimental setup, and Table 2 lists the relevant fluid properties.
Experimental setup data for both experiments
| Shape | Height | Maximum volume | Dout | Din | Dstirrer | Water height | |
|---|---|---|---|---|---|---|---|
| H2S exp. | O2 exp. | ||||||
| Cylinder | 80 cm | 52.8 L | 30 cm | 29 cm | 10 cm | 16 cm | 15 cm |
| Shape | Height | Maximum volume | Dout | Din | Dstirrer | Water height | |
|---|---|---|---|---|---|---|---|
| H2S exp. | O2 exp. | ||||||
| Cylinder | 80 cm | 52.8 L | 30 cm | 29 cm | 10 cm | 16 cm | 15 cm |
Fluid properties as defined for the CFD simulations
| Water phase | Air phase | |
|---|---|---|
| Density ρ (kg/m3) | 1,000 | 1.2 |
| Kinematic viscosity ν (m2/s) | 1.0 × 10−6 | 1.48 × 10−5 |
| Diffusion coefficient Dphys (H2S) (m2/s) | 1.8 × 10−9 | 1.74 × 10−5 |
| Diffusion coefficient Dphys (O2) (m2/s) | 2.42 × 10−9 | 1.98 × 10−5 |
| Water phase | Air phase | |
|---|---|---|
| Density ρ (kg/m3) | 1,000 | 1.2 |
| Kinematic viscosity ν (m2/s) | 1.0 × 10−6 | 1.48 × 10−5 |
| Diffusion coefficient Dphys (H2S) (m2/s) | 1.8 × 10−9 | 1.74 × 10−5 |
| Diffusion coefficient Dphys (O2) (m2/s) | 2.42 × 10−9 | 1.98 × 10−5 |
RESULTS
Domain discretization and convergence
To describe the computational domain, different meshes have been generated using a variety of mesh generation tools. The details for the meshes generated and the computational effort required for 1 s of simulation time are presented in Table 3. The table lists the number of cells and the respective time taken for 1 s of simulation. All the simulations were conducted using the same number of compute nodes and cores of the North-German Supercomputing Alliance (HLRN). One compute node with 96 cores was utilized in all the simulations as it was found that increasing the computation power further increased the computation time. For a structured grid to be constructed, the minimum length of the cell was equal to the width of the blade (1 mm or 0.001 m), which drastically increased the number of cells. Varying the cell size distribution with smaller cells around the blade and increasing the cell size for the surrounding reactor caused a significant gradient in cell size, which either took more time to simulate or resulted in an unstable simulation. Therefore, a uniform unstructured grid was considered for the mass transfer simulations. The initial and final residuals for the different parameters showed a good level of convergence, resulting in the conclusion that the solution can be seen as tightly converged. The final mesh is highlighted in Table 3 and was discretized using the open-source software Salome-Meca (acronym in French, Simulation numérique par Architecture Logicielle en Open source et à Méthodologie d'Évolution, translation: Numerical Simulation by Computing Architecture in Open Source and with Evolving Methodology).
Details for different meshes generated and their respective simulation times, final mesh in bold
| Tool | Mesh | Number of cells | Approximately computation time (s) required per 1 s simulation time | |
|---|---|---|---|---|
| Salome-Meca | Unstructured | Fine | 179,213 | 7,200 s |
| Unstructured | Moderate | 63,237 | 300 s | |
| GMESH | Structured | Fine | 243,612 | 3,600 s |
| snappyHexMesh | Structured | Coarse | 421,725 | 2,400 s |
| Structured | Moderate | 3,370,568 | 10,800 s | |
| Tool | Mesh | Number of cells | Approximately computation time (s) required per 1 s simulation time | |
|---|---|---|---|---|
| Salome-Meca | Unstructured | Fine | 179,213 | 7,200 s |
| Unstructured | Moderate | 63,237 | 300 s | |
| GMESH | Structured | Fine | 243,612 | 3,600 s |
| snappyHexMesh | Structured | Coarse | 421,725 | 2,400 s |
| Structured | Moderate | 3,370,568 | 10,800 s | |
H2S emissions
Comparison of measured data to results obtained from simulations for H2S for different stirring rates ((a) 300 rpm, (b) 400 rpm, and (c) 500 rpm).
Comparison of measured data to results obtained from simulations for H2S for different stirring rates ((a) 300 rpm, (b) 400 rpm, and (c) 500 rpm).
O2 reaeration
Comparison of measured data to results obtained from simulations for O2 for different stirring rates ((a) 300 rpm and (b) 500 rpm)).
Comparison of measured data to results obtained from simulations for O2 for different stirring rates ((a) 300 rpm and (b) 500 rpm)).
Model efficiency of simulations
Table 4 lists the calculated model efficiency criteria for the simulated results to the averaged measured data (O2 values at 400 rpm are not available due to the absence of corresponding experimental data). The resultant Nash–Sutcliffe model efficiency (NSE) for H2S (approximately 0.8) and O2 (approximately 0.9) shows that the solvers can produce an acceptable representation of mass transfer processes under turbulent conditions. The calculated root mean square error (RMSE) and normalized RMSE (normalized by the difference between minimum and maximum concentration) also indicate reasonable acceptance between the modeled and the experimental results. The correlation coefficients (Pearson's R) also represent a good fit between the two results (r between 0.876 and 0.995). R2 values between 0.77 and 0.99 show that the modeled data match the actual data to a satisfactory degree.
Different model efficiency criteria for the H2S and O2 solver with respect to the average measurements
| Stirring rate | NSE | RMSE | r | R2 | ||||||
|---|---|---|---|---|---|---|---|---|---|---|
| (RPM) | H2S (−) | O2 (−) | H2S (mg/L) | H2S normalized (−) | O2 (mg/L) | O2 normalized (−) | H2S (−) | O2 (−) | H2S (−) | O2 (−) |
| 300 | 0.66 | 0.91 | 0.002 | 0.274 | 0.004 | 0.128 | 0.932 | 0.995 | 0.87 | 0.99 |
| 400 | 0.85 | – | 0.001 | 0.118 | – | – | 0.949 | – | 0.90 | – |
| 500 | 0.83 | 0.93 | 0.002 | 0.161 | 0.012 | 0.102 | 0.876 | 0.994 | 0.77 | 0.98 |
| Stirring rate | NSE | RMSE | r | R2 | ||||||
|---|---|---|---|---|---|---|---|---|---|---|
| (RPM) | H2S (−) | O2 (−) | H2S (mg/L) | H2S normalized (−) | O2 (mg/L) | O2 normalized (−) | H2S (−) | O2 (−) | H2S (−) | O2 (−) |
| 300 | 0.66 | 0.91 | 0.002 | 0.274 | 0.004 | 0.128 | 0.932 | 0.995 | 0.87 | 0.99 |
| 400 | 0.85 | – | 0.001 | 0.118 | – | – | 0.949 | – | 0.90 | – |
| 500 | 0.83 | 0.93 | 0.002 | 0.161 | 0.012 | 0.102 | 0.876 | 0.994 | 0.77 | 0.98 |
In general, all the model efficiency criteria analyzed show that the solvers can be a useful tool for quantitatively predicting the mass transfer of H2S and O2 in highly turbulent cases.
MASS TRANSFER COEFFICIENTS
is illustrated in Figure 5, with the x-axis presented on a logarithmic scale for clarity. Comparison to data obtained from Springer et al. (2020) suggests a correlation between the Reynolds number and 
. This agrees with the findings of Springer et al. (2020), which suggested a correlation between the interfacial turbulent kinetic energy (TKE) and 
. This relationship has been investigated within the scope of this publication, but we have found that the relation between TKE and mass transfer coefficient KL appears to depend on additional factors such as the ratio between stirrer inundation and tank diameter. Further investigations would be needed to quantify the exact relationship. An exact comparison to the findings by Springer et al. (2020) is difficult due to the different orders of magnitude in the Reynolds numbers investigated. 
Influence of Reynolds number on mass transfer coefficient 
as published in Springer et al. (2020) and results of our CFD simulations.
Influence of Reynolds number on mass transfer coefficient as published in Springer et al. (2020) and results of our CFD simulations.
/
which equals to 0.70. This coefficient can be expressed as a function of the diffusion coefficients of the two media and an exponent n as described by Hvitved-Jacobsen et al. (2013):
/
values range between 0.64 (for slow-flowing conditions) and 0.80 (for turbulent conditions) (US EPA 1974). Hvitved-Jacobsen et al. (2013) reported a higher 
/
of 0.87, leading to 
/
between 0.87 and 0.93. Overall, the calculated values align with the observations of Carrera et al. (2017), who determined 
for mass transfer in a stirring tank. The ratio between mass transfer coefficients from Carrera et al. (2017) and our results are presented in Figure 6. The plot demonstrates a similar trend in the effect of stirring rate on mass transfer, consistent with Carrera et al.'s observations. Furthermore, the individual values for the ratio KL,2S/KL,O2 for each stirring rate, when plotted with the observation of Carrera et al. (2017), lie in the same range (represented as a band in Figure 6). 
Influence of stirring rate on 
as published by Carrera et al. (2017) and results of our CFD simulations.
Influence of stirring rate on as published by Carrera et al. (2017) and results of our CFD simulations.
DISCUSSION
The CFD simulations carried out suggest an overestimation of mass transfer for O2 and H2S. These can be explained by different aspects. First, due to the expensive computational cost for O2 and H2S mass transfer simulations (1 h on a supercomputer with 98 nodes for 1 s of simulation time) and the complex setup of the H2S experiments, only a limited amount of data is available. The comparison between simulated and measured O2 and H2S concentrations therefore focuses on the first 300 s without reaching saturation. Table 4 suggests a reasonable level of accuracy between simulations and experiments. Longer simulations in simpler case setups might enable further insight and better comparability of the KL-values. To reach this goal, a simpler impeller setup leading to different flow patterns could be used. Instead of the axial flow impeller used here, a radial impeller could lead to flow phenomena that could be captured by a simpler computational domain (Shewale & Pandit 2023).
Furthermore, Matias et al. (2018) have found that, under turbulent conditions, the gas phase resistance plays an important role for substances with a moderate Henry coefficient such as H2S. Ignoring this gas phase resistance can lead to an overestimation of mass transfer (up to 13%). For substances with a lower gas phase resistance such as O2, this effect is less relevant but might still lead to an overestimation of mass transfer.
CONCLUSIONS
The influence of turbulence on H2S mass transfer has been subject to intensive research in the past years. So far, a general understanding has been gained from experimental investigations, but the influence has not been simulated using a numerical three-dimensional model with the RANS turbulence model. This is the first publication to address the solver's ability to accurately simulate H2S or O2 mass transfer across the air–water interface under turbulent conditions. The results agree well with experimental measurements and findings from existing literature.
This investigation is relevant to mass transfer processes in urban water infrastructure. Being able to describe H2S emissions in a stirring tank with a CFD model cannot only enable us to investigate complex phenomena and develop design improvements in sewer systems but could also help avoid lab experiments that come with a health risk and decrease the number of experiments necessary (Carrera et al. 2017). The utilization of the two solvers (interH2SFoam and interO2Foam) in this study enables the quantification of mass transfer for both species within a complex turbulent system. Moreover, a comprehensive assessment of the influence of turbulence on the mass transfer phenomenon was successfully conducted. In comparison to the findings presented in previous publications, as well as laboratory experiments, a noteworthy level of agreement is observed. These results enhance the reliability of the current study and reinforce the potential of these solvers as suitable tools for investigating mass transfer phenomena in turbulent systems.
To enhance the robustness of the model and facilitate comprehensive comparisons with experimental findings, further use cases with sewer-specific geometries would be desirable. That is why future research will focus on mass transfer processes in drop structures.
ACKNOWLEDGEMENTS
The modeling work was performed using high-performance computers at the HLRN. The computational resources provided are gratefully acknowledged.
FUNDING
This research was funded by the German Research Foundation (DFG) within the Research Training Group ‘Urban Water Interfaces’ (GRK 2032-1/2). The funding is gratefully acknowledged.
DATA AVAILABILITY STATEMENT
All relevant data are available from an online repository or repositories, see Teuber (2020) and Dixit (2024).
CONFLICT OF INTEREST
The authors declare there is no conflict.
%20cropped.png?versionId=7151)