On pore-scale modeling and simulation of reactive transport in 3D geometries

Pore-scale modeling and simulation of reactive flow in porous media has a range of diverse applications, and poses a number of research challenges. It is known that the morphology of a porous medium has significant influence on the local flow rate, which can have a substantial impact on the rate of chemical reactions. While there are a large number of papers and software tools dedicated to simulating either fluid flow in 3D computerized tomography (CT) images or reactive flow using pore-network models, little attention to date has been focused on the pore-scale simulation of sorptive transport in 3D CT images, which is the specific focus of this paper. Here we first present an algorithm for the simulation of such reactive flows directly on images, which is implemented in a sophisticated software package. We then use this software to present numerical results in two resolved geometries, illustrating the importance of pore-scale simulation and the flexibility of our software package.


Introduction
Understanding and controlling reactive flow in porous media is important for a number of environmental and industrial applications, including oil recovery, fluid filtration and purification, combustion and hydrology.Traditionally, the majority of theoretical and experimental research into transport within porous media has been carried out at macroscopic (Darcy) scale.Despite the progress in developing devices to perform experimental measurements at the pore-scale, experimental characterization of the pore-scale velocity, pressure and solute fields is still a challenging task.Due the influence of the porescale geometry on the processes of interest, direct numerical simulation (DNS) promises to be a very useful computational tool in a wide range of fields, and in combination with experimental studies, can be used to determine quantities of interest that are not experimentally quantifiable (Blunt et al., 2013).
Significant progress over the past 10 -15 years in the pore-scale simulation of single phase flow has resulted in the computation of permeability tensors for natural and technical porous media becoming a standard procedure.A number of academic as well as commercial software tools, capable of processing 3D CT images in addition to virtually generating porous media, are available, for example Aviso, GeoDict and Ingrain.Most of those software tools have the additional ability of simulating two-phase immiscible flow at the pore-scale directly on a computational domain obtained through the segmentation of 3D CT images, often using the lattice Boltzmann (LB), the level set or volume of fluid methods.In contrast, substantially less work on the DNS of reactive flow has been performed, and only a few software tools with this capability exist.A limited number of computational studies examining reactive transport where the reactions only occur within the fluid phase (and not at a surface) exist (Molins et al., 2012;Shen et al., 2011).In contrast, the literature and computational tools examining full 3D reactive flow where the reactions occur at the pore wall (surface reactions), is sparse.The majority of existing studies and available numerical simulation packages describing pore-scale sorptive transport are based on pore-network mathematical models (see, for example, Raoof et al. (2010); Varloteaux et al. (2013) and literature therein), where the geometry needs to be converted into an idealized series of connected pores and throats to represent the porous medium.However, during this process, information on the morphology of the underlying media can be lost (see, for example, Raoof et al. (2010) and Lichtner and Kang (2007)).In this paper we present an algorithm and a sophisticated software package, called Pore-Chem, which uses cell centered finite volume (FV) methods to numerically solve 3D solute transport with surface reactions at the pore-scale.In particular the software package has the ability to solve the systems of equations modeling colloidal reactive transport on a geometrical domain obtained directly through imaging techniques, such as computerized tomography, which allows for a very accurate spatial description of the computational domain.
In this paper we describe the transport of a generic solute through the Navier-Stokes (NS) system of equations coupled to a convection-diffusion (CD) equation.The CD equation is complemented by boundary conditions which describe various types of surface reactions comprising a Robin boundary condition for the solute coupled to an ordinary differential equation (ODE) describing the dynamics of the adsorption at the pore wall.To represent the computational domain, a voxelized geometry is considered, where each individual voxel is either solid or fluid.The solute transport model is chosen due to its applicability to a broad range of problems.In addition, our goal is to describe the transport and reaction of sub-micron particles, for which inertial effects of the individual particles are negligible.Discrete models, where each particle is modeled as an individual entity, are necessary when considering larger particles, for example with a radius greater than one micron, for which inertial effects become important.Several commercial software packages, for example GeoDict (GeoDict, 2012(GeoDict, -2014)), solve a range of discrete mathematical models describing colloid transport and adsorption.However, numerical simulation of these models is often significantly more computationally expensive than a continuous mathematical model and accounting for different reaction kinetics is usually not possible in such packages.
Reactive flow in porous media is intrinsically a multiscale problem.The goal of our developments is to support problems where scale separation is possible and in cases where it is not possible.The first case, where the separation of scales is viable, is usually the focus of asymptotic homogenization theory.In the second case, when scale separation is not possible, numerical upscaling methods like multiscale finite element methods are often applied.During the homogenization procedure, when applicable, certain assumptions are imposed, allowing for the derivation of macroscopic (Darcy scale) equations from the microscale formulation, with effective parameters, such as the permeability and the effective reaction rate, obtained through solution of a cell problem (Hornung, 1997).A number of studies have employed homogenization theory to derive a macroscopic description of sorptive reactive transport for particular parameter regimes.The homogenization of solute transport in porous media in the presence of surface reactions has been performed for both high Péclet numbers (convection dominated regime) (Allaire et al., 2010b,a;Allaire and Hutridurga, 2012), and when the Péclet number is of order one (Hornung, 1994;Kumar et al., 2013;van Duijn, 1991).In addition to being able to solve cell problems in a number of settings, our software has the ability of solving a much broader class of problems at the pore-scale, without being restricted by the assumptions required during homogenization.Furthermore, it provides the possibility to study various different types of surface reactions described by different kinetics.
The remainder of the paper is organized as follows.In Section 2 the mathematical model is presented and cast into dimensionless variables.The method of achieving a numerical solution to the system of equations is outlined in Section 3, and illustrative results using this method are presented in Section 4. Finally, conclusions are drawn in Section 5.

Mathematical model
We now detail the mathematical model, which describes the transport and reaction of a generic solute at a 2D interface within a 3D pore-scale resolved geometry, where we assume that the number of solute particles is sufficiently large that representation within a continuum framework is valid.
Let us denote the spatial domain of interest by Ω, an open subset of R 3 .We assume that we can decompose Ω into a solid domain, denoted by Ω s , and a fluid domain, denoted by Ω f , such that Ω = Ω s ∪ Ω f .Denoting the external boundary of our domain, being the closure of Ω, by ∂Ω, we partition this into an inlet, ∂Ω in , an outlet, ∂Ω out , and external walls, ∂Ω wall , so that ∂Ω = ∂Ω in ∪ ∂Ω out ∪ ∂Ω wall . (1) We note that, although we here consider only one inlet and one outlet, extension to consider multiple inlets and outlets is straightforward.Finally, we denote the interfacial boundary between the fluid and solid portions of the domain by Γ = Ω f ∩ Ω s .To allow for different types of reactive boundary conditions to be described, we decompose Γ into N ≥ 1 different boundary types as follows where Γ i has distinct properties.Making such a decomposition enables us to attribute different reaction rates or kinetic descriptions to different portions of the domain, allowing for different types of solid material.
In order to describe the flow of the water through the membrane, by appealing to the conservation of momentum, the incompressible NS equations are used: where v(x, t) and p(x, t) are the velocity and pressure of the fluid respectively, while µ ≥ 0 and ρ ≥ 0 are the viscosity and the density of the fluid which we assume are constants (Bear, 1988;Acheson, 1995).Suitable boundary conditions on ∂Ω are given by where V in is the inlet velocity with V in > 0, P out is the pressure at the outlet, and n is the outward pointing normal to the boundary ∂Ω.Although here we have used no-slip and no-flux flow conditions for the external walls, symmetry or periodic boundary conditions can also be imposed which may be more appropriate depending on the problem to be solved.Further boundary conditions are required to be specified on the reminder of the boundary to Ω f , being the fluid-solid interface.To allow for the slip of flow along the fluid-solid interface, and the inclusion of additional effects such as charged fluids or matrices, we use the Navier-Maxwell slip conditions, given by where β i is the slip length on x ∈ Γ i for i = 0, 1, . . .N − 1 measured per unit length, t is any unit tangent to the surface such that t • n = 0, and the superscript T denotes the transpose (Lauga et al., 2007).In the case that β i = 0 for some i, then the standard no-slip and no-flux boundary conditions for the flow are enforced on Γ i .We specify initial conditions through where v 0 and p 0 are known functions.We discuss the choice for these in Section 3. We denote the concentration of the solute within the fluid by c(x, t), measured in particle number per unit volume.Appealing to the conservation of mass and assuming no fluid-phase reactions occur, the spatio-temporal evolution of the solute concentration is given by where D ≥ 0 is the solute diffusion coefficient which we assume to be scalar and constant.We assume a known concentration of the solute at the inlet, and prescribe zero flux of the solute at the outlet and on the external walls as follows: where c in > 0 is assumed to be constant.

Models for surface reactions
We are required to specify boundary conditions for c(x, t) on x ∈ Γ i , to describe the surface reactions occurring here.In general, there are two stages of adsorption of a particle from the bulk solution to the solid surface.The first stage involves the diffusion of particles from the bulk solution to the subsurface and the second stage then involves the transfer of particles from the subsurface to the surface.After the adsorption of a molecule at the interface, there is a reorientation of the colloid molecules at the surface, which results in a change in the surface tension (Kralchevsky et al., 2008).Assuming that both the rate of diffusion of the particle from the bulk to the subsurface, and the rate of the transfer of the particles from the subsurface to the surface are important in determining the overall rate of reaction, we use a mixed kinetic-diffusion adsorption description, given by Here m(x, t) is the surface concentration of the particle under consideration (Kralchevsky et al., 2008), measured in units of number per unit surface area, which contrasts with c(x, t) being measured in units of number per unit volume.The function f i (c, m) describes the kinetics of the rate of change of the surface concentration on the i th reactive boundary for i = 0, . . .N − 1 (Danov et al., 2002).Equation ( 9) states that the change in the surface concentration is equal to the flux across the surface, where the movement from the bulk to the surface is termed adsorption, and movement from the surface to the bulk is termed desorption.If Γ i is nonreactive, for some i, then f i = 0, so the adsorbed concentration on this boundary type remains constant, and a no-flux boundary condition for the solute concentration is prescribed.For reactive boundaries, the choice of f i and its dependence on c and m is highly influential in correctly describing the reaction dynamics at the solid-fluid interface.A number of different isotherms exist for describing these dynamics, dependent on the solute attributes, the order of the reaction, and the interface type.
The simplest of these is the Henry isotherm, which assumes a linear relationship between the surface pressure and the number of adsorbed particles, and takes the form Here κ a,i ≥ 0 is the rate of adsorption, measured in unit length per unit time, and κ d,i ≥ 0 is the rate of desorption, measured per unit time, at reactive boundary type i for i = 0, . . .N − 1. Equation ( 10) states that the rate of adsorption is proportional to the concentration of particles in solution at the reactive surface.As such, the rate of adsorption predicted does not saturate at higher surface concentrations.However, physically we expect the rate of adsorption to decrease as the quantity of adsorbed particles increases and the available surface area for adsorption decreases.Even though the Henry isotherm predicts no limit to surface concentration and does not model any interaction between the particles, it has been used in a large number of analytical studies due to its linearity.The Langmuir adsorption isotherm was the first to be derived mathematically, and is suitable to describe the adsorption of a monolayer of localized non-ionic non-interacting molecules at a 2D solid interface, and a derivation from statistical physics may be found in (Baret, 1969).It is also frequently used to describe the adsorption of molecules at a solid-liquid interface, and is described by Here m ∞,i > 0 is the maximal possible adsorbed surface concentration, measured in number per unit area, at reactive boundary type i for i = 0, . . .N − 1.In comparison to the Henry isotherm, the Langmuir isotherm predicts a decrease in the rate of adsorption as the adsorbed concentration increases due to the reduction in available adsorption surface.The Henry isotherm, given in Equation ( 10), is a linearization of Equation ( 11), explaining why it produces an accurate representation only at low surface concentrations.
More complex descriptions of the reaction kinetics exist to describe non-localized adsorption and particles which interact.For example, the Frumkin isotherm describes localized adsorption where particle interaction is accounted through an additional parameter: where k is the Boltzmann constant and T is the temperature in Kelvin.The parameter β ≥ 0 describes how cooperative the reaction is, and is related to the interaction energy between the particles.In the case that β = 0, the Langmuir isotherm is recovered.The Frumkin isotherm is infrequently used in the mathematical modeling of colloidal transport, which may be due its nonlinearity and the difficulties in determining the interaction energy between the particles.We make the assumption that the adsorption or desorption of our solute does not alter the position of the reactive boundary, which in the case of small volumes of particles being adsorbed is sufficiently accurate.By the conservation of mass, such an assumption implies any adsorption or desorption on the surface is represented by a corresponding increase or decrease in the density of the solid material through time.In some cases, for example when the molecules are big or the number being adsorbed is large, interface evolution needs to be considered and may be done in a similar manner to Tartakovsky et al. (2008); Roubinet and Tartakovsky (2013) and Boso and Battiato (2013).This is particularly important in the application of rock dissolution and precipitation, where large geometrical changes are observed.
To close the system of equations, we impose the initial conditions where c 0 and m 0,i are known for i = 0, 1, . . .N − 1.
Our problem is, therefore, described by two systems of equations with one-way coupling; the incompressible NS equations, described by ( 3) -( 6) and the CD equation described by ( 7) -( 8), with reactive boundaries conditions (9), initial conditions (13), and a description of the reaction kinetics, for example Equation ( 10), ( 11) or ( 12).We now cast the equations into dimensionless variables, before detailing the methods used to obtain a numerical solution.

Nondimensionalization
Using a caret notation to distinguish the dimensionless variable from its dimensional equivalent, we let where L > 0 is a typical length scale of the computational domain and V in = V in .As our computational domain consists of voxels, the relationship between each voxel and its material property is conserved upon nondimensionalization, while the length, area and volume of each voxel scales with L, L 2 and L 3 respectively.Given this, we let Ω, Ωs and Ωf , with boundaries ∂ Ωin , ∂ Ωout , ∂ Ωwall and Γi for i = 0, 1, . . .N −1 represent the dimensionless versions of the equivalent dimensional domains and boundaries, where the voxel size is scaled accordingly.In dimensionless variables, we, therefore, have where are the global Reynolds and Péclet numbers respectively, being the ratio between the inertial and viscous forces and the ratio between advective and diffusive transport rates respectively.Boundary conditions are given by v = n, x ∈ ∂ Ωin , p = 0 and ∇v with βi = β i L .In the case of the Henry isotherm, we have whereas, nondimensionalization of the Langmuir and Frumkin isotherms yields respectively, for i = 0, 1, . . .N − 1, where β = 2βc in L kT and m∞,i = m ∞,i c in L .In ( 17) -( 19) the Damköhler numbers are given by Da a,i = κ a,i V in , and Da and describe, for each i, the ratio of the rate of reaction (either adsorptive or desorptive) to the rate of advective transport.The initial conditions are given by v(x, 0) = v0 (x), p(x, 0) = p0 (x), ĉ(x, 0 where v0 We now proceed to discuss the numerical methods used to obtain an approximate solution to our dimensionless system of equations given by ( 14), with boundary conditions given by ( 16) -( 18) and initial conditions specified through (21).

Numerical methods
The full system of equations cannot be solved using analytical techniques, and so numerical methods need to be employed to calculate an approximate solution.We employ FV methods to numerically solve our system of equations, motivated by their local mass conserving properties.Other methods, for example LB or finite difference methods may also be used for solving the flow problem.We note that our CD solver is completely compatible with LB methods (the compatibility of our solver with finite difference methods depends on the grid selection).Although the authors are not aware of a detailed comparison of the performance of FV and LB methods in solving the NS equations at the pore-scale, some incomplete internal studies indicate that LB methods can be advantageous for geometries with a very low porosity and a high tortuosity, while FV methods are favorable in other cases.
Due to the one-way coupling between our two systems of equation, the velocity and pressure solutions are at steady state.For the sake of generality, we consider the unsteady equations, and begin by solving the system of equations describing the fluid flow, namely (14a), (14b) along with (16a) -(16d), to obtain a steady state numerical solution where ∂ v ∂ t = 0 is satisfied.This is achieved using a Chorin fractional timestepping method and we refer the reader to Čiegis et al. (2007) and Lakdawala (2010) for full details and for further references, where the methods used are described.
Once the solution of v is obtained, we proceed to solve the system of equations describing the solute transport and reaction, (14c) along with the boundary conditions (16e) -( 18) and initial conditions (21), using a FV method with a cell-centered grid.For the sake of brevity, full details of the numerical method employed are not given, but we refer the reader to, for example Causon et al. (2011), for more detailed information on FV methods.Firstly, dimensionless time is uniformly partitioned into Q time points, denoted by t0 , t1 , . . .tQ−1 with tk = k ∆ t , where ∆ t the dimensionless timestep size.Then the spatial domain, Ω, is split into P non-overlapping cubic finite volumes, B l for l = 0, 1, . . .P 1 , which span the 3D computational domain, such that Ω = ∪ P −1 l=0 B l .Considering a single representative finite volume, B l , we denote its six faces by F l,j with center xj where the subscript j = e, w, n, s, t, b denotes the east, west, north, south, top and bottom faces respectively.Integration of (14c) over the control volume, B l , and time interval [ tk , tk+1 ], upon application of the divergence theorem, yields where ∂B l is the boundary of B l , so that ∂B l = j=e,w,n,s,t,b F l,j .Denoting the center of the finite volume by x c and using the approximations φ(x, t)dS = Âj φ(x j , t) and for some scalar function φ for j = n, s, e, w, t, b, where Âj is the area of the face F l,j , we have Denoting ĉj (τ ) = ĉ(x j , τ ) and vj (τ ) = v(x j , τ ) for j = n, s, e, w, t, b, c, by first order finite difference methods we have where ĉk j = ĉ(x j , tk ) for j = n, s, e, w, t, b, c.In the case that one of the faces of the control volume lies on a boundary, the appropriate boundary conditions must be discretized; for the inlet, outlet and solid (or symmetry) boundaries this is straightforward due to the Dirichlet and zero Neumann boundary condition imposed there via (16f) and (16g).Therefore, we omit the details for the discretization of the boundary conditions on the external boundary ∂ Ω.The appropriate discretization of the reactive boundary conditions, prescribed on the fluid-solid interface through (16e) and the corresponding description of the reaction kinetics, here either ( 17), ( 18) or ( 19), is slightly more involved and deserves a more detailed discussion.
In a fully implicit and coupled discretization the resulting discrete equations are nonlinear and the Newton method needs to be used (Kelley, 1995).In a broad class of practically interesting problems we have considered to date, we have not faced very strong coupling between the dissolved and adsorbed concentrations.Therefore, a fully implicit and coupled discretization was not required and we have found that an operator splitting approach, or just a Picard linearization, has worked well.In this approach, the dissolved particle concentration is computed at t = tk+1/2 , and then the is used to compute the deposited mass at the time t = tk+1 .Runge-Kutta methods, or other methods for numerically solving stiff ODEs, are also straightforward to implement, and may be the subject of future studies if required (Kelley, 1995).In order to illustrate the method used in Pore-Chem, we describe the discretization for the Langmuir isotherm, which is achieved as follows.Firstly ( 18) is substituted into (16e), which is then split into a Robin boundary condition and an ordinary differential equation: In both of these cases, by ( 18), fi = 0 and so no reactions occur at the spatiotemporal point under consideration, in which case, by (24a), we have ∇ĉ Here B is a constant of integration which may be evaluated at t = tk to give Upon substitution into (25), we have for x ∈ Γi , t > 0, where we have approximated ĉ(x) by ĉ(x, tk ).This is done to prevent nonlinear terms in unknown variables from appearing in the discretized version of (16e).Discretization of the other two isotherms is implemented in a similar manner.Consequently, we may approximate the Robin boundary condition, (16e), on the reactive face F l,j ∈ Γi for j = e, w, n, s, t, b and i = 0, 1, . . .N − 1 using finite difference methods fully implicitly as follows: where n = ±1 is the direction of the outward pointing normal.Using (28), the appropriate terms are assembled, along with (23) minus the relevant diffusive flux term, for the finite volume on which the reactive surface lies into a matrix A k+1 and vector g k+1 , where A k+1 ĉk+1 = g k+1 and ĉk+1 is a vector of the dimensionless solutions ĉ(x, tk+1 ) at the discretized points of the computational domain.Once all the terms for all the finite volumes within the domain have been assembled into A k+1 and g k+1 , the equation A k+1 ĉk+1 = g k+1 is solved using a biconjugate gradient stabilized method to give the updated fluid concentration, ĉk+1 , at each discrete spatial point in Ωf , while the updated adsorbed concentration, mk+1 , is given by (26).Time is then updated, the next timestep considered, and we proceed in the usual manner until the final time point is reached.Numerical simulation of the equations is performed in Pore-Chem and a schematic, outlining the numerical algorithm used, is given in Figure 1.In the following section, the dimensionless equations are solved, but we present results in dimensional quantities.

Results
We present illustrative results, using the numerical method outlined in Section 3, on two separate computational domains, a real geometry and a virtually generated geometry, for two applications where surface reactions are important.The first set of simulations are performed on a portion of palatine sandstone, obtained using micro computerized tomography (µ-CT) by Frieder Enzmann at the Johannes Gutenberg University of Mainz (Becker et al., 2011).Surface reactions are highly important in many fields of the Earth sciences, including calcite growth, oxidation-reduction reactions, formation of biofilms, to name only a few (Steefel et al., 2005).The second set of results is performed on a computational domain virtually constructed to be representative of a commercially available microfiltration functionalized membrane.The use of functionalized membranes is a promising method for removing contaminants from water, and involves treating the pore walls of the membrane so that they adsorb certain microorganisms or drugs (Ulbricht, 2006).Such membranes have pore sizes on the sub-micron scale and the resolution provided by µ-CT imaging techniques is not high enough to give representative images, motivating the use of a virtually generated geometry.The two computational domains under consideration are shown in Figure 2.
For the numerical simulations presented, a dead-end setup is used, with a schematic illustrating the domain and boundary conditions shown in Figure 3. Layers of pure water voxels are added at the inlet and at the outlet, as shown in Figure 2 and illustrated in Figure 3, to allow free flow to develop.This results in a total computational domain size of 200 × 200 × 300 voxels for the sandstone geometry, and 100 × 100 × 120 voxels for the membrane geometry.For the results presented one type of reactive solid-fluid interface is considered, so that N = 1, and consequently we now drop the i subscript notation relating to the type of reactive boundary.

Fluid flow
As described in Section 3, we are first required to solve for the fluid flow.We assume that the fluid generates no slip as it passes over the membrane, and we set the slip length, β, to be zero so that, by (5), there is zero normal and tangential velocity at the fluid-solid interface, Γ.We set inflow velocities to be typical for the application under consideration, and use V in = 1 mm/s for the membrane geometry and V in = 1.5×10 −4 mm/s for the rock geometry.The parameters chosen for the inlet velocities, and the fluid density and viscosity, yield small Reynolds numbers; Re = 4.2 × 10 −7 in the case of the rock geometry and Re = 7.83 × 10 −3 for the membrane geometry.Therefore, the fluid flow in both computational domains is in a Stokes regime.Solving (14a) and (14b), along with the boundary conditions (16a) -(16d), numerically until steady-state is achieved, yields the solutions as reproduced in Figure 4. Due to our assumption that the maximal possible number of adsorbed particles is sufficiently small to ignore the effects of geometry modification, these remain constant through time.Examination of Figure 4 reveals the dependence of the local velocity field on the morphology of the computational domains.

Contaminant transport
We now turn our attention to solving the reactive flow.Illustrative parameters are used, and we employ the Henry isotherm, given by Equation ( 10), for the rock geometry, and the Langmuir isotherm, given by Equation ( 11), for the membrane geometry.We choose c in arbitrarily to be 2×10 4 number/mm 3 and we set the initial concentration of fluid contaminant to be equal to the inflow boundary condition, so that ĉ0 (x) = 1 for x ∈ Ωf , while the quantity of adsorbed contaminant is initially assumed to be zero, m0 (x) = 0. Furthermore, for the membrane geometry, we set the dimensionless maximal surface concentration of adsorbed contaminant, m∞ , to be 10 −4 , which equates to a dimensional value of m ∞ = 0.014 number/mm 2 .Due to the form of the equations, the choice of c in does not influence the dimensionless system of equations describing the transport and reaction of the contaminant given by (14c) along with the boundary conditions (16e) -( 19), except for through the parameter m∞ .Therefore, In the case that the dimensional system of equations is being solved, the appropriate dimensionless quantities of interest are replaced by the dimensional equivalent.increasing or decreasing c in for fixed m∞ purely scales the dimensional concentration (both fluid and adsorbed) by a constant factor.Figure 5 illustrates the concentration of the contaminant on the solid surface over time for the rock geometry, where we use Da a = 0.1, Da d = 0.001 and Pe = 2.0 × 10 −5 .Due to the slow rate of reaction and small Péclet number, the quantity of solute adsorbed over this time period is not large enough to significantly alter the fluid concentration over the time period under consideration.Furthermore, as the rate of reaction is much smaller than the mass transport rate, the adsorbed concentration remains spatially homogeneous.
In contrast, with the parameters Pe = 10 and Da a = Da d = 10 for the membrane geometry, we see significant spatial heterogeneity in both the dissolved and adsorbed concentrations, as illustrated in Figure 6.As time progresses, the fast rates of reaction result in a depletion of the dissolved concentration at the pore wall and a dependence of both the dissolved and adsorbed concentrations on the local membrane morphology.

Conclusions
We have presented an algorithm for solving solute transport at the pore-scale within a resolved porous medium, with reversible surface adsorption at the pore wall.A pore-scale description of reactive transport, as opposed to a description at the Darcy scale, allows for a very accurate representation of the processes of interest.The system of equations comprise the NS equations and a CD equation, with Robin boundary conditions coupled to an ODE accounting for the surface reactions.Assuming that each particle is sufficiently small in size not to alter the flow of the fluid within the computational domain, and that its reaction at the wall does not significantly alter the pore-scale geometry, there is a oneway coupling between the NS and CD systems of equations.Although, for simplicity, we consider one species of solute and examine only surface reactions, extension to several different species of solute with both volumetric and surface reactions is straight forward and implemented within our software package Pore-Chem.The algorithm presented employs a FV method, and in this paper we have particularly focused on the discretization method used to solve the reactive boundary conditions for the adsorption and desorption at the interface.Illustrative numerical results, using our software package Pore-Chem, are presented on two separate geometries.The first of these is a 3D µ-CT image of a piece of Palatine Sandstone rock, while the second geometry is virtually generated within GeoDict (GeoDict, 2012(GeoDict, -2014) ) to be representative of a commercially available functionalized membrane.The results demonstrate the potential of such a numerical package, with the ability to solve reactive transport directly on images and on virtually generated geometries, in further progressing the understanding of the interplay between the transport and reaction rates at the pore-scale.In a future publication, currently in preparation, we will investigate the influence of the computational domain morphology on the reaction dynamics, and investigate the effect of different parameter regimes on the numerical results and quantities of interest.
The advantages of using a pore-scale description are multiple.Firstly, it allows us to simulate reactive transport over a range of different parameter regimes, and in particular, outside the applicability region of the equivalent upscaled model.In highly disordered media the Péclet number can significantly vary locally, which poses problems for asymptotic upscaling methods, but not for a pore-scale description.Secondly, different kinetic models for the reactions can be used without the need to re-perform the upscaling procedure.In the future we plan to extend the algorithm in order to solve coupled multiscale problems using the heterogeneous multiscale method, in a similar manner to Battiato et al. (2011) and Iliev et al. (2014).Such a development will enable problems at larger spatial scales to be considered, which could aid further research into the influence of pore-scale processes in a number of highly interesting and important research applications.

Figure 1 :
Figure 1: Flow chart to illustrate how numerical computation of the transport with surface reactions is implemented within Pore-Chem, with main features indicated only.The light red ovals indicate input files, while the yellow ovals indicate output files.The green and the blue rectangular boxes indicate steps involved for the flow and efficiency solvers respectively, while the dark red boxes indicate steps involved in both.The white diamonds indicate decision making steps.In the case that the dimensional system of equations is being solved, the appropriate dimensionless quantities of interest are replaced by the dimensional equivalent.

Figure 2 :
Figure 2: The two computational domains under consideration, plotted in 3D on the left-hand side of the figure, and in 2D through a representative cross section on the right-hand side of the figure.Figure (a) shows the palatine sandstone geometry obtained through µ-CT, with a voxel size of 1.4×10 −6 m, as used in Becker et al. (2011).Figure (b) shows a virtually generated geometry which aims to reproduce the morphology of a commercially available microfiltration membrane, with a voxel length of 7.03 × 10 −8 m (3 s.f.)This geometry was created using the software GeoDict (GeoDict, 2012-2014), and a comparison to experimentally evaluated quantities is made in Di Nicolò et al. (2015).
Figure 2: The two computational domains under consideration, plotted in 3D on the left-hand side of the figure, and in 2D through a representative cross section on the right-hand side of the figure.Figure (a) shows the palatine sandstone geometry obtained through µ-CT, with a voxel size of 1.4×10 −6 m, as used in Becker et al. (2011).Figure (b) shows a virtually generated geometry which aims to reproduce the morphology of a commercially available microfiltration membrane, with a voxel length of 7.03 × 10 −8 m (3 s.f.)This geometry was created using the software GeoDict (GeoDict, 2012-2014), and a comparison to experimentally evaluated quantities is made in Di Nicolò et al. (2015).
Figure 2: The two computational domains under consideration, plotted in 3D on the left-hand side of the figure, and in 2D through a representative cross section on the right-hand side of the figure.Figure (a) shows the palatine sandstone geometry obtained through µ-CT, with a voxel size of 1.4×10 −6 m, as used in Becker et al. (2011).Figure (b) shows a virtually generated geometry which aims to reproduce the morphology of a commercially available microfiltration membrane, with a voxel length of 7.03 × 10 −8 m (3 s.f.)This geometry was created using the software GeoDict (GeoDict, 2012-2014), and a comparison to experimentally evaluated quantities is made in Di Nicolò et al. (2015).

Figure 3 :
Figure 3: Schematic to illustrate the computational domain.The solid microfiltration membrane is show in grey while the water is shown in white.Voxels are added to the top and bottom of the domain, at the inlet and outlet, to enable free-flow to develop.Modified, with permission, from Di Nicolò et al. (2015).

Figure 4 :Figure 5 :
Figure 4: Velocity magnitude, (a) and (c), and pressure, (b) and (d), fields, measured in mm/s and kPa respectively for both the computational domains considered.Due to our assumption that the adsorption of the solute does not alter the geometry, these remain constant throughout the experimental time frame.The black boxes shows the outline for the entire computational domain, where we exclude the additional voxels embedded at the inlet and the outlet for 3D plotting purposes.
Figure 6: Numerical results at 5 × 10 −5 s intervals for the fluid concentration, c, and the adsorbed concentration, m, of contaminant in the membrane geometry.
• n = 0 and a zero Neumann boundary condition, which is straightforward to implement.Otherwise, if Da I > 0, assuming that ĉ(x, t) is constant over the time period in question, namely t ∈ [ tk , tk+1 ], and equal to ĉ(x) at each spatial point, (24b) may be integrated to give m