Modular CFD Solver for Hydrogen Sulfide and Oxygen Transport: Implementation in OpenFOAM 11

13th Urban Drainage Modelling Conference, Innsbruck (Austria), September 2025 Modular CFD Solver for Hydrogen Sulfide and Oxygen Transport: Implementation in OpenFOAM 11 Katharina Teuber1*, Abhinav Dixit 2 & Reinhard Hinkelmann2 1Department of Civil Engineering, Geoinformation and Health Technology, Jade University of Applied Sciences, Ofener Str. 16/19, 26121 Oldenburg, Germany, ORCID: https://orcid.org/0000-0002-8959-8808 2Chair of Water Resources Management and Modeling of Hydrosystems, Technische Universität Berlin, GustavMeyer-Allee 25, 13355 Berlin, Germany *Corresponding author email: [email protected] This work is licensed under CC BY 4.0. To view a copy of this license, visit https://creativecommons.org/licenses/by/4.0/ Abstract Numerical modeling environments are subject to continuous evolution, leading to frequent software updates that often render custom solver modifications obsolete. This work presents the solver massTransferVoF, which is fully compatible with OpenFOAM 11. The solver describes mass transfer processes of hydrogen sulfide and oxygen across the water surface in two-phase flow simulations using a volume of fluid approach. The implementation process, challenges in maintaining backward compatibility, and the advantages of modular solvers over legacy structures are discussed. Validation results of three different test cases indicate that the new solver maintains consistency with previous implementations while offering access to new functionalities, improved maintainability and futureproofing for upcoming OpenFOAM versions. Highlights • Development of modular solver massTransferVoF compatible with OpenFOAM 11. • Capable of simulating H2S and O2 mass transfer. • Validation cases: 2D diffusion (single phase), 1D diffusion across interface (equilibrium and dynamic). Introduction Sewer networks are a key part of the urban underground infrastructure and enable a safe transport of wastewater towards wastewater treatment plants. Corrosion is the primary factor impacting the longevity of concrete sewers, globally leading to financial losses in the range of billions of dollars annually. In Flanders, Belgium, these costs even amount to roughly 10% of the total cost associated with sewer treatment, indicating large cost saving potentials for appropriate mitigation measures (Wani et al., 2025). Concrete corrosion processes are mainly driven by mass transfer occurring across the wastewatersewer atmosphere interface, namely hydrogen sulfide (H2S) emission and reaeration through oxygen (O2) mass transfer (Sharma et al, 2008, Hvitved-Jacobsen et al., 2013). Depending on the concentration of H2S in the sewer atmosphere, odorous smells or even health risks for sewer workers might occur. There is a wide range of models predicting H2S emissions and subsequent corrosion processes as well as ventilation effects. The well-known WATS model can be considered as state of the art model for prediction of concrete corrosion. This model considers 1D hydraulic processes (Hvitved-Jacobsen et al., 2013). Current research in the field focusses on improving simulation time (Shi et al., 2025) and the application of data-driven approaches (Yin et al., 2024, Wani et al., 2025). Other studies use 3D 13th Urban Drainage Modelling Conference, Innsbruck (Austria), September 2025 Computational Fluid Dynamics (CFD) to describe more complex hydraulic phenomena and their impact on emission and reaeration processes (Springer et al., 2020, Teuber et al., 2025). Recent publications show that the turbulent effects are a main driver of mass transfer processes in sewers. CFD models offer the possibility to investigate these turbulent processes and relevant turbulent properties (such as turbulent kinetic energy, turbulent dissipation etc.) in detail. Springer et al. (2020) derived a correlation between the mass transfer coefficient and the interfacial turbulent kinetic energy, while Teuber et al. (2025) propose a correlation between mass transfer coefficient and Reynolds number. The studies mentioned focus on a limited number of experimental and CFD investigations. A large number of structures would be beneficial to derive general conclusions. In order to answer these questions, suitable CFD solvers are necessary. When dealing with CFD simulations, OpenFOAM is an open source simulation suite, allowing for custom solver development. Teuber (2020) and Dixit (2024) published two OpenFOAM solvers (interH2SFoam and interO2Foam) for OpenFOAM versions 2.4.0 and 6, which are based on the volume of fluid (VoF) solver interFoam and target the simulation of H2S and O2 mass transfer processes across the water-air interface in sewers. The solvers are based on the framework developed by Haroun et al. (2010) and contain different customizations for applications in sewer systems. In contrast to closed source frameworks, these solvers within the OpenFOAM environment enable a simulation of the actual mass transfer occurring within the analysed systems. Since their publication, the OpenFOAM simulation suite has undergone major developments. One of the newest versions is OpenFOAM 11, released in July 2023. This version introduced a significant restructuring of solver architecture by transitioning to a modular framework. While this enhances solver flexibility and maintainability, it necessitates the adaptation of existing user-developed solvers to comply with the new structure. The transition presents both opportunities and challenges: on one hand, it simplifies future updates and keeps up with the newest developments, i.e. new dynamic meshing facilities, but on the other hand, it requires substantial modifications to existing custom solvers. In this study, we therefore present the adaptation of interH2SFoam and interO2Foam from their most recent implementation in OpenFOAM 6 (Dixit, 2024) to OpenFOAM 11 under the name massTransferVoF. Three validation cases show, that the solver is able to produce reasonable results. Methodology The previously implemented solvers interH2SFoam and interO2Foam are based on the two-phase solver interFoam. Within OpenFOAM 11, the former interFoam solver has been renamed to incompressibleVoF. To improve the solver’s usability, the previous solvers interH2SFoam and interO2Foam have been combined into one. The user may specify the species transfer model to use (H2S or O2) or define a user-specific Henry coefficient. Following the new nomenclature, the solver was renamed to massTransferVoF. The solution and implementation procedure are outlined in the following. Solution procedure Hydrodynamic simulations are based on a volume of fluid 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). The governing equations are solved depending on a phase fraction value α [-] which describes the phases present in each cell (α = 0.0 – air, α = 0.5 – water surface, α = 1.0 – water). Scalar transport of different species – such as H2S and O2 – is described using an advection-diffusion equation. Based on this, in solvers (ii) and (iii), mass transfer at the water surface is described by 13th Urban Drainage Modelling Conference, Innsbruck (Austria), September 2025 𝜕𝜕𝜕𝜕 𝜕𝜕𝜕𝜕 + 𝛻𝛻 ∙ �𝑼𝑼 � �  𝜕𝜕� =∇ ∙(�𝐷𝐷𝑝𝑝ℎ𝑦𝑦𝑦𝑦 + 𝐷𝐷𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡�∇𝜕𝜕 +𝜙𝜙) (1) 𝜙𝜙=−(𝐷𝐷𝑝𝑝ℎ𝑦𝑦𝑦𝑦 + 𝐷𝐷𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡)( 𝜕𝜕(1 −𝐻𝐻𝐻𝐻) 𝛼𝛼+𝐻𝐻𝐻𝐻(1 − 𝛼𝛼))∇𝛼𝛼 (2) Where C is the tracer concentration [mol/m³]; 𝑼𝑼 � �  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 [-]. Like density and viscosity in the VoF approach, the concentrations and diffusion coefficients are considered as single-phase properties depending on the phase fraction value α: 𝜕𝜕=𝛼𝛼𝜕𝜕𝐿𝐿+ 𝜕𝜕𝐺𝐺 (1 − 𝛼𝛼) (3) where the subscripts L and G denote the fluids water (L - liquid) and air (G - gas). The physical diffusivity Dphys is calculated by using a harmonic average: 𝐷𝐷𝑝𝑝ℎ𝑦𝑦𝑦𝑦 = ( 𝐷𝐷 𝑝𝑝ℎ𝑦𝑦𝑦𝑦,𝐿𝐿 𝐷𝐷 𝑝𝑝ℎ𝑦𝑦𝑦𝑦,𝐺𝐺 𝛼𝛼𝐷𝐷𝑝𝑝ℎ𝑦𝑦𝑦𝑦,𝐿𝐿+ (1 − 𝛼𝛼)𝐷𝐷𝑝𝑝ℎ𝑦𝑦𝑦𝑦,𝐺𝐺 ) (4) The diffusion coefficients for Dphys,L and Dphys,G are defined by the user. The temperature dependent Henry coefficient has been implemented using the van’t Hoff equation following Sander (2023). Implementation procedure Two primary approaches were considered for implementing massTransferVoF in OpenFOAM 11: - Backward compatibility and - Modular solver implementation. OpenFOAM 11 claims backward compatibility, but in practice, numerous dependencies, header file changes, and deprecated functions lead to significant compilation errors. Fixing these issues manually proved time-consuming and unreliable. Given the limitations of backward compatibility, the solvers were restructured using OpenFOAM's new modular solver framework. This transition allows for clearer separation of concerns within the code, enhanced debugging capabilities, and improved adaptability to future OpenFOAM updates. The chosen methodology prioritized rewriting the solver to align with modular solver best practices. This involved restructuring key components such as mesh setup, time handling, field definitions, and solver loops into distinct modular files. Case studies Three simple test cases have been chosen to validate the accuracy of the solver implementation. The results of each test case can be validated by comparison with reference cases. Each case setup will be introduced in the corresponding subsection in the results section. Results and discussion Attempts to retain backward compatibility led to numerous compilation errors, making it clear that a modular approach was the more effective solution. Our main results are presented in the following subsections. As can be seen, there is a good agreement between CFD results described in the previous section and reference solutions. 13th Urban Drainage Modelling Conference, Innsbruck (Austria), September 2025 Case 1: no flow, 2D diffusion, one phase 2D tracer spreading of a circular tracer injected on a rectangular plate has been simulated using the advection-diffusion equation implemented in the solver. Since there is no flow in this test case, the advection part of the equation is not considered. The rectangular plate has a length and width of 5 m, a numerical grid containing rectangular cells with a cell size of 0.01m x 0.01m has been chosen. As reference solution, the 2D diffusion equation has been solved using the finite-difference method using a FTCS scheme (forward differences in time, central differences in space) based on Hill (2020). The method has been implemented in Python. A diffusion coefficient of Dphys = 5 · 10-6 m²/s has been chosen. The simulation time has been set to t = 20,000 s, the results were compared at t = 0 s, t = 5,000 s, t = 10,000 s and t = 15,000 s. Figure 1 shows the 2D tracer spreading by diffusion over time. Figure 2 displays a comparison between FTCS scheme and the OpenFOAM solver massTransferVoF. The results agree well. It can be derived that diffusion processes are accurately described using the current solver implementation. Figure 1. Case 1 – 2D diffusion leads to tracer spreading over time. 13th Urban Drainage Modelling Conference, Innsbruck (Austria), September 2025 Figure 2. Case 1 - FDM and OpenFOAM solution produce similar results for each time step analysed illustrating the accuracy of the algorithm. Case 2: no flow, 1D diffusion across horizontal interface (equilibrium conditions) 1D diffusion in vertical direction is considered following Haroun et al. (2010). The diffusion occurs across a horizontal interface between a liquid and a gas phase which are both at rest. The interface is located at y = e = 0.5 m, the overall height of the domain is y = 1.0 m. The concentrations are imposed at y = 0.0 m in the liquid phase CL(y=0.0m) = CL0 and at y = 1.0 m in the gas phase CG(y=1.0m) = CG0. Concentration jumps are considered for different Henry coefficients (He = 0.1, He = 10, H2S and O2). The chosen diffusivities follow the species characteristics for H2S and O2 and follow a ratio Dphys,L/Dphys,G = 10 for the other simulations (Table 1). At equilibrium, the exact solution for the concentration is 𝜕𝜕𝐿𝐿= 𝜕𝜕 𝐺𝐺 0−𝐻𝐻𝐻𝐻𝜕𝜕 𝐿𝐿 0 𝐻𝐻𝐻𝐻 + 𝐷𝐷𝑝𝑝ℎ𝑦𝑦𝑦𝑦,𝐿𝐿 𝐷𝐷 𝑝𝑝ℎ𝑦𝑦𝑦𝑦,𝐺𝐺 𝑦𝑦 𝐻𝐻+𝜕𝜕𝐿𝐿 0 for 0.0m ≤ y ≤ e (5) 𝜕𝜕𝐺𝐺= 𝜕𝜕 𝐺𝐺 0−𝐻𝐻𝐻𝐻𝜕𝜕 𝐿𝐿 0 𝐻𝐻𝐻𝐻 𝐷𝐷𝑝𝑝ℎ𝑦𝑦𝑦𝑦,𝐿𝐿 𝐷𝐷 𝑝𝑝ℎ𝑦𝑦𝑦𝑦,𝐺𝐺 + 1 𝑦𝑦 − 2𝐻𝐻 𝐻𝐻+𝜕𝜕𝐺𝐺 0 for e ≤ y ≤ 2e (6) The initial concentrations for the simulations are CL = CL0 = 0.5 in the liquid phase (0.0 m ≤ y ≤ 0.5 m) and CG = CG0 = 1.0 in the gas phase (0.5 m ≤ y ≤ 1.0 m). A regular grid consisting of 150 cells in the ydirection is used for the computation. The time step follows Haroun et al.’s (2010) choice of ∆t = 4·10-5 e2/Dphys,L. 13th Urban Drainage Modelling Conference, Innsbruck (Austria), September 2025 Table 1. Case 2 – Diffusion coefficients and simulation time for each simulation. He = 0.1 He = 10 H 2 S (T = 298.15 K and T= 293.15 K) O 2 (T = 298.15 K) Dphys,G [m²/s] 10-5 10-5 1.50 · 10-5 1.98 · 10-5 Dphys,L [m²/s] 10-6 10-6 1.40 · 10-9 2.42 · 10-9 Simulation time [s] 500.000 500.000 225.000.000 225.000.000 Figure 3. Case 2 - Analytical and OpenFOAM solution produce similar results at equilibrium conditions for different Henry coefficients illustrating the accuracy of the algorithm: (a) He = 0.1, (b) He = 10, (c) H2S mass transfer at T = 298.15 K and T = 293.15 K, (d) O2 mass transfer at T = 293.15 K. Case 3: no flow, 1D diffusion across horizontal interface (dynamic mass transfer) The third test case aims at illustrating the solver’s ability to successfully fulfil the jump condition at the interface during non-equilibrium conditions. Nieves-Remacha et al (2015) used a similar setup as Haroun et al. (2010). They implemented the algorithm by Haroun et al. (2010) in OpenFOAM but also calculated a reference solution in Matlab to analyse the concentration changes in each domain during simulation. Here, the domain has the following characteristics: height: y = 0.01 m, location of interface: y = 0.005 m, He = 100, ∆t = 0.01 s. Two simulations with different diffusion coefficients in each phase have been carried out. Simulation 1 has the same diffusion coefficient in both phases, simulation 2 has different diffusion coefficients in each phase. Table 2 lists the diffusion coefficient in each phase for each simulation. 13th Urban Drainage Modelling Conference, Innsbruck (Austria), September 2025 Table 2. Case 3 – Diffusion coefficients for each simulation. DL [m²/s] DG [m²/s] Simulation 1 10-5 10-5 Simulation 2 10-7 10-5 Figures 4 and 5 illustrate the concentration profiles within the domain of the two different simulations. When comparing the concentration profiles along the domain, the OpenFOAM solver provides a good agreement with results obtained by Nieves-Remacha et al. (2015). For comparison, the results by Nieves-Remacha et al. (2015) have been retrieved using the WebPlotDigitizer tool, leading to the possibility of small uncertainties in the reference solution. They still enable a comparison between OpenFOAM solution by Nieves-Remacha et al. (2015) and their Matlab solution (both solutions were overlapping and are therefore displayed as one line) and our solution in OpenFOAM 11. The agreement between our solution and the results published by Nieves-Remacha et al. (2015) lead to the conclusion that the modular solver implementation in OpenFOAM 11 is able to produce accurate concentration profiles at different time steps. Figure 4. Case 3 – Reference (Nieves-Remacha et al., 2015) and OpenFOAM solution produce similar results during mass transfer and tracer spreading for DL = 10-5 m²/s and DG = 10-5 m²/s, illustrating the accuracy of the algorithm, illustrated times: 1 s, 5 s, 10 s, 20 s, 30 s, 40 s, 100 s. Figure 5. Case 3 – Reference (Nieves-Remacha et al., 2015) and OpenFOAM solution produce similar results during mass transfer and tracer spreading for DL = 10-7 m²/s and DG = 10-5 m²/s, illustrating the accuracy of the algorithm, illustrated times: 1 s, 2 s, 5 s, 10 s, 20 s, 50 s, 100 s. 13th Urban Drainage Modelling Conference, Innsbruck (Austria), September 2025 Conclusions and future work In this paper we have presented the new modular solver massTransferVoF, which is implemented in OpenFOAM 11. A restructuring of the OpenFOAM suite towards a modular structure deemed previously implemented solvers interH2SFoam and interO2Foam obsolete. The solver can be applied to mass transfer problems in urban water management and contains solver customizations that make the code especially suitable for H2S and O2 mass transfer problems in sewer systems (i.e. ventilation, H2S hotspots such as drop shafts). Future studies aim to further explore the role of turbulence on mass transfer rates. Apart from the relevance for the urban drainage community, this transition serves as a reference for other researchers and developers aiming to adapt legacy solvers to modern OpenFOAM architectures in general. Future work will focus on further numerical analysis, additional validation test cases, and exploring new functionalities enabled by OpenFOAM 11. Lessons learned: • Backward compatibility of custom OpenFOAM solvers prior to version 11 has limitations, resulting in the need for new modular implementation. • Integration into newer versions enables access to newest OpenFOAM features. • massTransferVoF can be used to describe temperature-dependent H2S and O2 mass transfer in sewers as well as custom Henry coefficients. Code availability The solver and setup files for the test cases described are available on github and Zenodo: - https://github.com/katharinaTe/massTransferVoF.git - https://doi.org/10.5281/zenodo.15819144 References Dixit, A. (2024). Extension and further validation of a 3D two-phase flow model for flow, transport and mass transfer in sewer systems. PhD thesis, TU Berlin, Berlin, Germany. https://doi.org/10.14279/depositonce-19857 Haroun, Y., Legendre, D. & Raynal, L. (2010) Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film. Chemical Engineering Science, 65(10), 2896-2909. https://doi.org/10.1016/j.ces.2010.01.012 Hill, C. (2020). Learning scientific programming with Python. Cambridge University Press. Hirt, C.W. & Nichols, B.D. (1981) Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics, 39, 201–225. https://doi.org/10.1016/0021-9991(81)90145-5 Hvitved-Jacobsen, T., Vollertsen, J. & Nielsen, A.H. (2013) Sewer processes: microbial and chemical process engineering of sewer networks. 2nd edn. Boca Raton, USA: CRC press. https://doi.org/10.1201/b14666 Nieves-Remacha, M.J., Yang, L. & Jensen, K.F. (2015) OpenFOAM computational fluid dynamic simulations of two-phase flow and mass transfer in an Advanced-Flow Reactor. Industrial & Engineering Chemistry Research, 54(26), 6649-6659. https://doi.org/10.1021/acs.iecr.5b00480 Rusche, H. (2003) Computational fluid dynamics of dispersed two-phase flows at high phase fractions. PhD thesis. University of London, London, United Kingdom. Sander, R. (2023). Compilation of Henry's law constants (version 5.0. 0) for water as solvent. Atmospheric Chemistry & Physics, 23(19). https://doi.org/10.5194/acp-23-10901-2023 Sharma, K. R., Yuan, Z., De Haas, D., Hamilton, G., Corrie, S., & Keller, J. (2008). Dynamics and dynamic modelling of H2S production in sewer systems. Water Research, 42(10-11), 2527-2538. https://doi.org/10.1016/j.watres.2008.02.013 Shi, T., Li, J., Ge, J., Watts, S., Wang, Y., Sharma, K., & Yuan, Z. (2025). An ODE-based swift and dynamic sewer airflow model. Water Research, 273, 123083. https://doi.org/10.1016/j.watres.2024.123083 Teuber, K. (2020) A three-dimensional two-phase model for flow, transport and mass transfer processes in sewers. PhD thesis, TU Berlin, Berlin, Germany. Available at: https://doi.org/10.14279/depositonce-9576 Teuber, K., Dixit, A., & Hinkelmann, R. (2025). CFD simulation of turbulent mass transfer of H2S and O2 in a stirring tank. Water Science & Technology, 91(1), 69-82. https://doi.org/10.2166/wst.2024.406 Wani, T.I., Rallapalli, S., Singhal, A. et al. Compounded fuzzy entropy-based derivation of uncertain critical factors causing corrosion in buried concrete sewer pipeline. npj Clean Water 8, 36 (2025). https://doi.org/10.1038/s41545-02500467-1 Yin, W. X., Lv, J. Q., Liu, S., Chen, J. J., Wei, J., Ding, C.,Yuan, Y., Bao, H. X., Wang, H. C. & Wang, A. J. (2025). Microbial-Guided prediction of methane and sulfide production in Sewers: Integrating mechanistic models with Machine learning. Bioresource Technology, 415, 131640. https://doi.org/10.1016/j.biortech.2024.131640