the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
P3DBRNS v1.0.0: a threedimensional, multiphase, multicomponent, porescale reactive transport modelling package for simulating biogeochemical processes in subsurface environments
Amir Golparvar
Matthias Kästner
Martin Thullner
The porous microenvironment of soil offers various environmental functions which are governed by physical and reactive processes. Understanding reactive transport processes in porous media is essential for many natural systems (soils, aquifers, aquatic sediments or subsurface reservoirs) or technological processes (water treatment or ceramic and fuel cell technologies). In particular, in the vadose zone of the terrestrial subsurface the spatially and temporally varying saturation of the aqueous and the gas phase leads to systems that involve complex flow and transport processes as well as reactive transformations of chemical compounds in the porous material. To describe these interacting processes and their dynamics at the pore scale requires a wellsuited modelling framework accounting for the proper description of all relevant processes at a high spatial resolution. Here we present P3DBRNS as a new opensource modelling toolbox harnessing the core libraries of OpenFOAM and coupled externally to the Biogeochemical Reaction Network Simulator (BRNS). The native OpenFOAM volumeoffluid solver is extended to have an improved representation of the fluid–fluid interface. The solvers are further developed to couple the reaction module which can be tailored for a specific reactive transport simulation. P3DRBNS is benchmarked against three different flow and reactive transport processes: (1) fluid–fluid configuration in a capillary corner, (2) mass transfer across the fluid–fluid interface and (3) microbial growth with a high degree of accuracy. Our model allows for simulation of the spatiotemporal distribution of all biochemical species in the porous structure (obtained from μCT images), for conditions that are commonly found in the laboratory and environmental systems. With our coupled computational model, we provide a reliable and efficient tool for simulating multiphase, reactive transport in porous media.
 Article
(3196 KB)  Fulltext XML

Supplement
(793 KB)  BibTeX
 EndNote
Subsurface environments (soils, aquifers, aqueous sediments) are (typically) porous media that host a multitude of biogeochemical processes and interactions and provide different versatile ecosystem functions (e.g. C sequestration; compound degradation; nutrient retention; provision of food, fibres and fuel; habitat for organisms; and water retention and purification; Baveye et al., 2016). These processes are controlled by various biological (e.g. microbial abundance and activity), chemical (e.g. distribution of dissolved and volatile species, mineral composition, and surface properties of the solid matrix) and physical (e.g. porous structure and permeability, water saturation) properties of the system. These features create a complex web of interactions, the magnitude and effectiveness of which change dynamically in space and time (Graham et al., 2014). Microbial communities, for example, and their metabolic capacity are considered to be directly related to energy and matter fluxes (Thullner et al., 2007), which are, in turn, governed by pore arrangements and their connectivity. Along with other environmental factors, this can also modify various properties of the porous media (e.g. by biomass accumulation on pore walls (Thullner, 2010) or mineral dissolution or precipitation (Meakin and Tartakovsky, 2009), which in turn are altering the conditions for biogeochemical processes, too.
In soils (or more generally the vadose zone), the dynamically varying distribution of the aqueous and gaseous phase leads to specifically complex and variable constraints for biogeochemical processes. In the past, obtaining (bio)chemical and microbiological information at the pore level was neither economically nor logistically a feasible option (Baveye et al., 2014). Also for the sake of applicability, traditionally, researchers had more tendency to look for Darcyscale solutions to tackle environmental issues (White and Brantley, 2003). The Darcyscale view (experimental, theoretical or a mixture of both) serves the purpose of practical applicability well (White and Brantley, 2003), but, for example, in the context of microbially mediated degradation processes in the vadose zone, it fails to provide insights on the driving mechanisms, as it overlooks important contributing factors, such as the tortuous porous structure/pathways open to the transport of (bio)geochemical species and the nonuniform distribution of water and air phases, as well as the nonlinear dependency between changes of the local nutrient and biomass concentrations and the bulk concentrations of (bio)geochemical species. Evidence at the microscopic level has shown that biological activity and evolution are more locally organized (Kuzyakov and Blagodatskaya, 2015) where Darcyscale studies lead to loss of crucial information. This has motivated the development of sophisticated physicsbased models implementing all aspects of hydrological, geochemical and biological processes involved in microbial growth and evolution.
Reactive transport models (RTMs) are a class of mathematical models that have been applied extensively to study biogeochemical systems for about four decades (Parkhurst and Appelo, 1999; Thullner et al., 2005; Thullner and Regnier, 2019; Meile and Scheibe, 2019). There is a long list of principal factors and mechanisms governing biogeochemical reactions at the pore scale. Numerically, these processes can be defined and solved either by fully (global) implicit approaches or by separating and solving different components one at a time. For the Darcy scale, a wide range of reactive transport models exist which allow for the simulation of biogeochemical processes (Steefel et al., 2015b). In turn, at the pore scale, models combining the simulation of flow, transport and reactive (biogeochemical) processes are scarce, and existing model developments are often driven by specific research questions and/or are subject to severe simplifications in the description of the pore space (Golparvar et al., 2021). Integrated models explicitly capturing the structural properties of the soil at the pore scale, the resulting multiphase flow and multispecies reactive transport simultaneously are hardly available (Tian and Wang, 2019).
Recently, new frontiers of porescale RTMs are emerging with the advances in computational power as well as with huge improvements in imaging techniques. The latter includes, for example, the static and dynamic scanning of porous structure as well as of fluids' distribution (Schlüter et al., 2019) or the detection of bacterial distributions in soil using the catalysed reporter deposition with fluorescence in situ hybridization (CARDFISH) technique (Schmidt et al., 2015). Direct numerical models (DNMs) are becoming the nexus of the next generation of RTMs as they represent the porous structure in a fully explicit manner (directly obtained from soil samples, digitized and fed into RTMs) in addition to offering a more flexible coupling of different components of reactive transport models (Baveye et al., 2018; Raeini et al., 2012; Li et al., 2010; Yan et al., 2016). Another advantage of using DNMs is that they offer a great deal of flexibility in considering settings and conditions that are experimentally impossible to impose (Tian and Wang, 2019).
In this work, we introduce the porescale RTM package P3DBRNS explicitly involving the structure and topology of the pore space, the coexistence/coflow of both the aqueous and the gaseous phase, the advective–diffusive transport of species in each phase, and an arbitrary set of reactive processes controlled by kinetic rate laws or thermodynamic constraints. The fluid flow field is updated via solving the Navier–Stokes (NS) equations (Patankar and Spalding, 1972). The volumeoffluid (VOF) approach is adopted to account for different phase distribution (Hirt and Nichols, 1981). The transport of chemical species is considered via solving the advection–diffusion–reaction equation (ADRE), where the concentration jump for soluble/volatile compounds across the fluid–fluid interface is modelled via the continuous species transfer (CST) method (Haroun et al., 2010). Reactive processes are defined and simulated externally via coupling the flow and transport model to the BRNS (Biogeochemical Reactions Network Solver) package (Regnier et al., 2002; Aguilera et al., 2005). The model structure is introduced, and the model performance is shown and compared with analytical counterparts. The model capabilities are depicted for a fully threedimensional case.
The entire numerical domain (Ω) can be decomposed to two main subregions: solid space (Ω_{S}) and void space (Ω_{ϑ}). The void space is further divided into aqueous phase (Ω_{ϑ,aq}) and gaseous phase (Ω_{ϑ,gs}), which are partitioned by the fluid–fluid interface (Ω_{I}). The overall domain is bounded externally between inlet (Ω_{in}) and outlet (Ω_{out}) boundaries, which allow for the inflow/outflow of different phases and chemical species, as well as noflux boundaries resembling physical walls, where nothing is allowed to leave or enter the domain (Ω_{wall}). The domain is limited internally by noflow boundaries where solid space and void space intersect (i.e. solid surface,$\phantom{\rule{0.125em}{0ex}}{\mathrm{\Omega}}_{\mathit{\vartheta}}\bigcap {\mathrm{\Omega}}_{\mathrm{S}}={\mathrm{\Omega}}_{\mathrm{wall}})$.
2.1 Fluid flow: governing equations
Evolution of a single/multiphase, isothermal, incompressible, immiscible fluid(s) can be expressed by basic conservation principles. These can be formulated into a singlefield formalism (Hirt and Nichols, 1981):
where u is the vector of velocity field; ρ is the fluid density; P is the pressure; g is the gravitational vector; and τ is the stress tensor, which can be defined as $\mathit{\tau}=\mathrm{2}\mathit{\mu}S=\mathrm{2}\mathit{\mu}\left(\mathrm{0.5}\left[\left(\mathrm{\nabla}\mathbf{u}\right)+{\left(\mathrm{\nabla}\mathbf{u}\right)}^{T}\right]\right)$, with μ as the fluid viscosity. F_{σ} denotes the interfacial tension force, which is only nonzero when two or more phases (excluding solid) are available. It is safe to use the incompressible form of the Navier–Stokes equation for low Mach and Reynolds numbers.
In the case of simultaneous flow of two different phases, their locations and distribution are represented via introducing an indicator function, α, taking values within the range [0, 1]. The first continuous fluid is marked as α=1, and the second fluid is denoted as $\mathit{\beta}=\mathrm{1}\mathit{\alpha}$, and for the transition from one fluid to the other (i.e. the interface, Ω_{I}), α varies between 0 and 1. Fluid density and viscosity in each grid cell is then calculated from a linear interpolation of this indicator function:
A mass conservative boundary condition at the fluid–fluid interface is written as
with ρ_{i} as the density of ith fluid, u_{i} as the velocity of ith fluid, w as the velocity of the interface, ${\mathrm{n}}_{{\mathrm{\Omega}}_{\mathrm{I}}}$ as the normal vector to the interface (Ω_{I}) pointing from the invading phase to the displaced one and the brackets showing the jump condition at the interface. Individual velocities, u_{i}, and the interface velocity w, are not directly calculated but furthermore averaged to derive the global mass conservation equation that is used for numerical discretization (for full derivation, consult Graveleau et al., 2017).
In the context of the finitevolume method (FVM), discretization of the physical domain produces a finite subset of discrete volumes (taking the shape of a polyhedral). The key implication of the volumeoffluid method is to define and solve for global variables, rather than having one equation for each variable in each phase. Hence, the idea is to transform the integrodifferential equations into their global versions by averaging them over each cell volume (Whitaker, 2013). For multiphase systems, after a few steps of linearization and approximation (see Hirt and Nichols, 1981, for a detailed derivation), the volumeoffluid formulation of the momentum equation (Eq. 2) is obtained as
with having $\stackrel{\mathrm{\u203e}}{\mathbf{u}}$ as the global averaged velocity vector. For the sake of simplicity, we drop the “average” notation from the global velocity vector (i.e. will refer to $\stackrel{\mathrm{\u203e}}{\mathbf{u}}$ as u), for the rest of this paper. The pressure gradients (and Reynolds numbers) considered in our simulations are in the range that render changes in the gas compressibility negligible.
Since we are dealing with only two fluids, index i in u_{i} takes only two values; α and β – one for each phase. A global, mass conservative, advection equation is used to describe the evolution of the indicator function:
where α indicates the volume fraction of phase 1, and ${u}_{\mathrm{c}}={u}_{\mathit{\alpha}}{u}_{\mathit{\beta}}$ is the vector of the compressive velocity, with u_{α,β} as the velocity vector of phase α and β right on the edge of the interface (detailed explanation on deriving Eq. 7 can be found in the Supplement). It is derived from the mass conservation equation written for phase α, which computationally helps with maintaining stiffness/sharpness of the interface. Sharpening the interface means having the interface span fewer computational grids. u_{c} is the vector of compressive velocity on each face of all computational grids. Since we do not solve for the velocity field of each phase individually, a direct calculation of u_{c} is not possible. However, we can rather take an indirect approach for computing u_{c} as follows:
In Eq. (8), c_{α} is a compression coefficient providing some level of control over how wide the interface spans. The max function operates on the magnitude of unit velocity vector calculated on all the faces of a computational grid. To counteract the numerical diffusion and avoid the spread of interface over several computational grids, values of c_{α}>1 provide an enhanced/sharper interface, whereas a value of c_{α}=0 gives no compression of the interface (Graveleau et al., 2017). In simulation scenarios introduced in this paper, c_{α} has been assigned the value of 1, unless stated otherwise.
To calculate the interfacial tension force, F_{σ}, Brackbill et al. (1992) have introduced a continuum surface force (CSF) which requires computation of the curvature of the interface:
with having κ as the mean interface curvature in each computational grid and ${n}_{{\mathrm{\Omega}}_{\mathrm{I}}}$ as the interface unit normal vector, defined as
Given the curvature, the interfacial tension force can be approximated as
where σ is the surface tension between two fluids (derivation of this approximation can be found in Brackbill et al., 1992).
2.2 Reactive transport: governing equations
Concentrations of mobile species are affected by advection (i.e. transport with the moving fluid), molecular diffusion and reactive transformation. Also, in the case of having two fluids simultaneously in the system, different species can cross the fluid–fluid interface, causing local fluctuations in concentration values. In general to account for all the changes in species concentrations, the ADRE for biogeochemical reactive components can be written as
where J_{d,i} is the molecular diffusive flux of component i, J_{m,i} is the mass flux of component i due to mass transfer across the fluid–fluid interface and R_{i} accounts for the changes in concentration of component i due to reactions. Molecular diffusion follows Fick's law:
where D_{i} is the diffusion coefficient of species i. At the interface, the assumption of thermodynamic equilibrium implies equality of chemical potentials. Given the condition that liquid concentration of component i is proportional to the partial pressure of the species in the secondary phase (e.g. gas, oil or minerals), a partitioning relationship such as Raoult or Henry's law (Danckwerts and Lannus, 1970) can be established to relate species concentrations on both sides of the interface:
with C_{i,α} as the concentration of species i in phase α, H_{i} as Henry's constant of species i and C_{i,β} as the concentration of species i in phase β. Depending on if a given compound's concentration in the aqueous phase or the gaseous phase is multiplied by Henry's coefficient (Eq. 14), the definition of Henry's constant switches between the solubility or volatility for that compound (i.e. ${H}_{i}^{\mathrm{soluility}}=\phantom{\rule{0.125em}{0ex}}\mathrm{1}/{H}_{i}^{\mathrm{volatility}})$ (Sander, 2015). Unless otherwise stated, the volatility concept of Henry's law is adopted in order to define the concentration relationship of a given compound across the fluid–fluid interface (Eq. 14). The concentration field around the fluid–fluid interface (where ∇α≠0) at equilibrium, for any values of H≠1, is discontinuous, which imposes the additional flux, J_{m,i}, to satisfy the concentration jump across the interface. Hence the mass transfer flux, J_{m,i}, can be derived within the VOF framework (i.e. CST) as follows (Haroun et al., 2010):
It is noteworthy that few assumptions and volume averaging methods are implemented to derive Eq. (15); readers are encouraged to check the references for more details. The diffusion coefficient is calculated from harmonic interpolation:
where ${D}_{i,\mathit{\alpha}\mathit{\beta}}$ is the diffusion coefficient of species i in phase α and β respectively.
Simulated reactions include kinetically as well as thermodynamically constrained reactions. For a kinetically constrained reaction j, the reaction rate ${r}_{j}=f({C}_{\mathrm{1}},\mathrm{\dots},{C}_{n})$ is needed, while for a thermodynamically constrained reaction k, the equilibrium conditions defined by a law of mass action ${M}_{k}=f({C}_{\mathrm{1}},\mathrm{\dots},{C}_{n})$ are needed, with M_{k} as the equilibrium constant. These equations can be of arbitrary form, and the resulting reaction network defines the term R_{i} in Eq. (12) (Aguilera et al., 2005; Regnier et al., 2002). For immobile species, concentration changes are only due to reactive processes.
2.3 Boundary conditions (BCs)
There are various types of boundary conditions, corresponding to the real physical conditions, most of which can be derived from two basic types:

the Dirichlet boundary (fixed value), which relates the value of a variable at a given geometric location to a constant value, for example, C_{i}=1M in Ω_{in}, meaning a constant 1 molar concentration of component i at the boundary,

the von Neumann boundary (fixed gradient), which provides the value of a variable's gradient at the face of the boundary cell, for example, ${\partial}_{n}\mathbf{p}=\mathrm{0}$ in Ω_{wall}, giving a zero pressure gradient on the wall.
In general, our model can apply any of these basic boundary conditions to any scalar or vector variables such as pressure, velocity/flux and concentration of volumeoffluid fields, but one needs to ensure that the imposed BC(s) are both compatible and that they reflect the correct physical boundary conditions. For example, for velocity/flux–pressure coupling, a Dirichlet (i.e. constant) boundary for flux at the inlet can be coupled with either (1) fixed discharge velocity/flux and zerogradient (i.e. von Neumann) pressure at both the inlet and outlet, (2) a constant pressure head at the inlet and atmospheric pressure at the outlet with zerogradient velocity/flux at both ends, or (3) fixed values of pressure and velocity/flux at one end and zero gradient at the other end. In the beginning of Sect. 2, the typical composition and configuration of an arbitrary computational domain are described. Inlet, outlet and impermeable boundaries are amongst the most common types that one might face. Inlet BC means for the direction of fluid flux to be pointing inwards (i.e. into the domain), while for the outlet, the direction of the flux should be outwards. Also for the impermeable wall, zeroorthogonal fluxes need to be satisfied. Either Dirichlet, von Neumann or a mixture of both can be used at a particular boundary. Mathematical translation and implementation of these boundaries are provided in the next section. Timedependent BCs (e.g. cyclic or seasonal water/species influx) are also readily available to be applied but have never been used in this work. Unless otherwise stated, boundary conditions that have been imposed on each section of the computational domain are described as follows.
At impermeable boundaries (Ω_{wall}). The physical wall implies no flux perpendicular to the normal vector to its surface. Noslip BC is an appropriate BC for the velocity field on the wall. In general, they can all be written as
For the velocity field, on the wall, a slip boundary condition is also available to be applied.
At inlet–outlet boundaries (Ω_{in}Ω_{out}). The concentrations of reactants, products and inert tracers are set to fixed values at inlet, while they are allowed to leave the domain at the outlet with zerogradient boundary conditions. A constant flow rate with zeropressure gradient is applied at the inlet, and an atmospheric pressure (fixed value) with zero velocity gradient is set at outlet. Also, in the case of twophase flow, the invading phase is set to enter from the inlet at a fixed value and exits from the outlet with zerogradient BCs. Mathematically, they can be expressed as
together with
with u_{n} as the normal velocity vector. While we have mostly applied Eqs. (18) and (19) for designing an inlet–outlet duo, other formats, such as defining a pressure head (and zerogradient velocity) on the inlet in combination with either constant exit pressure or constant discharge rate, are readily available to implement as well.
At the fluid–fluid–solid contact line $\left({\mathrm{\Omega}}_{{\mathrm{l}}^{\mathrm{3}}}\right)$. At the fluid–fluid–solid contact line, in the case of no interactions or no reaction of any chemical species with the solid, the boundary condition at the triple point is derived to be
with n_{s} as the normal vector to the solid surface (Graveleau et al., 2017). Also, the concept of contact angle is applied by making the following modification to the interface normal vector:
where n_{s} is the normal vector, and t_{s} is the tangential vector to the solid surface (Brackbill et al., 1992). At the triple point, i.e. fluid–fluid–solid interface, ${\mathrm{n}}_{{\mathrm{\Omega}}_{{\mathrm{l}}^{\mathrm{3}}}}$is used for normal vector to the interface. CSF, though, has been reportedly generating nonphysical spurious currents (Scardovelli and Zaleski, 1999). For this, many have tried to eliminate/mitigate this issue by explicit representation of the interface either via using the geometric VOF method (Popinet, 2009) or coupled levelset (LS) VOF functions (Albadawi et al., 2013). Geometric VOF is quite suitable for structured grids, but for porous structures with highly unstructured grids, the calculations can become quite complicated. Alternatively, Raeini et al. (2012) suggested filtering the capillary forces parallel to the interface, which can significantly reduce the nonphysical velocities. In short, the modifications they proposed and which are used here are (1) smoothing the indicator function to have a better measure of the interface curvature, (2) sharpening the indicator function for computation of the interfacial tension force, (3) filtering the capillary pressure force parallel to the interface and (4) filtering capillary fluxes based on the capillary pressure gradient (for a full description of each point, please consult with Raeini et al., 2012). The correction introduced for filtering capillary forces helps with eliminating some of the parasitic velocities parallel to the interface. To sum up what has been presented so far, we integrated (a) the original interFoam solver from the OpenFOAM library that only solves for the advection–diffusion transport of two phase flow with (b) the improvedinterfaceresolver library from Raeini et al. (2012) and (c) added a scalar transport solver on top of them. Finally, the fullscale advection–diffusion–reaction model of the biogeochemical species is attained by coupling this to an external reaction network solver, which is explained in Sect. 2.4.
2.4 Numerical formulation
The mass conservation (Eq. 1), momentum (NS – Eq. 2) and indicator function (Eq. 7) equations are all implemented within the opensource computational fluid dynamics (CFD) package, OpenFOAM (Greenshields, 2015). OpenFOAM utilizes the finitevolume methodology (FVM), a common choice for CFD problems as FVM works only with conservative flux evaluation at each computational cell's boundaries, making it robust in handling nonlinear transport problems. Also all the differential equations mentioned before are first written in their integral form over each cell volume and then converted to the surface summations using Green's theorem.
The original twophase (VOF) flow solver, i.e. interFoam, is modified to construct our biogeochemical reactive transport package. The momentum equation (Eq. 2) is linearized in a semidiscrete form as
where A_{d} holds the diagonal elements of the coefficient matrix, H(u) contains offdiagonal elements of the coefficient matrix including all source terms and F entails any body forces (interfacial tension force only in this case). Temporal discretization is handled via the firstorder Euler method, while spatial discretization is managed via secondorder finitevolume schemes. Convection terms of the momentum equation and indicator function (Eq. 7) are computed using a bounded selffiltered central differencing (SFCD) scheme (based on Gauss's theorem). Rearranging Eq. (22) for velocity and imposing the continuity equation (Eq. 1), the following linear pressure equation can be obtained:
S_{f} in Eq. (23) denotes the outward area vector of face f, the notation ${\mathrm{\nabla}}_{f}^{\perp}$ shows face normal gradients calculated right on the face centres, 〈〉_{f} shows the interpolated values of a facecentred parameter from its cellcentred counterpart and φ_{F,f} is the interfacial force flux term.
The velocity–pressure coupling of Eqs. (1) and (2) is solved using Pressure Implicit with Splitting of Operators (PISO) (Issa, 1986). PISO embodies a predictor–corrector strategy to simultaneously update pressure and velocity within each time step. The resultant system of equations is solved on the cell faces and then interpolated back to calculate velocities and pressure at the cell centres. The coupling of the indicator function (Eq. 7) and momentum equation is explicitly defined and solved right after the PISO step is finished. Within the same time step, transport and reaction of different species are then solved sequentially – using a sequential noniterative algorithm (Steefel et al., 2015a; Steefel and MacQuarrie, 1996). Time step size is controlled by introducing a Courant number. Time is discretized using either Euler or Crank–Nicolson methods, and spatial discretization is performed using the van Leer secondorder total variation diminishing scheme (TVD) (van Leer, 1979).
The reaction network is built separately and externally solved within the BRNS package – which employs firstorder Taylor series expansion terms and uses the Newton–Raphson method to iteratively solve the system of linear equations (Regnier et al., 2002). BRNS utilizes MAPLE programming language to construct the Jacobian matrix (which contains the partial derivatives of unknown parameters, i.e. concentrations) and other problemrelated data such as rate parameters and translating them to a FORTRAN package. The FORTRAN code is then compiled to generate shared objects (*.so file) that can be dynamically called later from the transport solver (Centler et al., 2010). The significance of having dynamically shared object files is more apparent when running computationally demanding cases/scenarios while decomposing and running the application in parallel. BRNS is invoked once the new concentrations are computed from the transport solver. The updated concentrations from the BRNS library (i.e. updating concentrations from redox reactions) are then fed back into the transport solver before moving to the next time step. This process repeats until the final time is reached. This coupling scheme has been successfully used for other RTM approaches before (Centler et al., 2010; Gharasoo et al., 2012; Nick et al., 2013). As the reactions are localized, the reaction solver is modular, and OpenFOAM inherently provides parallelized simulations (via domain decomposition). The P3DBRNS can easily be used to model larger systems. To achieve higher performance, it is recommended to utilize physical cores rather than hyperthreading. The parallelization of our model strengthens its scalability in the sense of the size (pore scale or Darcy scale) of the simulated system. However, in terms of upscaling (e.g. from the pore scale to the Darcy scale), an intermediate step would be required depending on the complexity of the processes that are involved and on the size of the domain.
Prior to run simulations, the physical settings of the domain are required to be specified; i.e. the physical geometry of the pore space with proper boundaries and the meshing scheme should be designed. OpenFOAM provides a basic utility for defining boundaries as well as mesh generation which are translatable by the OpenFOAM engine. Any other meshing software/freeware can be freely used as long as an OpenFOAMcompatible format for the meshed file can be created. The overall workflow required to build and run a case/scenario is summarized in Fig. 2.
The presented reactive transport model is designed (1) to capture realworld pore structures in up to three dimensions, (2) to explicitly simulate the transient distribution of a gas and a liquid phase within the entire pore space, and (3) to simulate a full set of advection–diffusion–reaction mechanisms. Capturing the correct curvature of the interface depends heavily on the grid resolution. For a fixed velocity magnitude, a higher grid resolution enforces a shorter time step size (from Courant number) for the numerical simulations to converge. Also, to validate different features of the model various, simplified scenarios were used which allow the use of analytical expressions as reference for the numerical results. Here we show three representative test scenarios addressing different features of the model (twophase flow, mass transfer across the fluid–fluid interface and reactive transport) individually. Subsequently, the model capabilities are depicted in a final biodegradation scenario, making use of the various model features simultaneously.
3.1 Fluid configurations
In order to test our model's performance in simulating twophase flow, we have zoomed into a twodimensional porous structure and isolated only one single corner, taking the shape of an equilateral triangle. A triangulated mesh is adopted that naturally conforms to the overall shape of the domain. Initially, two immiscible fluids (one wetting and one nonwetting, for example, water and air) are placed in such a way that their interface forms a straight line (Fig. 3a). The side length of the triangle is 1 mm with a mesh size of 1 µm. Under thermodynamic equilibrium conditions, the force exerted by the pressure difference between two phases is countered by the interfacial tension force. This, along with the contact angle of the nonwetting phase at the wall surface in presence of the wetting phase (e.g. water), determines the topology of the fluid–fluid interface. For a given corner halfangle, the distance that the wetting phase spreads over the solid surface from the corner vertices (the highlighted section in green in Fig. 3b), b, can be calculated as
with r as the radius of the interface's curvature, θ as the contact angle and β as the corner halfangle (Blunt, 2017). In order to reach thermodynamic equilibrium, we performed transient, twophase flow simulations to compute velocity, pressure and indicator function fields until the triple contact line $\left({\mathrm{\Omega}}_{{\mathrm{l}}^{\mathrm{3}}}\right)$ is static. For this, we first divided the equilateral triangle in half, as the problem is symmetrical along the height of the triangle. The symmetrical plane implies that there is no gradient (of any scalar or vector field) perpendicular to its surface, while the tangential components (of all fields) remain the same. To find the fluid configuration at equilibrium, we simulated the twophase flow scenario in two steps. First, we applied a closedboundary condition on the bottom domain by setting u=0 together with ${\partial}_{\mathbf{n}}p=\mathrm{0}$. Also a closed boundary is imposed on the topmost part of the domain, which follows the same BC as the bottom. This way, the interface is able to reconfigure and reorient itself in order to recreate the imposed contact angle with the wall, and at the same time, pressure is allowed to build up in both phases to support the shape of the interface. Then, in order to obtain an equilibrium curvature for the interface, bottom and top domains are opened. This is achieved by setting the (1) pressure in Ω_{in} to the average pressure within the nonwetting phase and (2) pressure in Ω_{out} to the average pressure within the wetting phase together with (3) ${\partial}_{n}\mathbf{u}=\mathrm{0}$ on both Ω_{in} and Ω_{out}. At this stage, we applied a special BC for the indicator function to allow the fluids to enter or leave the domain at both ends, so that the interface could freely transition to its static shape. At the inlet (Ω_{in}), the BC for α is set to switch between ${\partial}_{n}\mathit{\alpha}=\mathrm{0}$, if the fluid flux is pointing outwards, and α=0 if the fluid flux is directed into the domain. Also at the outlet (Ω_{out}), the BC for α switches between α=1, if the fluid flux is inwards, and ${\partial}_{n}\mathit{\alpha}=\mathrm{0}$, if the fluid flux is outwards. This ensured that appropriate fluids entered the domain from either inlet or outlet boundaries. The radius of curvature can be also evaluated from the Young–Laplace equation (${P}_{\mathrm{c}}=\frac{\mathit{\sigma}}{r})$. With a pressure difference of 255.33 (kg m^{−1} s^{−2}) obtained from the last step and a surface tension of 0.07 (kg s^{−2}), the radius of curvature is calculated to be $\mathrm{2.17}\times {\mathrm{10}}^{\mathrm{4}}$ m. In a different approach, once the interface attains stationarity, we calculated r for Eq. (24) as the reciprocal of the interface's mean curvature ($\mathrm{2.77}\times {\mathrm{10}}^{\mathrm{4}}$ m). For a contact angle of 10^{∘} and a corner halfangle of 30^{∘}, the analytically calculated value for the length, b, of the section in contact with the wetting phase is 375 µm, while the numerical solution yields 370 µm. With a relative error of 1.21 %, this shows a reasonable match between numerical and analytical solutions in modelling twophase flow.
3.2 Mass transfer across the fluid–fluid interface
Mass transfer of dissolved species between different phases is particularly of importance for various biogeochemical processes in unsaturated subsurface environments.
Model performance in simulating mass flux across the fluid–fluid interface is validated via a numerical experiment in which two immiscible stationary fluids (an aqueous – α – and a gaseous – β – phase, u=0 in Ω_{ϑ}) are horizontally (to remove buoyancy effects) residing in a onedimensional tube of 10 mm length with mesh size of 100 µm. The general partial differential equation (PDE) of Eq. (12) takes the form of a simple diffusive transport as
with C_{tr} as the concentration of a volatile tracer and D_{tr,i} as the diffusivity of the tracer in phase i. Each phase is set to occupy half of the total volume (Fig. 4a). The system is initialized with a volatile chemical species of concentration of 1 mol m^{−3} in Ω_{ϑ,aq} and 0 mol m^{−3} in Ω_{ϑ,gs}. At the inlet and the outlet boundary, tracer concentration equals that of the nearest solution such that, in a short simulation time, it yields no concentration gradient into or out of the domain. The tendency of the dissolved chemical component to cross the fluid–fluid interface is expressed using a constant Henry coefficient. Tracer diffusivity is set to be $\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$ m^{2} s^{−1} in both phases. The analytical solution for Eq. (25) can be found in Bird (2002).
Three scenarios with low, neutral and high affinity of the volatile compound towards the gaseous phase are considered with corresponding Henry coefficients of 0.01 (lowvolatility, similar to naphthalene), 1 (moderatevolatility, for example, vinyl chloride) and 100 respectively (highvolatility, for example, heptane). For a low value of H (H=0.01 – Fig. 4b), little (almost no) tracer is crossing the interface, while at neutral condition (H=1 – Fig. 4c), tracer diffusion is invariant to the phase it is occupying. Evidently for high values of H (H=100 – Fig. 4d), significant reduction in the tracer amount within the liquid phase can be detected, which is accompanied by a notable change in concentration across the interface. This complies fully with the concentration jump for such highly volatile components between the liquid and the gas phase. The numerical results are ubiquitously identical to the results of the analytical solution (Fig. 4b, c, d). The effect of grid size/resolution is also investigated for this scenario. With 10 times higher grid resolution, the total CPUelapsed time is increased from ∼650 to ∼3500 s. The concentration profile remains unchanged, but, the average residual of the numerical solution of the concentration field, calculated at the end of the simulation, is increased from $\mathrm{3.1}\times {\mathrm{10}}^{\mathrm{10}}$ to $\mathrm{6.5}\times {\mathrm{10}}^{\mathrm{10}}$ (meaning increasing resolution does not necessarily help with numerical convergence).
3.3 Microbial growth and reactive transport
Our modelling framework can parameterize any type of reactions; we put the main focus in this subsection on microbially driven redox transformations (i.e. a type of reactions commonly encountered in soils and other porous media environments) and on the implementation of the corresponding mathematical formulation. To validate our model with a scenario in which bacterial biomass is allowed to evolve (i.e. to grow and to decay), we adapted a conceptual biodegradation scenario from Cirpka and Valocchi (2007) in which a fully watersaturated, twodimensional channel is subjected to a constant flux of two different components, ED (electron donor, for example, hydrocarbon) and EA (electron acceptor, for example oxygen). The imposed uniform flow field is assumed to be constant over time and only has the x component. The bacteria residing in the channel facilitate the reaction between ED and EA, which can be written in an abstract form as ${f}_{a}\mathrm{ED}+{f}_{b}\mathrm{EA}\stackrel{\mathrm{biomass}}{\u27f6}\phantom{\rule{0.125em}{0ex}}{f}_{c}$Prod, where biomass is the microbial biomass, Prod is the product(s) (e.g. metabolites such as carbon dioxide), and f_{a}, f_{b} and f_{c} are stoichiometric coefficients. Assuming a doubleMonod kinetics for expressing microbial growth and the microbially driven reaction rates, as well as assuming none of the reactants nor products are involved in secondary reactions, the ADRE (Eq. 12) for each chemical species can then be written as
where u is the velocity (which has only a constant xcomponent); D_{t} is the transverse dispersivity; C_{ED}, C_{EA}, C_{Met} and C_{bio} are concentrations of ED, EA, Prod and biomass respectively; K_{ED} and K_{EA} are halfsaturation constants for respective compounds in the biomass growth term; Y is the yield coefficient; μ_{max} is the maximum bacterial growth rate; and λ is the bacterial decay rate. Using these equations, Cirpka and Valocchi (2007) developed an analytical solution for steadystate conditions, which in the version of Cirpka and Valocchi (2009) is used as reference for the numerical results.
The numerical experiment is designed to have ED and EA, occupying 25 % and 75 % of the inlet respectively, and, simultaneously, invading the domain under a constant uniform velocity field, with concentration of ${C}_{\mathrm{ED}}^{\mathrm{inlet}}$ and ${C}_{\mathrm{EA}}^{\mathrm{inlet}}$. In a realworld scenario, this can be seen as a plume of a contaminant (i.e. a hydrocarbon as ED) being carried into the domain within an oxygenated stream, and essentially we are interested in knowing the final concentration/distribution of all biochemical species within the domain. The parameters used in this scenario are summarized in Table 1. Transient reactive transport simulations are performed until a steady state is achieved. For validation, we analyse all concentration profiles along the y axis at a fixed distance of x=2 m and compare them with the analytical solutions. The analytical and numerical results show an almost perfect agreement (Fig. 5b–e).
3.4 Theme: demonstrating model capabilities
The scenarios described above are designed to serve as the sole purpose of creating a baseline for validating the numerical toolbox – simple enough where analytical solutions could exist. Unlike the simplicity introduced in previous sections, simulating soil processes with all of the complexities, though, would require having all the modelling elements to be present. We thus present here a scenario with an unsaturated soil hosting the facultative anaerobic bacteria Agrobacterium tumefaciens, which performs aerobic respiration under oxic conditions but switches to denitrification using nitrate, nitrite or nitric oxide under anoxic conditions (Kampschreur et al., 2012). This example allows us to show our model capabilities, as it involves (1) the actual microstructure of the soil, (2) unsaturated conditions and (3) an enzymatic reaction network with limiting/inhibition terms. The microstructure is obtained via subsampling from a larger µCT image with a voxel size of 6 µm (see Supplement). A twophase simulation is then performed on the voxelized subsample to obtain the fluxes and phase distribution of air and water within the pore space. For this, the entire domain is initially filled with water and subject to injecting air from the top boundary with a constant flux of 0.013 mL h^{−1}. An important note to make here is with a relatively high influx, advection transport acts as the bottleneck for numerical time steps. Hence, reactions are performed at quite a slower pace (i.e. larger time steps roughly estimated around 10 h). This separation of processes helps improve the overall run time of the simulations. Generally, the time step sizes are automatically enforced by the Courant number from the transient advective–diffusive transport equation (with order of 10^{−5} s). The biomass is assumed to be nonmotile, meaning it sticks to the solid surface and shows no planktonic behaviour. Fluids are allowed to leave the domain from the bottom part (kept at atmospheric pressure), while all the remaining sides are set to be impermeable walls. Once fluid configurations in the domain are stationary, their distribution, along with the velocity profile, is used as a basis for the reactive transport simulations (phase distributions can be found in the Supplement). Using succinate (${\mathrm{C}}_{\mathrm{4}}{\mathrm{H}}_{\mathrm{4}}{\mathrm{O}}_{\mathrm{4}}^{\mathrm{2}}$) as the organic carbon substrate to be degraded, a metabolic reaction network is constructed with four microbial degradation pathways each following Monodtype kinetics: (1) aerobic respiration with a nitric oxide (NO) inhibitory term; (2) nitrate (NO${}_{\mathrm{3}}^{}$) reduction; (3) nitrite (NO${}_{\mathrm{2}}^{}$) reduction; and (4) NO reduction, with having oxygen (O_{2}) as the inhibitory element for all denitrification conversions (Eq. 27). Also three additional equations are considered for the synthesis of the three different enzymes required for degradation processes (Eq. 28). We consider only one single strain of bacteria (Agrobacterium tumefaciens) which has the benefits of performing both aerobic respiration and denitrification. Bacteria are considered to be nonmotile with an initial concentration of 0.25 mol m^{−2} and uniformly covering the entire grain surface area. Succinate has its initial concentration in the aqueous phase set at 0.2 mM (0 mM in the gaseous phase), while all other species have their initial concentrations of 0 mM in both aqueous and gaseous phases. The boundary condition for all concentration fields on all boundaries is set to zero gradient, except for the inlet boundary (fully saturated with air), where for oxygen it is set to 0.03567 mM, and for all others it is set to 0 mM. In order to avoid depletion of the nitrate in the system, a nitrate concentration of 0.1 µM (as initial condition) is provided. The complete reaction network can be written as follows (Kampschreur et al., 2012):
Several assumptions are made for preparing the kinetics of the reactions: (1) reaction rates are limited by the maximum specific uptake rate of succinate and are hence independent of its concentration (Beun et al., 2000); (2) a sufficient amount of buffer is added to the solution to keep the pH level constant; (3) three nitrogen reductase enzymes (ξ_{sat, NOR} for NO reduction, ξ_{sat, NIR} for nitrite reduction and ξ_{sat, NAP} for nitrate reduction) can have saturation values varying between 0 (i.e. nonexisting) and 1 in a bacterial cell; and (4) inhibitory oxygen limits the reduction of NO, ${\mathrm{NO}}_{\mathrm{2}}^{}$ and ${\mathrm{NO}}_{\mathrm{3}}^{}$. Reaction rates are designed to have a dependency on the enzymes' level and biomass concentration with proper limiting/inhibiting terms. Equation (12) is used to describe the evolution of each biochemical species. The final system of advective–diffusive–reactive equations is adapted from Kampschreur et al. (2012):
The full list of modelling parameters used for this study can be found in the Supplement.
Reactive transport simulations were performed until a quasisteady state was achieved. This was characterized by all chemical species concentrations reaching a steady state, as determined by the degradation activity of the given distribution of microorganisms. Since microbial growth takes place at much larger timescales than the porescale transport processes, no significant growth takes place during the simulated time period, and results shown are nearly identical to the initial conditions. Simulation results show that the presence of air in this twophase system affects the distribution of biochemical species. Air, in the nonwetting phase, occupies the central part of the pore space, while the aqueous phase is expected to cover the corners and crevices (Fig. 6a). For oxygen with ${H}_{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{31}$, a higher concentration is observed in the air compared to that of the adjacent aqueous phase (Fig. 6d). An analysis of how the volatility of a tracer compound may affect its residence time in the porous medium is given in the Supplement. Since the biomass is only present on the grain surfaces (Fig. 6b), oxygen, nitrate and succinate deplete as the microbially mediated reactions only occur at these microlocations. Fresh oxygen and nitrate thus need to diffuse from the bulk (either from the aqueous phase or the air) to the reactive sites. The regions with high (i.e. not degraded) succinate concentrations are compatible with low concentration regions of oxygen and nitrate; i.e. the reactions are limited by the bioavailable oxygen and nitrate (Fig. 6b–e). Finally, all three enzymes have an increased abundance in anaerobic regions with an active biomass (saturation map of nitrate reductase enzyme is shown in Fig. 6e). While the saturation of nitrate reductase enzyme grows linearly with time (until 0.25 s), the rate at which the nitrite and NO reductase enzymes (ξ_{sat, NIR} and ξ_{sat, NOR} respectively) grow is rather slow for the very beginning of the simulation (until ∼0.2 s), but it surges exponentially afterward. A spatially integrated assessment of the degradation processes showed that for the presented example, 99 % of the total succinate degradation is attributed to aerobic respiration, while a trivial amount is attributed to the three anaerobic processes (nitrate reduction, nitrite reduction and NO reductions).
The presented results highlight the ability of the model to combine a highresolution simulation of multiphase flow and transport processes with the simulation of complex biogeochemical processes. This allows for a realistic simulation of the porescale distribution of reactive processes and for the derivation of an accurate aggregated description of these processes.
As it can be seen from Fig. 6, our model can be used (among other options) to identify clusters in which succinate is most and least depleted. This would ease the process of analysing the results by isolating the parameters that are boosting/limiting the degradation of the carbon source. A 3D visualization of the oxygen and succinate distributions can be found in Golparvar (2022).
In this paper, we have presented a newly developed modelling framework for simulating reactive transport processes in real porous soil structures obtained from µCT images under unsaturated conditions. The successful application of various benchmark test showed the model's accuracy in the simulation of (1) the movement of water and air phase in variably saturated conditions via the enhanced algebraic volumeoffluid method (Raeini et al., 2012) coupled with the Navier–Stokes equation; (2) the transport of different species in both phases by the full advective–diffusive transport equation; and finally (3) use of the operator splitting technique, an arbitrary set of biogeochemical reactions solved externally by the Biogeochemical Reaction Network Simulator and communicated back into the main solver.
The presented model provides a novel and unique combination of porescale simulations of twophase flow, transport of dissolved and volatile species, and their reactive transformations. This makes it an accurate and powerful tool for the simulation of soil systems or other unsaturated porous media and of the reactive transport processes therein. While developed with the aim for simulating biogeochemical processes in soils, the model is equally applicable for simulating other abiotic reactive processes coupled to the dynamics of flow and transport in variably saturated pore structures of arbitrary geometry. Our modelling framework is properly designed for simulating biogeochemical processes such as carbon–nitrogen–sulfur–phosphorus cycles in soil as well as mixing and migration of contaminants in both unsaturated soil and water aquifers. It comes with the benefit of explicit recognition of the soil structure (i.e. using the 3D structure as close to the original shape as possible with the least number of simplifications/modifications), phase dynamics/distributions, and the capability of designing the complete redox reactions necessary for a given process in a straightforward fashion. It is best suitable for running prepilot tests as feasibility scenarios where the stakes for the success of the project are high. Also, our model provides the best tool for designing hypothetical experiments that are hard (if not impossible) to implement experimentally (e.g. a specific distribution of biomass–reactants within the domain, or variation of specific properties of reactive compound and/or the porous matrix). Furthermore, the highresolution modelling results provided by this model support the upscaling of reactive transport process description from the pore to the Darcy scale and from the process to the observation scale respectively.
Although the current version of our numerical model is already covering a wide range of biophysiochemical properties of the soil constituents, for having a more realistic representation of multiphase, multicomponent reactive transport in partially saturated porous media, a few more factors still might be considered in future developments of the model: (1) shrinkage/expansion of the air/aqueous phase due to mass transfer of chemical species across the fluid–fluid interface; (2) the accounting of gas compressibility by adding an equation of state for tracking changes of air volume/density under flowing conditions; (3) translation of accumulated biomass on the grain surfaces into new flowresistance components, which would potentially change the velocity streamlines (i.e. bioclogging); (4) changes of the grain surface structures and of the associated solid–liquid interface due to mineral precipitation/dilution or due to accumulation/depletion of solid organic material; (5) chemotactic behaviour of the microbial species; and (6) osmotic forces and electromigration. Due to the modular structure of P3DBRNS, such features can be relatively easily included in future upgrades of our model.
The source codes, benchmark and demonstration cases, along with instruction for installing and running each case that are presented in this paper, are archived at https://github.com/amirgolp/P3DBRNS (last access: 27 February 2022); https://doi.org/10.5281/zenodo.6301317 (amirgolp, 2022).
The supplement related to this article is available online at: https://doi.org/10.5194/gmd178812024supplement.
AG was responsible for model/software curation, validation and visualization. Conceptualization and methodology development were managed by AG, MK and MT. The writing of the original manuscript was handled by AG, while all authors contributed to the revision and curation of the final draft. The entire work is supervised by MT.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This work was funded by the Helmholtz Association via the integrated project “Controlling Chemicals' Fate”.
This research has been supported by the Helmholtz Association (Controlling Chemicals' Fate).
The article processing charges for this openaccess publication were covered by the Helmholtz Centre for Environmental Research – UFZ.
This paper was edited by Richard Neale and reviewed by two anonymous referees.
Aguilera, D. R., Jourabchi, P., Spiteri, C., and Regnier, P.: A knowledgebased reactive transport approach for the simulation of biogeochemical dynamics in Earth systems, Geochem. Geophy. Geosy., 6, Q07012, https://doi.org/10.1029/2004gc000899, 2005.
Albadawi, A., Donoghue, D. B., Robinson, A. J., Murray, D. B., and Delauré, Y. M. C.: Influence of surface tension implementation in Volume of Fluid and coupled Volume of Fluid with Level Set methods for bubble growth and detachment, Int. J. Multiphas. Flow, 53, 11–28, https://doi.org/10.1016/j.ijmultiphaseflow.2013.01.005, 2013.
amirgolp: amirgolp/P3DBRNS: v0.1.0beta (v.1.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.6301317, 2022.
Baveye, P. C., Palfreyman, J., and Otten, W.: Research Efforts Involving Several Disciplines: Adherence to a Clear Nomenclature Is Needed, Water Air Soil Pollut., 225, 1997, https://doi.org/10.1007/s1127001419977, 2014.
Baveye, P. C., Baveye, J., and John Gowdy, J.: Soil “Ecosystem” Services and Natural Capital: Critical Appraisal of Research on Uncertain Ground, Front. Environ. Sci., 4, 41, https://doi.org/10.3389/fenvs.2016.00041, 2016.
Baveye, P. C., Otten, W., Kravchenko, A., BalseiroRomero, M., Beckers, É., Chalhoub, M., Darnault, C., Eickhorst, T., Garnier, P., and Hapca, S.: Emergent properties of microbial activity in heterogeneous soil microenvironments: different research approaches are slowly converging, yet major challenges remain, Front. Microbiol., 9, 1929, https://doi.org/10.3389/fmicb.2018.01929, 2018.
Beun, J., Verhoef, E., Van Loosdrecht, M., and Heijnen, J.: Stoichiometry and kinetics of polyâhydroxybutyrate metabolism under denitrifying conditions in activated sludge cultures, Biotechnol. Bioeng., 68, 496–507, 2000.
Bird, R. B.: Transport phenomena, Appl. Mech. Rev., 55, R1–R4, 2002.
Blunt, M. J.: Multiphase flow in permeable media: A porescale perspective, Cambridge University Press, https://doi.org/10.1017/9781316145098.012, 2017.
Brackbill, J. U., Kothe, D. B., and Zemach, C.: A continuum method for modeling surface tension, J. Comput. Phys., 100, 335–354, https://doi.org/10.1016/00219991(92)90240Y, 1992.
Centler, F., Shao, H., De Biase, C., Park, C.H., Regnier, P., Kolditz, O., and Thullner, M.: GeoSysBRNSA flexible multidimensional reactive transport model for simulating biogeochemical subsurface processes, Comput. Geosci., 36, 397–405, https://doi.org/10.1016/j.cageo.2009.06.009, 2010.
Cirpka, O. A. and Valocchi, A. J.: Twodimensional concentration distribution for mixingcontrolled bioreactive transport in steady state, Adv. Water Resour., 30, 1668–1679, https://doi.org/10.1016/j.advwatres.2006.05.022, 2007.
Cirpka, O. A. and Valocchi, A. J.: Reply to comments on ”Twodimensional concentration distribution for mixingcontrolled bioreactive transport in steady state” by H. Shao et al., Adv. Water Resourc., 32, 298–301, 10.1016/j.advwatres.2008.10.018, 2009.
Danckwerts, P. V. and Lannus, A.: Gasliquid reactions, J. Electrochem. Soc., 117, 96–352, https://doi.org/10.1149/1.2407312, 1970.
Gharasoo, M., Centler, F., Regnier, P., Harms, H., and Thullner, M.: A reactive transport modeling approach to simulate biogeochemical processes in pore structures with porescale heterogeneities, Environ. Modell. Softw., 30, 102–114, 2012.
Golparvar, A., Kästner, M., and Thullner, M.: Porescale modeling of microbial activity: What we have and what we need, Vadose Zone J., 20, e20087, https://doi.org/10.1002/vzj2.20087, 2021.
Golparvar, A.: Movies, figshare, https://doi.org/10.1002/vzj2.20087, 2022.
Graham, E., Grandy, S., and Thelen, M.: Manure effects on soil organisms and soil quality, Emerging Issues in Animal Agriculture, Michigan State University Extension, 2014.
Graveleau, M., Soulaine, C., and Tchelepi, H. A.: PoreScale Simulation of Interphase Multicomponent Mass Transfer for Subsurface Flow, Transport Porous Med., 120, 287–308, https://doi.org/10.1007/s1124201709211, 2017.
Greenshields, C. J.: OpenFOAM user guide, OpenFOAM Foundation Ltd, version, 3, 47, 2015.
Haroun, Y., Legendre, D., and Raynal, L.: Volume of fluid method for interfacial reactive mass transfer: Application to stable liquid film, Chemi. Eng. Sci., 65, 2896–2909, https://doi.org/10.1016/j.ces.2010.01.012, 2010.
Hirt, C. W. and Nichols, B. D.: Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys., 39, 201–225, https://doi.org/10.1016/00219991(81)901455, 1981.
Issa, R. I.: Solution of the implicitly discretised fluid flow equations by operatorsplitting, J. Comput. Phys., 62, 40–65, 1986.
Kampschreur, M. J., Kleerebezem, R., Picioreanu, C., Bakken, L. R., Bergaust, L., de Vries, S., Jetten, M. S., and Van Loosdrecht, M.: Metabolic modeling of denitrification in Agrobacterium tumefaciens: a tool to study inhibiting and activating compounds for the denitrification pathway, Front. Microbiol., 3, 370, https://doi.org/10.3389/fmicb.2012.00370, 2012.
Kuzyakov, Y. and Blagodatskaya, E.: Microbial hotspots and hot moments in soil: Concept & review, Soil Biol. Biochem., 83, 184–199, https://doi.org/10.1016/j.soilbio.2015.01.025, 2015.
Li, X. Y., Huang, H., and Meakin, P.: A threedimensional level set simulation of coupled reactive transport and precipitation/dissolution, Int. J. Heat Mass Tran., 53, 2908–2923, https://doi.org/10.1016/j.ijheatmasstransfer.2010.01.044, 2010.
Meakin, P. and Tartakovsky, A. M.: Modeling and simulation of porescale multiphase fluid flow and reactive transport in fractured and porous media, Rev. Geophys., 47, 2008RG000263, https://doi.org/10.1029/2008RG000263, 2009.
Meile, C. and Scheibe, T. D.: Reactive Transport Modeling of Microbial Dynamics, Elements, 15, 111–116, https://doi.org/10.2138/gselements.15.2.111, 2019.
Nick, H. M., Raoof, A., Centler, F., Thullner, M., and Regnier, P.: Reactive dispersive contaminant transport in coastal aquifers: numerical simulation of a reactive Henry problem, J. Contam. Hydrol., 145, 90–104, https://doi.org/10.1016/j.jconhyd.2012.12.005, 2013.
Parkhurst, D. L. and Appelo, C. A. J.: User's guide to PHREEQC (Version 2): a computer program for speciation, batchreaction, onedimensional transport, and inverse geochemical calculations, Publisher: U.S. Geological Survey, Report 994259, https://doi.org/10.3133/wri994259, 1999.
Patankar, S. V. and Spalding, D. B.: A calculation procedure for heat, mass and momentum transfer in threedimensional parabolic flows, Int. J. Heat Mass Tran., 15, 1787–1806, https://doi.org/10.1016/00179310(72)900543, 1972.
Popinet, S.: An accurate adaptive solver for surfacetensiondriven interfacial flows, J. Comput. Phys., 228, 5838–5866, https://doi.org/10.1016/j.jcp.2009.04.042, 2009.
Raeini, A. Q., Blunt, M. J., and Bijeljic, B.: Modelling twophase flow in porous media at the pore scale using the volumeoffluid method, J. Comput. Phys., 231, 5653–5668, 2012.
Regnier, P., O'Kane, J. P., Steefel, C. I., and Vanderborght, J. P.: Modeling complex multicomponent reactivetransport systems: towards a simulation environment based on the concept of a Knowledge Base, Appl. Math. Model., 26, 913–927, https://doi.org/10.1016/S0307904X(02)000471, 2002.
Sander, R.: Compilation of Henry's law constants (version 4.0) for water as solvent, Atmos. Chem. Phys., 15, 4399–4981, https://doi.org/10.5194/acp1543992015, 2015.
Scardovelli, R., and Zaleski, S.: Direct numerical simulation of freesurface and interfacial flow, Annu. Rev. Fluid Mech., 31, 567–603, https://doi.org/10.1146/annurev.fluid.31.1.567, 1999.
Schlüter, S., Zawallich, J., Vogel, H.J., and Dörsch, P.: Physical constraints for respiration in microbial hotspots in soil and their importance for denitrification, Biogeosciences, 16, 3665–3678, https://doi.org/10.5194/bg1636652019, 2019.
Schmidt, H., Vetterlein, D., Köhne, J. M., and Eickhorst, T.: Negligible effect of Xray ìCT scanning on archaea and bacteria in an agricultural soil, Soil Biol. Biochem., 84, 21–27, 2015.
Steefel, C. and MacQuarrie, K.: Approaches to modeling of reactive transport in porous media, Rev. Miner. Geochem., 34, 85–129, 1996.
Steefel, C. I., Appelo, C. A. J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., Lagneau, V., Lichtner, P. C., Mayer, K. U., Meeussen, J. C. L., Molins, S., Moulton, D., Shao, H., Šimùnek, J., Yabusaki, S. B., and Yeh, G. T.: Title Reactive transport codes for subsurface environmental simulation, Comput. Geosci., 19, 445–478, 2015a.
Steefel, C. I., Appelo, C. A. J., Arora, B., Jacques, D., Kalbacher, T., Kolditz, O., Lagneau, V., Lichtner, P. C., Mayer, K. U., Meeussen, J. C. L., Molins, S., Moulton, D., Shao, H., Šimùnek, J., Spycher, N., Yabusaki, S. B., and Yeh, G. T.: Reactive transport codes for subsurface environmental simulation, Comput. Geosci., 19, 445–478, https://doi.org/10.1007/s105960149443x, 2015b.
Thullner, M., Van Cappellen, P., and Regnier, P.: Modeling the impact of microbial activity on redox dynamics in porous media, Geochim. Cosmochim. Ac., 69, 5005–5019, https://doi.org/10.1016/j.gca.2005.04.026, 2005.
Thullner, M., Regnier, P., and Van Cappellen, P.: Modeling Microbially Induced Carbon Degradation in RedoxStratified Subsurface Environments: Concepts and Open Questions, Geomicrobiol. J., 24, 139–155, https://doi.org/10.1080/01490450701459275, 2007.
Thullner, M.: Comparison of bioclogging effects in saturated porous media within one and twodimensional flow systems, Ecol. Eng., 36, 176–196, https://doi.org/10.1016/j.ecoleng.2008.12.037, 2010.
Thullner, M. and Regnier, P.: Microbial Controls on the Biogeochemical Dynamics in the Subsurface, Rev. Miner. Geochem., 85, 265–302, https://doi.org/10.2138/rmg.2019.85.9, 2019.
Tian, Z. W. and Wang, J. Y.: Lattice Boltzmann simulation of biofilm clogging and chemical oxygen demand removal in porous media, Aiche J., 65, e16661, https://doi.org/10.1002/aic.16661, 2019.
van Leer, B.: Towards the ultimate conservative difference scheme. V. A secondorder sequel to Godunov's method, J. Comput. Phys., 32, 101–136, https://doi.org/10.1016/00219991(79)901451, 1979.
Whitaker, S.: The method of volume averaging, Springer Science & Business Media, 2013.
White, A. F. and Brantley, S. L.: The effect of time on the weathering of silicate minerals: why do weathering rates differ in the laboratory and field?, Chem. Geol., 202, 479–506, https://doi.org/10.1016/j.chemgeo.2003.03.001, 2003.
Yan, Z. F., Liu, C. X., ToddBrown, K. E., Liu, Y. Y., BondLamberty, B., and Bailey, V. L.: Porescale investigation on the response of heterotrophic respiration to moisture conditions in heterogeneous soils, Biogeochemistry, 131, 121–134, https://doi.org/10.1007/s1053301602700, 2016.