Lattice Boltzmann-Based Approaches for Pore-Scale Reactive Transport

Important geoscience and environmental applications such as geologic carbon storage, environmental remediation, and unconventional oil and gas recovery are best understood in the context of reactive flow and multicomponent transport in the subsurface environment. The coupling of chemical and microbiological reactions with hydrological and mechanical processes can lead to complex behaviors across an enormous range of spatial and temporal scales. These coupled responses are also strongly influenced by the heterogeneity and anisotropy of the geologic formations. Reactive transport processes can change the pore morphology at the pore scale, thereby leading to nonlinear interactions with advective and diffusive transport, which can strongly influence larger-scale properties such as permeability and dispersion. Therefore, one of the greatest research challenges is to improve our ability to predict these processes across scales (DOE 2007). The development of pore-scale experimental and modeling methods to study reactive processes involving mineral precipitation and dissolution, and biofilm dynamics allows more fundamental investigation of physical behavior so that more accurate and robust upscaled constitutive models can be developed for the continuum scale.

A pore-scale model provides fundamental mechanistic explanations of how biogeochemical processes and pore-scale interfacial reactions alter flow paths by pore plugging (and dissolving) under different geochemical compositions and pore configurations. For example, dissolved CO2 during geological CO2 storage may react with minerals in fractured rocks, confined aquifers, or faults, resulting in cementation (and/or dissolution) and altering hydrodynamics of reactive flow. This can be observed in a natural analogue where primary porosity in sandstone is cemented by carbonate precipitates, affecting dissolved CO2 flow paths at the Little Garde Wash Fault, Utah (e.g., Fig. 1a–b). Several other examples demonstrating macroscopic characteristics of calcium carbonate (CaCO3) precipitation in Figure 1 include an elongated concretion along the groundwater flow direction, CaCO3 precipitation along the vertical pathway sealed …


INTRODUCTION
Important geoscience and environmental applications such as geologic carbon storage, environmental remediation, and unconventional oil and gas recovery are best understood in the context of reactive fl ow and multicomponent transport in the subsurface environment.The coupling of chemical and microbiological reactions with hydrological and mechanical processes can lead to complex behaviors across an enormous range of spatial and temporal scales.These coupled responses are also strongly infl uenced by the heterogeneity and anisotropy of the geologic formations.Reactive transport processes can change the pore morphology at the pore scale, thereby leading to nonlinear interactions with advective and diffusive transport, which can strongly infl uence larger-scale properties such as permeability and dispersion.Therefore, one of the greatest research challenges is to improve our ability to predict these processes across scales (DOE 2007).The development of pore-scale experimental and modeling methods to study reactive processes involving mineral precipitation and dissolution, and biofi lm dynamics allows more fundamental investigation of physical behavior so that more accurate and robust upscaled constitutive models can be developed for the continuum scale.
A pore-scale model provides fundamental mechanistic explanations of how biogeochemical processes and pore-scale interfacial reactions alter fl ow paths by pore plugging (and dissolving) under different geochemical compositions and pore confi gurations.For example, dissolved CO 2 during geological CO 2 storage may react with minerals in fractured rocks, confi ned aquifers, or faults, resulting in cementation (and/or dissolution) and altering hydrodynamics of reactive fl ow.This can be observed in a natural analogue where primary porosity in sandstone is cemented by carbonate precipitates, affecting dissolved CO 2 fl ow paths at the Little Garde Wash Fault, Utah (e.g., Fig. 1a-b).Several other examples demonstrating macroscopic characteristics of calcium carbonate (CaCO 3 ) precipitation in Figure 1 include an elongated concretion along the groundwater fl ow direction, CaCO 3 precipitation along the vertical pathway sealed by a low-permeability layer, and calcite-cemented hanging-wall damage/mixed zone impacting permeability and fl ow paths.Continuum-scale reactive transport modeling has been studied for mineral precipitation and its impact on permeability (e.g., Steefel andLasaga 1990, 1994).The numerical studies suggest that diffusion-controlled precipitate can armor fractures and limit the amount of mixing and subsequent reaction (Steefel andLichtner 1994, 1998).Recently, mineral precipitation and its impact on fl ow alteration at the pore scale was evaluated under representative fl ow and mixing conditions in porous media (Tartakovsky et al. 2008;Zhang et al. 2010a;Yoon et al. 2012;Boyd et al. 2014).In these works, experimental and numerical pore-scale reactive transport investigations demonstrated that both thermodynamic and kinetic constraints affect precipitation rates, the distribution of mineral polymorphs, and the corresponding extent of porosity occlusion.Elemental analysis also suggested that crystal growth rates are affected by solution chemistry and location of reaction zones.In addition, other experimental results revealed the important effect of mineral formation on the estimation of diffusion coeffi cient (Navarre-Sitchler et al. 2009), overall reactive surface area (Navarre-Sitchler and Brantley 2007), and reactive surface area within connected pores (Landrot et al. 2012).Hence, accuracy of large-scale rate estimates depends on the type of reaction, pore geometry, reactive surface characteristics, reaction kinetics, macroscopic concentration gradient, and interaction between reaction and fl ow paths.As suggested in Battiato et al. (2011) and Steefel et al. (2013), Darcy-scale simulation with effective parameters may not account for nonlinear reactive transport processes associated with localized precipitation and dissolution of mineral phases.
Another notable example of porosity and permeability alteration due to reactive transport is microbial growth and the associated biogeochemical reaction products, such as biominerals.Microbial growth including extracellular polymeric substances (EPS) can signifi cantly reduce porosity and permeability by pore plugging (Davery et al. 1998;Zhang et al. 2010b;Kirk et al. 2012;Yoon et al. 2014).Additionally, biomineralization induced by urea hydrolysis and denitrifi cation processes (Nemati and Voordouw 2003;Dejong et al. 2006;Martin et al. 2013;Zhu et al. 2013) altered porosity and fl ow patterns signifi cantly.On the other hand, organic acids (e.g.acetic acid, lactic acid, etc.) that are produced during microbial growth can dissolve rocks, resulting in increasing in porosity and permeability (Adkins et al. 1992).It is generally believed that most subsurface microorganisms are attached onto the surface of the solid phase as a biofi lm, a complex mixture of multiple active microbial species, inert biomass, and extra-cellular polymeric substances (Rittmann 1993;Rittmann and McCarty 2001;IWA Task Group 2006;Vilcáez et al. 2013).Although volume-averaging theory can be used to derive continuum-scale models, idealized geometrical assumptions about the porous media and the biofi lm are necessary (Golfi er et al. 2009).Instead, pore-scale numerical models capable of handling complex dynamic solid-liquid boundaries have been developed to couple fl ow, solute transport, and biofi lm growth at the pore scale.Since biofi lm dynamics at the pore scale can affect macroscopic properties such as porosity, permeability, and dispersivity (Taylor and Jaffé 1990;Vandevivere and Baveye 1992;Seifert and Engesgaard 2007;Seifert and Engesgaard 2012;El Mountassir et al. 2014;Yoon et al. 2014), the fundamentals of biogeochemical processes at the pore scale must to be understood if we are to address practical fi eld-scale problems.
Recent advances in multiscale imaging techniques for the analysis of complex pore structures and mineralogies have revolutionized our ability to characterize geomaterials quantitatively at the pore scale (e.g., Gualda et al. 2010).By combining a suite of imaging techniques, it is now largely possible to understand how surface chemistry and pore topology impact reactive transport processes affecting the change of porosity and permeability (e.g., Landrot et al. 2012;Yoon et al. 2012).Similar recent advances in computational power and methods have led to the development of a number of pore-scale models for fl ow and reactive transport (see reviews in Meakin and Tartakovsky 2009;Steefel et al. 2013).To summarize, these models include lattice Boltzmann (LB) (Kang et al. 2003(Kang et al. , 2006(Kang et al. , 2010)), smooth particle hydrodynamics (SPH) (Tartakovsky et al. 2007(Tartakovsky et al. , 2008(Tartakovsky et al. , 2009)), direct numerical simulation (DNS) (Flukiger and Bernard 2009;Molins et al. 2012Molins et al. , 2014)), hybrid LB-DNS (Yoon et al. 2012) and pore network models (Li et al. 2006;Kim et al. 2011).In this article we review approaches based on the lattice Boltzmann method (LBM) for pore-scale reactive transport models relevant to the impact of biogeochemical processes on the dynamic permeability evolution.
The LBM, which was developed as an extension to lattice gas models for fl uid fl ows (Qian and Orszag 1995;Stockman et al. 1997;Chen and Doolen 1998), is well suited for studying hydrodynamics in complex geometries of porous and fractured media with a variety of reactive surface boundaries.Unlike the conventional computational fl uid dynamics at the continuum scale, the LBM employs a mesoscopic equation from a kinetic theory to determine macroscopic fl uid dynamics (Succi 2001).A basic introduction to the LBM for geoscientists and beginners including example (pseudo-) codes can be found in Sukop and Thorne (2007), Mohamad (2011), and open-source domains (Lattice Boltzmann methods 2015).LBM considers fl ow as a collective behavior of pseudo-particles residing at a mesoscopic level whose evolution is described by a discrete Boltzmann equation.Specifi cally, the collision term in the Boltzmann equation is approximated by the simple Bhathagar-Gross-Krook (BGK) method (Bhatnagar et al. 1954) or single relaxation time (SRT) model.The lattice Boltzmann Equation (LBE) can be obtained from the continuum Boltzmann equation with a BGK collision term by proper discretization in time, space, and momentum (Chen and Doolen 1998).The Navier-Stokes equation can be recovered in an incompressible limit using the Chapmann-Enskog multi-scale expansion technique (Wolf-Gladrow 2000).LBM can easily deal with complex boundaries using a simple bounce-back or modifi ed bounce-back scheme (He et al. 1997) in the fl uid particle distributions, and can be applied to multi-phase fl ow with phase transition (Chen and Doolen 1998) and particulate suspension fl ows (Ladd 1994;Joshi and Sun 2009).Since the LBE is an explicit time-stepping equation for the distribution function, it can be easily programmed on a computer.Based on the local nature of the kinetic reaction terms or processes, LBM is also suited for parallel computing.Recently a comparison of representative LB models for benchmark problems demonstrated a strong parallel scalability (Groen et al. 2013).
A number of studies have presented reactive transport models using LBM, focusing on crystallization processes (Miller et al. 2001;Kang et al. 2004;Lu et al. 2009), dissolution (Kelemen et al. 1995;Kang et al. 2002;Huber et al. 2008;Parmigiani et al. 2011;Chen et al. 2014b), and biofi lm dynamics (Sullivan et al. 2005(Sullivan et al. , 2006;;Pintelon et al. 2009).In particular, Kang and co-workers have developed multicomponent reactive transport models assuming homogeneous reactions are instantaneous, and all heterogeneous reactions at the mineralwater interface are kinetically controlled (Kang et al. 2002(Kang et al. , 2003(Kang et al. , 2006(Kang et al. , 2007)).LBM now has expanded capabilities for (reactive) transport processes by coupling with continuumbased or DNS methods (e.g., fi nite-volume method (FVM), fi nite-element method (FEM), fi nite-difference method (FDM)).In this article, the lattice Boltzmann-based approaches are referred to as the LBM coupled with other techniques such as LB for fl ow and FVM for reactive transport.The main objective of this article is to present the current state-of-the-art in using LB-based approaches to simulate coupled fl ow and transport affecting pore structure change and fl ow feedback.In the following subsections we briefl y introduce reactive transport modeling of biogeochemical processes.LBM for reactive transport models will be presented in detail.We also describe recent advances in LB applications to multiphysics systems useful for multiphase fl ow and electrical transport, and coupled LB-DNS approaches for reactive transport models.Representative LB-based simulations will be reviewed to demonstrate its usefulness in pore-scale reactive transport processes and future research directions will be addressed.

REACTIVE TRANSPORT MODELING OF BIOGEOCHEMICAL PROCESSES AT THE PORE SCALE
Multicomponent reactive transport models have been well outlined for a basic continuum theory for reactive transport (Lichtner 1985(Lichtner , 1996) ) and biofi lm modeling in water and wastewater treatment (Rittmann and McCarty 2001;IWA Task Group 2006).Recently, pore-scale reactive transport modeling has been actively developed and applied, which were fundamentally overviewed for multiphase fl ow and reactive transport models (Meakin and Tartakovsky 2009) and pore-scale processes during geologic carbon storage (Steefel et al. 2013).In this section, the mathematical equations for reactive transport modeling at the pore scale will be briefl y described for the rest of article and key biogeochemical processes relevant to pore structure dynamics will be introduced.The primary aspects of a reactive transport model include (1) governing equations for fl uid fl ow, (2) transport equations for multiple species, (3) (bio)geochemical reaction models describing homogeneous reactions (reactions within a single phase) and heterogeneous reactions at the interface of different phases, and (4) the motion of the heterogeneous interface due to the reactions.

Fluid fl ow
The fl uid fl ow in porous media is infl uenced by inertial, viscous, and body forces.In reservoirs and natural porous media systems where the compressibility of the fl uid is assumed small, the incompressible Navier-Stokes equations for the fl uid fl ow can be described by the conservation of momentum and mass: where u is the fl uid velocity vector,  is the fl uid density, p is the pressure,  is the dynamic viscosity of the fl uid, and F is the body-force density including all the effective external body forces.In Equation (1) the left-hand side represents the inertial terms and the right-hand side is written in terms of the full rate-of-deformation (stress) tensor owing to, for example, concentration-dependent viscosity.This dependence couples the fl ow system to the reactive transport system where the time-dependent part of the acceleration term is retained because of the concentration dependence of the viscosity (e.g., Davidson et al. 2012).In this article we will treat the concentration-(or reaction-) dependent viscosity as a special case.
If the viscosity is constant or independent of the species concentrations, the fl ow problem is uncoupled from the species transport.At low Reynolds numbers Re= |u|l / «1 where l is a characteristic length, the inertial terms can be neglected and under the steady conditions, the Stokes equation can be used to obtain the velocity fi eld at the pore scale: The quasi-static formulation for the fl ow implies that the geometry of the solid phase changes slowly compared to the time scale for reactive transport unless the reactions will alter the geometry of the solid phase.We will separately treat the effect of biogeochemical reactions on the change of solid (including biofi lm) phase later.

Multicomponent reactive transport
In the species transport equations at the pore scale, the detailed inclusion of fl uid fl ow at the pore scale allows us to model the diffusion by a simple Fickian form in contrast to a dispersion term at the continuum Darcy scale.For the i-th (ion) species, the reactive transport equation at the pore scale for an incompressible fl ow fi eld has a general form (Lichtner 1996) where C i is the concentration of the i-th species, J i is the species fl ux, R i is the reaction rates at which the i-th species is produced, consumed, and decayed by (bio)chemical reactions.The total fl ux J i consists of advection, diffusion, and electrochemical migration terms and can be described using the Nernst-Planck equation (Levich 1962 where D i is the diffusion coeffi cient of the i-th species, z i is the charge of the i-th species, F is the Faraday constant, R is the ideal gas constant, T is the temperature,  is the electrical potential.The multispecies diffusive fl uxes including the electrochemical migration can be important in dilute solutions (Van Cappellen and Gaillard 1996;Li et al. 2007;Lichtner and Kang 2007;Liu et al. 2011;Molins et al. 2012;Zhang and Klapper 2014).In this article, we will briefl y introduce the lattice Poisson-Boltzmann method to account for the electroosmotic fl ows as a special case.For the sake of simplicity, it is assumed that all ion species have the same diffusion coeffi cient to ensure local charge balance (Lichtner 1985) and avoid complicated charge-separation and charge-induced electrostatic interaction.The pore-scale advection-diffusion-reaction equation is where N R is the number of possible homogenous and heterogeneous reactions, v ir is the stoichiometric coeffi cient, and R r is the reaction rate.
To solve for a chemical system with many species, it is practically accepted to assume local partial equilibrium where both local equilibrium and kinetic reactions take place simultaneously.In this article, the homogeneous reactions are assumed to be in instantaneous equilibrium and all reactions at the heterogeneous interfaces are kinetically controlled.This approach is comprehensively discussed by Lichtner (1985Lichtner ( , 1996) ) and the reactive transport equations in terms of the total concentrations ( j ) are given by where N m is the number of kinetically controlled reactions, v jk is the stoichiometric coeffi cient in the reaction k with primary species j, R k is the reaction rate for the k-th mineral reaction at the mineral-fl uid interface, R bio is the reaction rate due to biologically mediated reactions (biomass growth and biological production), C j is the concentration of primary species j, N eq is the number of secondary species (subject to equilibrium reaction), and C i is the concentration of secondary species i.The C i is given through the mass action law as functions of the primary species concentrations: where  i is the activity coeffi cient and K i is the equilibrium constant.

Mineral precipitation and dissolution
The rate laws for mineral precipitation and dissolution at the pore scale can be described by transition state theory (TST) (Lasaga 1981) where the rate law is a function of the saturation state of the solution.The reaction rate can be expressed for all minerals (Lasaga 1998) with a general form by eq 1 , where a m is the specifi c surface area of the m-th mineral, k m is the reaction rate constant, Q s is the ion activity product, K eq is the equilibrium constant of the reaction, and n is the order of the reaction.The limitations of the applicability of the TST rate expressions were recently reviewed by Steefel et al. (2013) to address the issues on the reversibility implied in the TST formulation, the bulk TST approaches compared to atomistic reaction processes, the lack of time-dependent mineral surface reactivity, and the limitations on transport-controlled reactions.Nonetheless, recent pore-scale modeling results (Molins et al. 2012;Yoon et al. 2012;Hiorth et al. 2013) show that the TST approaches in a form of Equation ( 10) can match experimental results relatively well if a fi ne spatial resolution (at the order of micron scale) is used together with well-characterized reaction surfaces.
As stated earlier, all mineral precipitation and dissolution will take place at the mineralfl uid interface.In addition, the rate of the change of the interface due to Equation ( 10) is much slower than the scale of fl uid velocities, and a quasi-steady approximation of the fl uid fl ow can be used.Assuming no-slip boundary condition at the mineral surface, heterogeneous reactions at mineral interfaces can be represented through the internal boundary conditions (Kang et al. 2002;Li et al. 2008;Meakin and Tartakovsky 2009).The boundary condition for the total concentration is given by (Kang et al. 2006) 1 eq 1 , where n is the direction normal to the interface pointing toward the fl uid phase.Recently, Huber et al. (2014) where V m , m V , and m a are the volumetric fraction, molar volume, and specifi c surface area of the m-th mineral, respectively.Solute diffusion in the solid phase is neglected.The volume is updated at each time step explicitly according to the equation where dt is the time interval for the mineral update.

Biofi lm dynamics
A great deal of understanding biofi lm dynamics has been gained from study of water and wastewater treatment processes where the bacteria are attached to a solid media forming an immobile biofi lm.Modeling biofi lms in engineered reactors is well described by Rittmann and McCarty (2001) and IWA Task Group (2006).Following these works, the biofi lm in porous media is typically assumed to be impermeable to water fl ow, and therefore the total pore space can be partitioned into a fl uid and biofi lm phase.Dissolved substrates and solutes are allowed to diffuse within the biofi lms where the diffusion coeffi cient within the biofi lm phase (e.g., Eqn.6 or 7) is generally taken as roughly 80% of its value in water (Rittmann and McCarty 2001).In addition, several works have allowed the biofi lm to be slightly permeable to water fl ow to improve model predictions (Thullner and Baveye 2008;Pintelon et al. 2012;Deng et al. 2013).Assuming all the bacteria are attached to the solid as a biofi lm, a local mass balance equation applies, similar to the case of mineral precipitation and dissolution growth decay where M represents the mass of bacteria per unit volume of fl uid and r growth and r decay represent the bacteria growth and decay rates, respectively.
The Monod equation is widely used to link bacteria growth to consumption of dissolved substrates where Y is the biomass yield coeffi cient (mass of bacteria produced per mass of substrate consumed), C is the dissolved substrate concentration (mass per volume fl uid), k max is the maximum specifi c utilization rate (mass of substrate per mass of bacteria-time), and K c is the half-maximum-rate concentration.For the sake of simplicity it is assumed that there is a single limiting substrate for bacteria growth; however, several studies have used dual Monod kinetics to investigate situations where two substrates (e.g. an electron donor and acceptor) are limiting (Knutson et al. 2005;Tang et al. 2013).
The growth rate is directly coupled to substrate consumption, so formally r growth = -R bio in Equation ( 7).As biomass grows, it can spread out and occupy an increasing volume of the pore space.This biofi lm growth is very similar to growth of solid phase due to net precipitation, which will be discussed in the LB section.Growing bacteria requires energy for basic cell maintenance, resulting in a continuous biomass decay, termed endogenous decay (Rittmann and McCarty 2001).It is typically accepted to use a simple fi rst-order rate law r decay = k d M where k d is a decay coeffi cient (1/time).In some multispecies models, a portion of the biomass that decays accumulates into a pool of inert biomass.There are other mechanisms for biofi lm loss that result from discrete events over time, including cell lysis and detachment.These will be briefl y discussed later.Typically, two different strategies have been used for biomass spreading, namely continuous and discrete models (Fig. 2).The fi rst approach treats biomass as a separate continuum phase (the fl uid being the other phase); the bacteria grow in response to substrate in the fl uid phase and this growth leads to an expansion into the fl uid domain.In discrete models, the biomass is represented as discrete quantities, either assigned to numerical grid blocks or to particles.In the discrete approach, bacteria consume substrate and grow, but the spreading process uses rules to spread the biofi lm.Attachment and these spreading approaches will be described in the next section.Latt ice Boltz mann-Based Approaches for Pore-Scale Reactive Transport 401

LATTICE BOLTZMANN METHODS FOR FLOW AND REACTIVE TRANSPORT LBM for fl uid dynamics
In LBM, a fl uid is treated as mesoscale particles moving in a lattice domain (space discretization) with discrete lattice velocities (momentum discretization) at discrete time steps (time discretization).The fl uid variables (such as density and velocity) can be obtained through the moments (summation) of distribution functions over all discrete velocities.Typically the two-dimensional nine-velocity (D2Q9) and the three-dimensional nineteen-velocity (D3Q19) models for the 2-D and 3-D space, respectively, are employed (Fig. 3).The most popular LB fl uid fl ow model is the Bhatnagar-Gross-Krook (BGK) or single relaxation time (SRT) model (Bhatnagar et al. 1954).However, a viscosity-dependent permeability is usually obtained when adopting SRT-LBM for simulating fl uid fl ow in porous media (Pan et al. 2006).In order to overcome such defect, the multi-relaxation-time (MRT) model has been proposed for fl ow simulation (d 'Humières et al. 2002;Pan et al. 2006).For a general purpose, we choose a 3-D multiple-relaxation time (MRT)-LBM in this article, because the MRT-LBM model has been shown to be superior to the SRT method (Pan et al. 2006;Hilpert 2011).
Multiple Relaxation Time LBM.The MRT-LBM model transforms the distribution functions in the velocity space of the SRT-LBM model to the moment space by adopting a transformation matrix.In the SRT-LBM model, the evolution equation for the distribution functions (i.e., LBE) can be written as eq ( , where f i (x,t) is the i-th density distribution function at the lattice site x and time t.S is the relaxation matrix.For the D3Q19 model with N = 18, the discrete lattice velocity e i is given by 0 0 ( 1,0,0), (0, 1,0), (0, 1,0) 1 6 ( 1, 1,0), (0, 1, 1), ( 1,0, 1)  7 18.

402
Yoon, Kang & Valocchi f eq is the i-th equilibrium distribution function and is a function of local density and velocity with the weight coeffi cient w i as w i = 1/3 for i = 0; w i = 1/18 for i = 1, 2,…, 6; w i = 1/36, for i = 7, 8,…,18. 1 / 3 s c  is the speed of sound.By multiplying both sides of Equation ( 17) with a transformation matrix Q, an (N + 1) × (N + 1) matrix, the evolution equation in the moment space can be expressed as where m and m eq are the velocity moments and equilibrium velocity moments, respectively.Q -1 is the inverse matrix of Q, both of which are given in d 'Humières et al. (2002).The transformation matrix Q is constructed based on the principle that the relaxation matrix  , S an (N + 1) × (N + 1) matrix, in moment space can be reduced to the diagonal matrix, namely There is a large degree of freedom in choosing the relaxation parameters.According to Pan et al. (2006), one optimal set of relaxation parameters is given as: where  is related to the fl uid viscosity by 2 0.5, which is called the two-relaxation-time (TRT) model (Ginzburg and d'Humières 2003).
The detailed equilibrium velocity moments m eq can be found in Pan et al. (2006) and density and momentum (j) are determined by , The expression in Pan et al. (2006) includes the mean density of the fl uid ( 0 ), which is employed to reduce the compressibility effects of the model.Equations (20-25) can recover the Navier-Stokes equations using Chapman-Enskog multiscale expansion under the low Mach number limitation.In contrast to the SRT model, the collision step in the MRT model is implemented in the moment space, while the streaming step is carried out in the velocity space.
Boundary conditions.As described earlier, complex geometries (boundaries) of porous media can be easily treated.A no-slip boundary condition at the solid surface can be realized through the "bounce-back" condition, which mimics the phenomenon that a particle refl ects its momentum when colliding with a solid surface.It has been recognized that the combination of the SRT-LBM model with a standard bounce-back (SBB) scheme for fl uid-solid boundaries may lead to viscosity-dependent permeability (Pan et al. 2006), causing errors in boundary locations.Recently, several schemes have been developed to more accurately represent such boundaries using spatial interpolations, including the linearly interpolated bounce-back (LIBB) scheme, the quadratically interpolated bounce-back (QIBB) scheme, and the multirefl ection (MR) scheme (Bouzidi et al. 2001;Ginzburg and d'Humières 2003;Yu et al. 2003).These methods have been incorporated separately into SRT-LBM and MRT-LBM, and a systematic quantitative comparison of the numerical accuracy and convergence rate of these methods for complex porous medium systems has been reported by Pan et al. (2006).In their work, the MRT models signifi cantly reduce the viscosity dependence of permeability in porous medium simulations.In the MRT models, the location where the exact fl ow boundary conditions are satisfi ed can be essentially viscosity-independent by specifying an appropriately constructed collision operator, but it is impossible to achieve this with an SRT model.

LBM for multi-component reactive transport
Advection-diffusion reaction for a single-component system.The LB models for chemically reacting fl uid fl ows were fi rst introduced by Kingdon and Schofi eld (1992) and Dawson et al. (1993).In their models, the LB equations for transport have a similar form to the fl ow equation, with the addition of a source/sink term representing chemical reactions.The chemical reactions used in Kingdon and Schofi eld (1992) represent the Selkov model.In a more general case, homogeneous chemical reactions taking place in an aqueous fl uid can be written as a mass action law as in Equation ( 9).If the concentrations of the aqueous species are assumed to be suffi ciently low so that their effect on the density and velocity of the solution is negligible, then the reactive transport of solute species can be described using another set of distribution functions, k g  , which satisfi es a similar evolution equation as where I r is the reaction rate of the r-th reaction, C k is the solute concentration of the k-th species, and  k is the relaxation time related to the diffusivity.
It is noted that the advection-diffusion equation obeyed by the concentration is linear in velocity u.This indicates that the equilibrium distributions for scalar transport need only to be linear in u; and thus lattices with fewer vectors are suffi cient for scalar transport.Thus, a reduced D2Q5 lattice model in which the discrete velocities (i = 5, 6, 7, 8) in the D2Q9 model (Fig. 3) are not considered, is often used for solute transport.The D2Q5 model is combined with an equilibrium distribution function that is linear in u (Sullivan et al. 2005 where the concentration C k and a constant J k are given by The rest fraction J 0 can be selected from 0 to 1.In the literature, different forms of equilibrium distribution functions are adopted.The equilibrium distribution function given by Equation ( 29) is a general formula, which becomes the one used in Huber et al. (2008) if  J 0 = 1/3 and Kang et al. (2007) if J 0 = 0.The accuracy and effi ciency of the reduced D2Q5 model has been confi rmed in the literature (Kang et al. 2007;Huber et al. 2008;Chen et al. 2012Chen et al. , 2013aChen et al. , 2014a)).Equation ( 29) covers a wide range of diffusivity by adjusting J 0 , which is a prominent advantage of such an equilibrium distribution (Sullivan et al. 2005;Chen et al. 2013b).Using the Chapman-Enskog expansion technique, one can prove that Equation ( 26) recovers the pore-scale advection-diffusion-reaction equation for an incompressible fl ow fi eld as in Equation ( 6).The physical diffusion coeffi cient using the LB parameters is given by Therefore, one can solve the reactive transport problem by solving Equation ( 26) for each species, assuming all reaction rate constants are known.
LBM for reactive species has been used for various problems.Dawson et al. (1993) simulated pure diffusion, homogeneous chemical reactions and pattern formation due to Turing instability.Their LB simulation results agreed well with theoretical predictions and captured the basic physics as predicted by the macroscopic reaction-diffusion equations.Weimar and Boon (1996) simulated an athermal fl ow that advects reactant species.They considered a nonlinear reactive system described by the Brusselator model.Yan and Yuan (2001) simulated the Belousov-Zhabotinskii reaction and showed that the LBM could capture the well-known chemical clock of the diffusion-reaction system.Recently it has also been used to simulate bacterial chemotaxis (Hilpert 2005) and bacterial growth in porous media (Zhang et al. 2010b).

Advection-diffusion reaction for multi-component system.
As described earlier, solving Equation ( 26) for each species corresponding to Equation ( 6) may be very computationally expensive.For multi-component systems, Kang et al. (2006) have derived the following LB equation for the total primary species concentrations, j  (Eqn.7) with reactions written in canonical form ,( 1,..., ), where 1,..., C j N  and C N is the number of primary species, j G  is its distribution function along the  direction, eq j G  is the corresponding equilibrium distribution function and aq  is the dimensionless relaxation time for all the aqueous species.With the homogeneous reactions in the canonical form, formulating a LB equation for total concentration j  and replacing the rates of these reactions with mass action equations, the number of unknowns and evolution equations is reduced from N C + N R to N C .The reduction can be signifi cant for a system with many aqueous species.For example, the chemical system of Na + -Ca 2+ -Mg 2+ -H + -SO 4 2--Cl --CO 2 with the reaction of calcite to form dolomite and gypsum has (N C + N R ) equal to 23, but N C equal to only 7 (Kang et al. 2010).Here only speciesindependent diffusion is considered, guaranteeing conservation of charge in the aqueous phase as described earlier.Different diffusion coeffi cients can be obtained by varying x  or t  in Equation (30).More information on LBM simulation of electrochemical systems that includes species-dependent diffusion can be found in He and Li (2000) and will be briefl y introduced later.Equation ( 31) can recover the pore-scale advection-diffusion equation without reaction terms in Equation ( 7).

LBM for chemical reactions at the fl uid-solid interfaces.
Early work on chemical reactions at solid surfaces was based on LGA.Wells et al. (1991) pioneered an LGA model that couples solute transport with chemical reactions at mineral surfaces and in pore networks.In their model, chemical reactions considered at solid surfaces included precipitation/dissolution, sorption and catalytic reaction.Dissolution and precipitation reactions were simulated by allowing wall nodes to serve as sources or sinks for mass of a dissolved component.Whenever a particle collides with a wall, a unit of mass may be exchanged, thus increasing or decreasing the local concentration in solution depending upon the saturation state of the fl uid.Later, Sullivan et al. (2005Sullivan et al. ( , 2006) simulated 2-D and 3-D packed bed reactors using a LBM.They accounted for the local fl uid velocity by including reaction via concentration changes in the equilibrium distribution.Recently, Huber et al. (2014) treated heterogeneous reactions at the fl uid-mineral interfaces which are explicitly part of the computational domain.So the internal fl uid-solid boundaries are treated in a similar way as any fl uid site in the computational domain, and the algorithm to solve for surface reactions is independent of the surface shape and orientation of the grains.
Typically heterogeneous reactions in the LBM are treated through boundary conditions.Kelemen et al. (1995) extended the LBM with a dissolution boundary condition such that the fl uid particle colliding with the wall has a probability of detaching a solid particle.Verhaeghe et al. (2006) designed boundary conditions to impose a concentration or a fl ux on a solid interface for use in multicomponent LBMs.Other approaches for concentration boundary treatment focused on deriving the unknown distribution functions at the fl uid-solid interface based on local conditions for different boundary condition types.To discuss the concentration boundary conditions, a schematic illustration of the fl uid-solid interface is presented in Figure 4. Following Zhang et al. (2012), the boundary condition at this interface is  Yoon, Kang & Valocchi Kang et al. (2007) proposed a new expression of the concentration distribution function in terms of the corresponding concentration and its gradient .
After each streaming step in Figure 4, g 2 is unknown and g 4 is known; and g 1 and g 3 do not affect the fl uid domain and hence are not needed to calculate their values.In this case, the reactive wall is static.They pointed out that since g 2 enters the domain and g 4 leaves the domain, the following equation based on Equation ( 33) is obtained at the reactive node R 2 4 1 , where c = dx/dt and C can be the total concentration for the multi-component system.The fact that the non-equilibrium portion of the distribution functions in opposite directions takes the opposite sign for a static wall derives another equation eq neq eq neq eq eq 2 4 Finally, the unknown distribution g 2 can be obtained by combining Equations (34-35).
Instead of using Equation ( 33), Zhang et al. (2012) proposed to directly solve Equation ( 32) to obtain the concentration at the boundary using the following difference scheme where C R and C F are the concentrations at interface node R and adjacent fl uid node F, respectively, and Δx is the distance between nodes F and R. In Equation (36) C R is the only unknown and thus can be directly solved.The unknown g 2 can be solved with Equation (35).Following Zhang et al. (2012), Chen et al. (2013b) improved the accuracy of boundary treatment by adjusting the bounce-back particle distribution functions according to the velocity and concentration values at the midpoint of a boundary lattice link.
In a recent study, Chen et al. (2013b) proposed to use Equation ( 33) to solve the unknown distribution function after C R is obtained from Equation (36).For the boundary node shown in Figure 4, the new boundary condition can be written as Equation ( 37a) is used to solve C R , followed by using Equation (37b) to solve g 2 .This boundary condition treatment can handle moving boundary conditions.The boundary condition described by Equation ( 37) can be uniformly used to treat both a reactive fl uid-solid boundary and a zero-fl ux fl uid-fl uid boundary.For the static fl uid-solid boundary, a precipitation or dissolution reaction leads to the following reduced form of Equation ( 37) Latt ice Boltz mann-Based Approaches for Pore-Scale Reactive Transport 407 For the moving fl uid-fl uid boundary, there is no reaction, and the concentration gradient of solute is zero.Thus Equation (37) reduces to

Update of solid-pore geometry
To model accurately reactive transport involving signifi cant mass transfer between solids and fl uids due to dissolution and/or precipitation, it is necessary to account for the time evolution of the solid phase (pore geometry).For tracking the fl uid-solid interfaces, the phase fi eld (PF) method (Chen 2002) and the cellular automaton (CA) methods (Sun et al. 2009) are widely used.Other methods such as Volume of Fluid (VOF) and Level Set (LS) commonly used for multi-fl uid phase fl ow, can be also adopted for this purpose (Li et al. 2008).Verberg andLadd (2000, 2002) designed an algorithm for simulating chemical erosion in rough fractures using an optimized LB scheme.A continuous bounce-back scheme allows for the boundary to be located anywhere between two grid nodes.The new solid structure is determined based on the local fl ux of tracer particles across the solid surface, assuming an instantaneous reaction, so dissolution is diffusion-controlled.In Verhaeghe et al. (2006), the amount of species injected in the system was calculated from the difference between populations leaving and entering the system.Li et al. (2008) studied a similar problem using LS to track the fl uid-solid interfaces.In addition, Luo et al. (2012) implemented a model using a diffuse interface method to track the fl uid-solid interface.Most recently, Huber et al. (2014) combined LBM and PF method for dissolution-precipitation processes involving single or multiple species.
To track the moving fl uid-solid interfaces caused by dissolution and precipitation, Kang et al. (2003Kang et al. ( , 2006) ) developed the volume of pixel (VOP) method, a type of CA method.Each fl uid-solid interface node represents a control volume (a control area in the 2-D case) with a size of 1 × 1 × 1 (in lattice units) and is located at the center of this volume.As in Figure 4, the node R is the center of the control volume surrounded by dashed lines.The control volume is assigned with a certain amount of mineral mass (volume) that changes with time due to chemical reaction (Eqn.12); then the volume is updated using Equation (13).In these studies, both a time interval and specifi c surface area (a m ) equal unity in lattice units.When V m reaches certain threshold values, the pore geometry needs to be updated.Here we will briefl y describe the volume of pixel method developed by Kang and coworkers.
Single mineral VOP method.In this method, the grid size is assumed to be small enough to represent each node by one mineral at one time and the effect of both dissolution and precipitation is recorded at that node through Equation ( 13).For dissolution, the solid node associated with V m can be simply removed (i.e., changed to a pore node) when V m reaches zero (Kang et al. 2002(Kang et al. , 2014;;Chen et al. 2014a).For precipitation, when V m reaches a certain threshold value, a nearby neighboring fl uid node (child node) can be chosen to become a solid node.There may be multiple such nodes, and different ways of choosing the child node may lead to different crystal growth patterns.Therefore, the single-mineral VOP model has an extra degree of freedom which allows for incorporating different crystal growth mechanisms such as a random-growth method for single species (Kang et al. 2003(Kang et al. , 2004) and multi-component systems (Kang et al. 2006).In this random-growth method, the growth has no preference in any particular direction and the method has been shown to be lattice-effect free.Figure 5 shows crystal growth patterns from a supersaturated solution of a single species using the randomgrowth method.It is clear that crystal shape in this case (d) is fairly round.The VOP has many advantages: a clear physical concept, simple and stable arithmetic, easy implementation of various surface reactions, and fl exible coupling with different nucleation and crystal growth mechanisms.The method has been used successfully to predict many moving solid-fl uid interface phenomena, such as crystal growth (Kang et al. 2004), rock dissolution due to acid injection (Kang et al. 2002(Kang et al. , 2003)), Liesegang bands or rings (Chen et al. 2012), and dissolution and precipitation involved in CO 2 sequestration (Yoon et al. 2012).
Growth in preferred directions can be achieved by exploiting the extra degree of freedom in the VOP model and incorporating corresponding physics in the growth rules.For examples, polygonal and dendritic crystals with symmetry have been produced by aligning the direction of growth to that of the maximum concentration gradient and different morphologies have also been obtained by varying the probability of adding a solid node in that direction (Lu et al. 2009).Another growth mechanism is epitaxial growth in which the precipitates quickly spread and cover the surface of the primary mineral, thus only a small amount of precipitates can completely stop the dissolution (Prieto et al. 2003).
Multiple mineral VOP method.In reality, changes of solid morphology can involve scales much smaller than the lattice or pore size used in the simulations.Therefore, it is reasonable to assume that each node can be represented by multiple minerals with the initial total volume fraction equal to unity.In this case, the volume fraction of each mineral is still updated by Equation ( 13) and dissolution recorded in the solid node.Mass accumulation due to precipitation, however, is recorded at the neighboring pore nodes.As seen in Figure 4, when the node F becomes a solid node composed of multiple minerals.In contrast to previous methods (Kang et al. 2003(Kang et al. , 2004(Kang et al. , 2006) ) for updating pore geometry, the multiple-mineral VOP method can account for coexistence of multiple minerals at the same node and the growth of minerals is deterministic rather than random.Figure 6 shows crystal structures at Da = 600 with the single-mineral VOP (random-growth) method (left) and multiple-mineral VOP approach (right) to update the solid phase (Kang and Lichtner 2013).At this high Da number, the process is diffusion-controlled.As a result, the crystal structure is not compact.However, while the crystal on the left is an open cluster-type structure, consistent with that of the multi-particle, diffusion-limited aggregation in the simulation of solidifi cation structures of alloy melt, the crystal structure on the right is highly symmetric and regular.The multiplemineral VOP approach was also used in the pore-scale study of reactive transport involved in CO 2 sequestration (Kang et al. 2010).Clearly, when combined with appropriate methods for the underlying physical properties and processes to update the solid phase, the LBM can simulate a variety of reactive fl ow problems.

LBM for biofi lm dynamics
The system of equations described in Equations ( 3), ( 6) and (14-16) includes three dynamic processes, which occur at very different time scales: biofi lm development on the order of hours or days; transport and reaction of chemical species on the order of minutes; and fl uid fl ow on the order of seconds (Picioreanu et al. 2000;Pintelon et al. 2009).Therefore, when the fl ow or substrate concentration suddenly changes, the biofi lm can be assumed unchanged (i.e., in a "frozen" state).Based on this assumption, the three processes are commonly decoupled and solved sequentially.As described earlier, the biofi lm consists of multiple active microbial species, inert biomass, and EPSs.Biofi lm spreading is typically computed using a continuum or discrete model as shown in Figure 2. Since a majority of pore-scale biofi lm models in literature adopt a discrete representation of the biomass, this approach will be explained in Figure 6.Crystal structures and solute concentration obtained using single-mineral VOP approach (random growth) (left) and multiple-mineral VOP approach (right).[Used by permission of Hindawi, from Kang QJ, Lichtner PC (2013) A lattice Boltzmann method for coupled fl uid fl ow, solute transport, and chemical reaction.In: Progress in Computational Physics.Ehrhardt M (ed) Bentham Science Publishers, more detail below.For continuous model of biofi lm spreading, the reader may refer to key literature for various variants (Alpkvist and Klapper 2007;Wang and Zhang 2010;Lindley et al. 2012;Zhang and Klapper 2014).The LBM earlier described for hydrodynamics and reactive transport can be applied for biogeochemical reactive transport.In fact, the LBM for advection-diffusion-reaction (Eqns.27-30) developed by Sullivan et al. (2005) for packed bed reactors was directly applied for biofi lm dynamics (von der Schulenburg et al. 2009).In Sullivan et al. (2005), the diffusion coeffi cient within intra-particle space was adjusted by changing the relaxation time ( k ) in Equation ( 30), which can be adapted for diffusivity within biofi lm.Knutson et al. (2005) simulated biofi lm growth in porous media to account for biofi lm morphology observed in microfl uidic experiments using the LB-FVM with a CA algorithm for biofi lm dynamics.von der Schulenburg et al. ( 2009) presented the fi rst numerical 3-D pore-scale model of biofi lm growth in porous media, based on an LB simulation platform complemented with an individual-based biofi lm model (IBM), where the biofi lm is represented as a collection of non-overlapping hard spheres.Following this work, Pintelon et al. (2009) added a biomass detachment technique using a fast-marching level-set method that modeled the propagation of the biofi lm liquid interface with a speed proportional to the adjacent velocity shear fi eld.
Attachment.The fi rst stage in biofi lm development is its attachment to the solid surface.This is a complex process that involves deposition of suspended bacteria onto the surface followed by physiologic changes (Breyers 1988).The deposition process has been quantifi ed using colloid fi ltration theory (e.g., Harvey and Garabedian 1991; see other citations in Rittmann 1993), but most of these studies address migration of bacteria in porous media.For pore-scale biofi lm dynamics, some initial spatial distribution of attached biomass is typically assumed.For example, in the approach of Knutson et al. (2005) and Tang et al. (2013), the initial mass of microbes was randomly distributed among grid cells along the liquid-solid interface, a distribution which simplifi ed the inoculation of biomass in the microfl uidic pore network.A similar approach to inoculation is used by von der Schulenburg et al. ( 2009), whereas Bottero et al. ( 2013) modeled a fi xed rate of attachment by converting grid cells from fl uid to biomass at the liquid-solid interface at each time step; they used a random spatial pattern, but screened out solid boundaries on particles in stagnant zones.Picioreanu et al. (1998), a popular approach to a discrete representation of the biomass is to assign biomass to fi xed grid cells (normally squares in 2-D and cubes in 3-D).A mass balance equation applies to each biomass cell (Eqn.14).The biomass in each cell grows due to substrate utilization, and when the total biomass in the cell exceeds a critical value, the excess biomass is spread (transported) to adjacent grid cells according to some rules such as CA as described for mineral precipitation.Many different CA rules for biomass spreading have been used in the literature; the three most common are briefl y discussed here (Tang et al. 2013).In the fi rst rule (e.g., Kreft et al. 2001), when the biomass concentration in a grid is higher than the critical value, the biomass in the source grid divides in half, similar to bacterial reproduction.One half stays in the grid, the other half displaces a randomly chosen nearest-neighbor, which in turn displaces another randomly chosen neighbor, and so on until a free liquid grid (i.e., a destination grid on the biofi lm surface) is available.In the second rule (e.g., Noguera et al. 1999), the excess biomass (i.e., biomass that is greater than the maximum) is distributed directly to one of the nearest free liquid grids (i.e., destination grid) non-adjacent to the source grid cell.In the third rule (e.g., Tang and Valocchi 2013), excess biomass spreads by pushing a line of grid cells that are on the shortest path from the source grid cell to the destination grid cell, causing changes to the fractions of different biomass species in the grid cells on the path.The three different rules yield similar simulated results for a single species biofi lm, but different results for multi-species biofi lm (Tang and Valocchi 2013).These rules can be modifi ed in many ways; for example, instead of distributing excess biomass to a randomly chosen nearest neighbor, the destination grid having the highest substrate concentration can be selected.Knutson et al. (2005) and Tang and Valocchi (2013) exclude destination cells where fl uid velocity or fl uid shear stress is above a certain threshold, and found that this exclusion was needed to qualitatively match biofi lm spreading observed in microfl uidic experiments (Nambi et al. 2003;Zhang et al. 2010b).

Discrete model for biofi lm growth and spreading. Following the pioneering work by
Even with simple CA rules (e.g., Picioreanu et al. 1998), the interplay between mass transfer and biofi lm growth can lead to a diverse biofi lm patterns as shown in Figure 7. Initially, 20% of the surface has a randomly spaced inoculum; a limiting substrate diffuses from the upper boundary to the lower surface resulting in biofi lm growth.The irregular "mushroom shaped" colonies form when the rate of diffusion is small relative to the rate of degradation at 50 days (Fig. 7).In the other discrete model (i.e., IBM), however, spherical particles grow by consuming substrate, resulting in an increased radius.When the particle mass exceeds a threshold, it divides into mother-daughter particles that spread out; a "shoving" algorithm is used to eliminate overlapping radii of particles (see Kreft et al. 2001 andPicioreanu et al. 2004 for details).Since the IBM is not constrained to grid coordinates, biofi lm morphologies tend to be more rounded than those resulting from the CA approach.The latter also leads to biofi lms that are artifi cially constrained to the underlying grid.The IBM is fl exible because it can include different bacteria species along with inert biomass and other components (Picioreanu et al. 2004).
Biofi lm loss processes.In addition to the endogenous decay process described earlier, there are several other mechanisms for loss of mass from biofi lms, including lysis, protozoan grazing, and detachment.Detachment can include several mechanisms, including erosion of small pieces from the outer surface of the biofi lm and sloughing of larger sections (IWA Task Group 2006).A recent trend is to use a relatively rigorous bio-mechanical model of biofi lms to study detachment and related processes (Dupin et al. 2001;Klapper and Dockery 2010;Lindley et al. 2012).Several studies adopting a discrete biofi lm model use mechanical principles to simulate detachment in porous media systems.Picioreneau and co-workers treat the biofi lm as an elastic solid and use classical mechanical equilibrium and compatibility conditions to determine "failure," i.e., detachment (e.g., Picioreneau et al. 2001).This approach requires input of basic biofi lm mechanical properties such as the tensile strength and elastic modulus, and is relatively straightforward to implement for 2-D CA biofi lm models.von der Schulenburg et al. ( 2009) neglected to include detachment due to relatively lower shear rates, but Pintelon et al. (2009) used the IBM for biofi lm spreading and used a fastmarching level-set (FMLS) method to simulate biomass detachment.Bottero et al. ( 2013) used a similar detachment model based on level sets, but allowed the cohesive strength of biofi lm to decrease as biomass decays.It is one of the few studies that explicitly considers the cell lysis mechanism; biofi lm lysis occurs when the fraction of inert biomass in a grid cell exceeds a certain value.In all models, including Bottero et al. ( 2013), all detached biomass is assumed to be instantly removed from the system, an assumption which may not be accurate in porous media.Biofi lm detachment can also be considered in pore-network models (e.g., Stewart and Kim 2004;Thullner and Baveye 2008), but these models are not covered in this article.

Multiphase reactive transport
There has been little work on numerical simulation of chemical processes coupled with multiphase fl ow at the pore scale (Meakin and Tartakovsky 2009;Parmigiani et al. 2011;Chen et al. 2013b).Although models for both single-phase multi-component reactive transport and multiphase fl ow have been developed, integrating the two is challenging because new physical processes absent in their respective individual system will arise due to the coupling between them.For example, the displacement patterns that range from planar fronts to capillary and viscous fi ngers will determine the fl uid-fl uid interfacial area.This interfacial area will affect chemical reactions at the interface and will change the chemical heterogeneity of the system.This chemical heterogeneity, together with the physical heterogeneity of the porous media, in turn, will infl uence the intrinsically complex contact line/contact angle dynamics, and hence infl uence the displacement patterns.
Recently, several works attempted to investigate chemical processes coupled with multiphase fl ow using the LBM.Parmigiani et al. (2011) studied a non-wetting fl uid fl owing into a wetting fl uid coupled with dynamic evolution of the solid geometries due to chemical reactions.Although it was suggested that their model can be used for reactive transport with dissolution or precipitation, the demonstration was limited to only melting of the solid phase.Chen et al. (2013b) combined the LBM and the VOP to simulate multiphase reactive transport with phase transition and dissolution-precipitation processes.This pore-scale model can capture coupled non-linear multiple physicochemical processes including multiphase fl ow with phase separation, mass transport, chemical reaction, dissolution-precipitation, and dynamic evolution of the pore geometries.The model was used to study the thermal migration of a brine inclusion in a salt crystal.However, the pore-scale model was limited to the liquid-gas twophase fl ow within a single-component water-vapor system, and the transport of only a single solute was considered.To overcome these problems, Chen et al. (2014b) developed a twophase multi-mixture model based on the pseudo-potential LB model (Shan and Chen 1993), the mass transport LB model and the VOP model for multi-component multiphase reactive transport with dissolution-precipitation processes.After the two-phase system is constructed using a unifi ed liquid component and a unifi ed gas component, the mass transport of each component in the corresponding phase is solved using the mass transport LB model.The dynamic evolutions of the liquid-solid interfaces due to dissolution-precipitation reactions are captured using the VOP method.

Electrokinetic transport
Electrokinetic fl ow, which involves multiple processes including fl uid fl ow, electrostatic interaction, species diffusion, and sometimes energy transfer, is pervasive both in nature and in engineered systems (Li 2004;LinkMasliyah 2006).Although numerous theories and models for large-scale electrokinetic fl ows have been developed for almost a century (Bancroft 1926;Li 2004), only in recent decades has the electrokinetic transport been applied at micro-and nano-scales.The Poisson-Boltzmann (PB) model, composed of a Navier-Stokes equation for fl uid fl ow and a PB equation for electric potential distribution, has been widely used for analysis and prediction of electrokinetic fl ows in microchannels (Li 2004;Karniadakis et al. 2005; Latt ice Boltz mann-Based Approaches for Pore-Scale Reactive Transport 413 Wang and Kang 2009;Wang et al. 2010a,b).The Boltzmann distribution assumption of ions in electric double layer (EDL) leads to the decoupling between the fl uid fl ow and the ion distribution, simplifying the predictions.The lattice spacing can be based on the Debye length representing an electrical double layer thickness (~1-20 nm for typical aqueous solutions), which was suffi cient to represent the main electrolyte fl ow within the microchannel (e.g., Wang et al. 2008).Many algorithms and analyses based on the PB model have been proposed for understanding electrokinetic fl uid mechanics and optimizing electrokinetic microdevices (Holst et al. 2000;Li 2001;Guo et al. 2005;Wang et al. 2006;Chen et al. 2007;Shi et al. 2008).However, the Boltzmann distribution is an equilibrium model with questionable applicability when ions have macroscopic motion (advection or diffusion), surface is heterogeneously charged, or the bulk solution boundary condition is not satisfi ed in the place suffi ciently far away from the charged surface.In contrast, the dynamic model uses the Nernst-Planck equation (Eqn.5) for ion transport instead of the hypothesized Boltzmann distribution (Levich 1962).In the dynamic model, each transport process is governed by the corresponding fundamental dynamic equation.The dynamic model has been used to validate the PB model by a few researchers (Qu and Li 2000;Yang et al. 2001;Fu et al. 2003;Kang and Suh 2006;Baldessari 2008), yet it was reported that the PB model failed in the case of overlapped double layers (Qu and Li 2000;Kang and Suh 2006).
Because the transport processes infl uence each other, the governing equations in the dynamic model are coupled together, posing a great challenge to the numerical solution.Wang and Kang (2010) developed a numerical framework to solve the dynamic model to study the applicability of the PB model for electrokinetic fl ows in microchannels using coupled LBMs.The governing equation for each transport process is solved by an LBM and the entire process is simulated iteratively.Later, Yang et al. (2014b) developed an MRT-LBM for the Nernst-Planck (NP) equation to solve the double-diffusive-convection equation having the cross diffusion effects in a channel of a 0.8 μm width, and demonstrated that the MRT-LBM is more accurate than an SRT-LBM for the NP equation.

Coupled LBM-DNS for multicomponent reactive transport processes
Due to accurate, fi ne resolution representation of pore structures and fl ow fi eld, porescale reactive transport models including LBM are often computationally expensive.The computational expense of these pore-scale models can be highlighted when the reaction occurs locally, compared to the size of model domain of interest.To overcome the computational expense of pore-scale models, a variety of hybrid methods have been developed by coupling different scale models (e.g., Weinan et al. 2007).Here we briefl y introduce the recent development of coupled methods involving LBM.
LBM now has expanded capabilities for pore-and continuum-scale coupling with continuum-based methods (e.g., FVM, FDM, FEM).In the multiscale formulation, the key to improved numerical stability, accuracy, and computational effi ciency is a coupling strategy to transfer information between two neighboring regimes of different methods.Luan et al. (2010) developed an analytical expression (i.e., reconstruction operator (RO)) to transfer the velocity information at the continuum scale obtained using FVM to the single-particle distribution function for LBM within the overlapping region.Furthermore, Luan et al. (2011) coupled a local grid refi nement of LB (multi-block LBM) with FV for multiscale simulations and successfully demonstrated the accuracy and effi ciency of the multi-block LBM for fl uid fl ows around complex solid geometries.A coupled LB-FV method also has been adopted for different physical processes.For example, Zarghami et al. (2012) developed a cell-centered FV method for discretizing the SRT to improve the stability of the fl ow fi elds over a range of Reynolds numbers (Re).More recently, coupled LB-FV methods have been developed for reactive transport processes for multiphysics problems.Chen et al (2013a) used a multiscale modeling framework with a coarse (FVM region)-fi ne (LBM region) grid system to simulate electro-chemical transport reaction in a membrane fuel cell system.They adopted RO methods (Luan et al. 2010(Luan et al. , 2011) ) to transfer density, velocities, and concentration at the continuum scale (FVM) to distribution functions for density and concentration (e.g., Eqn.17, 26) at the pore scale in the LBM and achieved a maximum grid size ratio of 10:1 between FVM and LBM regions.Following this work, a generalized RO methodology was proposed for fl uid fl ow, heat transfer, and reactive transport processes (Chen et al. 2013b).A unifi ed coupling scheme for coupled LB-FVM was recently developed for unsteady fl uid fl ow and heat transfer by deriving a generalized form of RO (Tong and He 2015) where the coupling fl exibility was generalized for multi-physics problems.
Hybrid approaches coupling LB for fl ow and FD or FV for reactive transport have also been developed for mineral precipitation and dissolution in porous media (Yu and Ladd 2010;Yoon et al. 2012).Yu and Ladd (2010) demonstrated the movement of the solid-fl uid interface due to dissolution by constructing a piecewise continuous surface which can remap the surface to the regular gird system.Yoon et al. (2012) also showed that the concentration fi eld can be solved as a steady solution and FVM can be computationally effi cient for the concentration fi eld.This hybrid approach validated the update algorithm of the solid phase against experimental data obtained in a micromodel and it demonstrated that coupling LB-FVM for reactive transport has the potential of using advanced reactive transport codes such as PFLOTRAN (Lichtner et al. 2014) with a proper interface module.Recently Patel et al. (2014) coupled the LBM for solute transport with the generic geochemical code, PHREEQC (Parkhurst and Appelo 2013).The heterogeneous reactions at the mineral-fl uid interface are treated as pseudohomogeneous reactions by incorporating an additional source (i.e., collision) term computed from the PHREEQC to the fl uid node adjacent to a mineral node.This approach was validated with a set of benchmarks including the impact of pH, surface area, and arrangement of mineral grains on dissolution rate and pore structure evolution.In addition, Boek et al. (2014) recently developed a new hybrid LB and particle tracking method to study fl ow and reactive transport in porous media.The hybrid method used a graphical process unit (GPU) algorithm for large scale LB calculations to enhance computational effi ciency.Particle advection was then solved with LB fl ow fi elds using a second-order predictor-corrector scheme, and particle diffusion with a random walk followed by heterogeneous reaction was solved.

THREE-DIMENSIONAL CHARACTERIZATION OF PORE TOPOLOGY
Over the past two decades multiscale imaging methods for geomaterials such as X-ray computed microtomography (micro-CT), laser scanning confocal microscopy (LSCM), focused ion beam-scanning electron microscopy (FIB-SEM), and energy dispersive X-ray spectroscopy (EDS) have been tremendously improved.Major advances in 3-D imaging techniques can be traced back to early 1990s when methods for characterizing geomaterials at the microscopic scale using confocal microscopy and micro-CT was established (Spanne et al. 1994;Fredrich et al. 1995;Coker et al. 1996;Coles et al. 1998).The most predominant 3-D imaging method used for geomaterials is a non-destructive micro-CT technique.Micro-CT has been used to measure the geometric properties of pore topology and multiple fl uid phases at an order of microns.Theoretical and detailed technical aspects of micro-CT have been extensively reviewed by Wildenschild and Sheppard (2013) who also highlighted issues associated with the imaging process.Recently, 3-D micro-CT images of a real soil at 25 μm resolution were used to identify a mixed medium containing both pore and porous regions where small pores below the image resolution were treated as porous regions (Yang et al. 2014a).The mixed medium was then used to validate a unifi ed multi-scale model for porescale fl ow simulations in a mixed medium by applying the Stokes-Brinkman equations for variably saturated soils.
LBM has been extensively applied for a variety of complex fl uid fl ow problems, including single and multi-phase fl ow and reactive transport in complex geometries obtained from micro-CT imaging (Ladd 1994;Martys and Chen 1996;Keehm 2004;Fredrich et al. 2006;Boek et al. 2014;Jiang and Tsuji 2014;Gao et al. 2015).A major advantage of LBM is the relative ease of accounting for very complex pore-geometries from segmented image data due to its discrete dynamics (Chen and Doolen 1998).This method enables pore-scale simulations of fl uid fl ow and transport phenomena to enhance fundamental insights in porous media.Chen et al. (2007) and Zhang (2011) provide interesting reviews and fundamentals for single and multi-phase fl ows, respectively.For example, 3-D rendering of the pore phase in natural Castlegate sandstone at 1.67 μm resolution from micro-CT images is shown in Figure 8a (Fredrich et al. 2006).Permeability predictions from LB simulations for 1 to 3 mm 3 image volumes show a good agreement, over 4 orders of magnitude in scale, with values measured experimentally in the core samples (25.4 mm in diameter and 50.8 mm in height).This fi nding revealed the representative volume and length scales that are necessary to characterize geometrically complex porous media and predict fl uid transport properties at the macroscale.More recently, LB-based approaches were used to investigate the porosity dynamics due to precipitation and dissolution in 3-D pore structures from micro-CT images (Jiang and Tsuji 2014;Gao et al. 2015).These recent studies using micro-CT images can be combined with other imaging techniques (e.g., backscattered electron mapping coupled with X-ray based methods) for mineralogical distribution (Peters 2009;Landrot et al. 2012) to investigate coupled fl uid fl ow and geochemistry, as well as the emergence of coupled behaviors due to biogeochemical reactions.Furthermore, 3-D digital representations of rock from micro-CT images and physics modeling based on these pore structures provide the opportunity to further advance our quantitative estimates of permeability, elastic moduli, and electrical conductivity of several representative 3-D rock samples (Arns et al. 2002(Arns et al. , 2009;;Andrä et al. 2013a,b).
Recently, Yoon and Dewers (2013) rigorously quantifi ed the existence of a statistical representative elementary volume (SREV) using FIB-SEM imaging of nano-porous chalk using LB simulations with 3-D reconstructed chalk data as shown in Figure 8b.They highlighted the quantitative analysis of nano-pore structural features impacting pore-scale fl ow and transport properties and confi rm that the size of the SREV for a sample of Gulf Coast chalk can be established at ~10 m based on anisotropic permeability, tortuosity, and specifi c surface area.With such clear scale separation, upscaling averaging techniques, such as homogenization, permits the application of continuum physics principles to coupled thermohydro-mechanical-chemical (THMC) modeling.However, porous materials with a high degree of physical and chemical heterogeneity exist with no clear scale separation (e.g.Bilger et al. 2005).The problem of scale relating to the traditional FIB-SEM with gallium (Ga) ion milling (e.g., sample volumes in the order of ~10 × 10 × 10 μm 3 ) is now being confronted, which makes the FIB-SEM analysis more qualitative than quantitative.Newly emerging techniques such as the FIB with an inductively coupled plasma ion source (plasma-FIB) (e.g., Smith et al. 2014) and X-ray microscopy (XRM) or nano-CT (e.g., Gelb et al. 2011) will enable us to overcome sample size limitations and to characterize multi-scale hierarchical pore topology.

Pore cementation/dissolution and fl ow feedback
Here we summarize a few of LB-based works that emphasize coupling among fl uid fl ow, reaction, and fl ow feedback.Kang and co-workers have applied an LB pore-scale reactive transport model to study heterogeneous chemical reactions in a simple pore geometry under a variety of conditions (Kang et al. 2002(Kang et al. , 2006)).For example, Kang et al. (2006) demonstrated that the competition between diffusive and advective transport can alter permeability during the dissolution of different sizes of fracture aperture.More recently, Kang et al. (2014) applied a reactive transport LB model to study the change of porosity and permeability due to dissolution over a wide range of Péclet (Pe, the ratio of advective to diffusive transport) and Damköhler (Da, the ratio of reaction rate to diffusive mass transfer rate) numbers.They demonstrated that the permeability-porosity relationship is affected by different dissolution regimes in terms of Pe and Da numbers, but also by the specifi c porous medium structure.In particular, the relationship for the more geometrically complex porous medium shows more complexity than that for the simple fractured medium (Fig. 9).For the complex medium, a combination of a high Pe and Da (i.e., fast-fl ow and fast-reaction regime) results in wormholing and exhibits the fastest permeability increase.At a moderate Pe and high Da, a transition from a transportlimited dissolution regime to wormholing is observed for the complex medium as shown in Figure 9. Huber et al. (2014) applied a newly developed LB model with a new boundary scheme for heterogeneous reactions to study the permeability change of a porous medium due to dissolution and precipitation.They highlighted that relative permeability changes differ under dissolution and precipitation regimes.The permeability changes during dissolution were controlled by both reaction and transport (as the transport regime transitions from diffusion to advection).Increases in permeability correlate positively with Pe, but negatively with Da for a given amount of porosity increase.On the other hand, increases in precipitation correlate negatively with Pe, but positively with Da.This behavior suggests the presence of a hysteresis in the permeability and porosity relationship during cycles of dissolution and precipitation.Jiang and Tsuji (2014) developed a LB-FVM model for exploring the porosity and relative permeability changes caused by carbonate precipitation in 3-D digital microstructures obtained from micro-CT imaging of a Berea sandstone sample.Changes of pore structures caused by carbonate precipitation are shown in Figure 10.For the sake of simplifi cation, chemical reactions in aqueous solution were assumed in equilibrium and only Ca 2+ concentration was solved by the advection-diffusion model.All pore structures were assumed reactive surfaces for carbonate precipitation.The porosity was reduced from 0.21 to 0.103 after ~ 20 days of precipitation (Fig. 10).The porosity change corresponds to the reduction of permeability from 769 mD to 11 mD.In particular, simulations of multiphase fl ow in the altered pore structure demonstrated that reduced pore space (a fi nal case in Figure 10) affects the relative permeability of the nonwetting phase more signifi cantly than that of the wetting phase.This fi nding implies that mineralization may signifi cantly impact the injectivity of CO 2 in water wetting conditions during geologic CO 2 storage.

Microfl uidic experiments for pore cementation and fl ow blocking
Recent developments of in situ measurement techniques for pore-scale reactive transport experiments provide a unique opportunity to test and validate pore-scale modeling approaches.In particular, a coupled pore-scale reactive transport model was rigorously validated against well-controlled microfl uidic experimental results to demonstrate the importance of realistic pore confi gurations, fl ow and transport physics and geochemistry, and to predict how mineral precipitation alters fl ow paths by pore plugging for different chemical conditions (Yoon et al. 2012;Fanizza et al. 2013;Boyd et al. 2014;Oostrom et al. 2014).The coupled model consists of LBM for fl uid fl ow and FVM for reactive transport with mineral precipitation and dissolution.Heterogeneous reactions were solved using the boundary condition given by Equations (10-13).Pore-scale experiments led by Werth and co-workers (Zhang et al. 2010a;Fanizza et al. 2013;Boyd et al. 2014) were used.The micromodel consisted of a 2 × 1 cm pore network containing a pattern of cylindrical posts (300 m in diameter) with 180-m pore spaces and 40-m pore throats, a depth of ~20-40 m, and a porosity of ~0.39.The micromodel surface was treated thermally to have a thin layer (~100 nm's) of silicon dioxide, mimicking natural materials.Representative experimental observations are shown in Figure 11.for pore blockage were measured.Specifi cally, the coupled LB-FVM approach for pore-scale reactive transport problems includes calcium carbonate precipitation and dissolution kinetics (Yoon et al. 2012), uranyl phosphate precipitation to account for reaction mechanisms under different chemical conditions (Fanizza et al. 2013), and the impact of magnesium on calcium carbonate precipitation (Boyd et al. 2014).In particular, tracer tests (e.g., Fig. 11b) indicate that mineral precipitation along the center mixing zone led to substantial pore blockage under all experimental conditions.The observed changes in crystal shape, size, polymorph, and precipitation rate under different chemical conditions suggest that proper chemistry systems be considered in determining the impact of pore blockage on mixing and reaction.
In Yoon et al. (2012) a chemical system for CaCO 3 precipitation and dissolution includes fi ve primary components: Na + , Ca 2+ , H + , CO 3 2-, and Cl -, and a number of secondary species: OH -, HCO 3 -, and H 2 CO 3 (aq).Simulation results are shown in Figure 12.The reactive transport model included the impact of pH upon carbonate speciation and calcite dissolution using the form of the calcite rate law (Chou et al. 1989): where k m is the reaction rate constant in Equation ( 10), k 1 , k 2 , and k 3 are the empirically obtained reaction rate constants, and a i is the activity of species i.The calcite reaction rate is pH-dependent (pH < 7), but becomes constant above a pH of ~ 7.5.A similar chemical system for CaCO 3 precipitation and dissolution in the fl owing system was used in the literature (Emmanuel and Berkowitz 2005;Tartakovsky et al. 2008;Flukiger and Bernard 2009;Kang et al. 2010;Molins et al. 2012).In the above formulation, the key discrepancy between pore and continuum scales is the reactive surface area.As described earlier, heterogeneous reactions at mineral interfaces can be treated through the boundary conditions (Eqn.11).In particular, Yoon et al. (2012) developed an algorithm to update a specifi c surface area over time in the micromodel, depending on the orientation of CaCO 3 growth and dissolution.Updating the surface area properly was critical to accounting for the overall reaction rates as shown in Figure 12a.Hence, the specifi c surface area was updated at each iteration step.It was also found that the proper estimation of the dissolution rate was necessary to adequately simulate the decrease of precipitate area (i.e., enhanced dissolution) after ~ 18 min.In this particular case, the dissolution rate increased by a factor of 100.Based on Raman spectroscopy and image analysis, a higher dissolution rate was attributed to stability of CaCO 3 polymorphs including amorphous calcium carbonate (ACC) and vaterite at early stages.ACC is the thermodynamically least stable among all calcium carbonate (CaCO 3 ) polymorphs, followed by vaterite, aragonite, and calcite.Under highly supersaturated conditions, ACC nanoparticles with a size of ~ 30 nm are formed within a minute; subsequently, these particles grow and transform with time to more stable crystalline forms such as vaterite crystals (Pouget et al. 2009;Liu et al. 2010), and eventually calcite (Ogino et al. 1987;Rieger et al. 2007).These polymorphs became unstable once the central part of pore body was blocked by precipitates, and then may have dissolved and re-precipitated under fl owing condition.Unlike static conditions, dissolved components can be transported away in the fl owing conditions, which can enhance dissolution rate dramatically.In addition, the surface features (e.g., pits) at nanoscale may enhance the dissolution.It was later computed that under highly supersaturated conditions used in the experiment the saturation ratio of ACC is greater than unity (Boyd et al. 2014), indicating that initial precipitation is likely to be initiated through homogeneous nucleation.
The coupled pore-scale reactive transport model was also used to account for the observed CO 2 seeps at the Grand Wash Fault fi eld site (Crystal Geyser, Utah), a natural analog to geologic CO 2 sequestration (Altman et al. 2014).Field observations along the Grand Wash fault suggest that CaCO 3 precipitation alters fl ow paths by fault zone cementation, resulting in shifts in preferential fl ow paths of the upward migrating CO 2 -saturated-brine.Pore blocking was observed to be faster in a fracture near the main fault and propagated further away from the main fault, which qualitatively refl ects the observed pattern at the Grand Wash Fault.This demonstrates that inclusion of more realistic pore confi gurations and geochemical composition may enhance our fundamental mechanistic explanations of how calcite precipitation alters fl ow paths by pore plugging.

Example results illustrating feedback between fl ow and biofi lms
The pore-scale biofi lm models have been used to investigate several important processes in porous media, including clogging and contaminant transformation, microbially-induced calcium carbonate precipitation (MICP), and bioengineering.We summarize a few of the many interesting works that emphasize coupling among fl uid fl ow, mass transfer, and biofi lm growth.For practical applications, bioclogging has been studied for many years using small porous medium reactors in the laboratory where it is generally only possible to observe bulk changes in permeability and other macroscopic properties (Vandevivere and Baveye 1992;Seifert and Engesgaard 2012).However, most pore-scale modeling studies of bioclogging have been in 2-D.Recently, von der Schulenburg et al. ( 2009) have conducted LB simulations in a 3-D bead packing.The porous medium geometry is derived from images taken by magnetic resonance imaging of saturated packing of random glass beads with mean diameter 400 m and overall porosity 0.59.The biofi lm develops over time, substrate is consumed, and the heterogeneity of the velocity fi eld increases.The LB simulation results show that the computed permeability was smaller than the value computed from a commonly used correlation equation (i.e., the Ergun equation).Comparison of 2-D and 3-D simulations revealed that 3-D simulation is required to accurately model bioclogging.Pintelon et al. (2012) extended the work of von der Schulenburg et al. ( 2009) by including biofi lm detachment and allowing water fl ow through the biofi lm to investigate the impact of biofi lm permeability and strength upon macroscopic properties such as the Da number and permeability.The level-set method of Pintelon et al. (2009) was used to include shear induced detachment.Biofi lm permeability was included by setting the fl uid viscosity within the biofi lm to a high value, similar to the approach used by Thullner and Baveye (2008) in their pore network model.MRT-LBM was used to solve for the fl ow, and the different viscosity values for the fl uid and biofi lm nodes were achieved by adjusting the relaxation rates.A 3-D bead pack porous medium is considered, similar to that of von der Schulenburg et al. (2009).Results from a 2-D slice transverse to the mean fl ow direction are shown in Figure 13.The top row is for a biofi lm with high material strength while the bottom is for a weak biofi lm.The parameter X equals the viscosity ratio of the biofi lm to fl uid cells, so from left to right the biofi lm region is more permeable to water fl ow.The biofi lm with the weaker cohesive strength has signifi cantly less overall accumulation than the stronger biofi lm; biofi lm permeability has a secondary infl uence because fl uid shear stress is smaller for the more permeable biofi lm.The permeability reduces over time as the biomass accumulates.It was observed that the biofi lm permeability needs to be relatively high (low X) before there is a signifi cant impact upon the overall permeability of the system.
Unlike the cases reported above with a single limiting substrate, biofi lms can develop in the situation where separately fl owing electron donors and acceptors meet and mix.This situation can occur during in situ bioremediation along the fringes of a contaminant plume from a persistent source where oxygen present in the ambient groundwater mixes with a dissolved organic pollutant resulting in a "reaction zone" with biomass along the plume fringe (e.g., Cirpka and Valocchi 2007;Valocchi 2012).Werth and co-workers (Zhang et al. 2010b;Yoon et al. 2014) have used 2-D microfl uidic pore networks to investigate these transverse reaction processes experimentally as described in the precipitation section.The electron donor and acceptor were injected through two separate inlets, allowing two substrates to mix around the central mixing zone.Initially, a pure bacterial strain was inoculated and the resulting biofi lm developed along the mixing zone where both substrates are present.Following the experiments in Zhang et al. (2010b) where the biofi lm developed with the aerobic degradation of the herbicide (i.e., R-2,4-DP), Yoon et al. (2014) introduced another substrate to induce adaptation of the strain to degrade a new substance on which an unadapted strain is not able to grow.Over a long term operation (~ 230 days), relatively thick and sturdy biofi lm was developed along the mixing line.Due to the presence of thick and non-uniform biofi lm, transverse mixing increased dramatically, which was also signifi cantly affected by fl owrates.Tang et al. (2013) developed a 2-D pore-scale model that uses LB-FDM for the fl ow fi eld and the advection-diffusion-reaction equations (with dual Monod kinetics), and the CA model for biofi lm spreading.There was a good qualitative comparison between the growth pattern in the model and micromodel experiment, for a close up of the system in the mixing zone (Fig. 14).The biofi lm develops along the downstream side of the cylindrical post and extends out into the pore body.It is important to note that the model biofi lm would completely bridge the pore body unless the CA rules are modifi ed to prevent biomass spreading to zones with fl uid shear above a certain threshold (Knutson et al. 2005).In the model of Tang et al. (2013), the biomass is assumed to be impermeable, and hence the growing biofi lm reduces the velocity in the mixing zone which results in decreased transverse mixing between the donor and acceptor (Fig. 14).The importance of this effect was further evaluated by repeating their simulations using the velocity fi eld computed without any biomass.They then computed the mass of donor degraded from the inlet to the outlet and compared that to what was measured in the experiment.The results clearly show that the feedback among biofi lm growth, velocity alteration, and mass transfer is important.Neglecting the dynamic changes in the pore space and instead using a constant velocity results in an overestimation of mixing and overall reaction.
There is an increasing knowledge about the mechanical properties of biofi lms, which will lead to more complete and accurate models of attachment and detachment processes.Although several groups have incorporated some of these processes into discrete CA or IB models of biofi lms (e.g., Pintelon et al. 2012 ;Bottero et al. 2013), physics-based models using Latt ice Boltz mann-Based Approaches for Pore-Scale Reactive Transport 423 mixture theory (e.g., Lindley et al. 2012) or dissipative particle dynamics (Xu et al. 2011) may be more fl exible.These more complete mechanical models could address the transport and possible re-attachment of biomass that detaches from upstream.

FUTURE RESEARCH DIRECTIONS
We have summarized LBM models for fl uid fl ow and reactive transport processes involving pore structure evolution and fl ow feedback due to mineral precipitation and dissolution and biofi lm dynamics.In particular, recent algorithms for heterogeneous reactions and biofi lm growth at mineral surfaces were thoroughly reviewed.Several key hybrid approaches combining LB with other methods were briefl y introduced.Due to the growing body of literature over diverse disciplines, many ongoing research areas are only briefl y cited and summarized.Despite the signifi cant amount of works, upscaling coupled pore-scale phenomena including heterogeneous biogeochemical reactions, chemical-mechanical-hydrological coupling, and links between pore-and continuum-scale rates are still frontier research questions.
Recent progress in imaging, experimental and computing technologies has resulted in an increased adoption of digital visualization or direct measurement methods in pore-scale reactive transport studies.Using multi-scale imaging methods, we are now able to directly measure how surface chemistry, biogeochemical reactions, and pore topology impact reactive transport processes at the pore scale.Likewise, recent increases in parallel computing and computational methods allow us to simulate reactive transport processes in great detail.For example, it will be possible within the near future to simulate reactive transport experiments directly on the reconstructed 3-D image of porous and fractured media from nano to core scales.Although we focused on reviewing fl ow and reactive transport processes under saturated conditions, unsaturated porous media have also been recently studied (e.g., Raoof 2011;Yang et al. 2014a).Coupled pore-scale fl ow and reactive transport under variably saturated conditions will be an important future research area due to its signifi cance for soil biogeochemical processes in terrestrial systems for carbon cycle and the interaction between hydrological and biogeochemical processes for water and mass fl uxes.
Many of the important engineering challenges that motivate the study of biofi lms in porous media (e.g., bioclogging, microbial enhanced oil recovery, bioremediation) involve coupled biotic and abiotic processes.For example, injection of electron donors for biostimulation to reduce dissolved metals in groundwater at some US DOE facilities can require modeling a complex network of biotic and abiotic reactions (e.g., Scheibe et al. 2006).However, porescale modeling has been limited because there are conceptual and mathematical challenges for modeling simultaneous biofi lm spreading and mineral precipitation.A recent attempt by Zhang and Klapper (2014) using the phase-fi eld model is an encouraging sign of developing a pore-scale model for coupled biotic and abiotic processes.Some researchers have begun to explore so-called in silico models of microbial metabolism based on genomic knowledge instead of the phenomenological Monod model for substrate utilization (e.g., Fang et al. 2011) and to incorporate microbial motility characteristics into cell-scale and porous media models (e.g., Kusy and Ford 2007).Inclusion of more realistic microbial phenomena in pore-scale biofi lm models will signifi cantly improve model predictions of coupled biotic and abiotic processes.
The LB-based approaches and experimental and numerical results reviewed in this article are leading to better understanding of complex problems of evolving pore morphology due to tightly coupled transport and nonlinear reactions.However, the ultimate goal is to improve predictive capability at larger continuum and fi eld scales.As conventional homogenization and volume-averaging techniques are not valid due to process nonlinearity and the lack of scale separation, pore-scale modeling is enabling new methods of upscaling, including multi-scale, multi-physics approaches (Chu et al. 2013;Tomin and Lunati 2013;Varloteaux et al. 2013;Scheibe et al. 2015).A key assumption of the methods reviewed here is that all processes as well as the pore geometry can be quantifi ed unambiguously at the smallest scale (i.e., pore scale).There are new emerging problems (e.g., unconventional resources in shale and carbonate rocks, caprock for CO 2 sequestration) for which it may be necessary to model key processes at the nanometer scale, but the practical scale of observation is generally limited to a micro-scale (10's-100's microns).New approaches based on mesoscopic theories, which connect the microscopic and macroscopic descriptions of the dynamics, provide a promising outcome (DOE 2012).Key relationships and parameter values at the macroscopic scale can be constructed from microscopic-scale simulations including sub-micron scale, which account for a range of reaction and fl ux regimes.For such problems, the LB-based modeling framework reviewed in this article as well as other multiscale simulation methods can help us improve our understanding of reactive transport processes by elucidating knowledge gaps and providing a tool in predicting coupling relevant processes across scales.

Figure 1 .
Figure 1.Examples of various mineral precipitation characteristics.Microphotographs of (a) unaltered and (b) altered sandstone in the vicinity of a natural CO 2 seepage conduit (Little Grand Wash Fault, Utah) [Used by permission of the American Chemical Society, from Altman et al. (2014) Journal of Physical Chemistry C, Vol.118, Fig. 2, p. 15106].Primary porosity in blue epoxy in (a) is altered by carbonate precipitates shown in (b).(c) Uniform elongate concretion from the Sierra Ladrones Formation (ancestral Rio Grande sediments, see Mozley and Davis (2005) for further information).(d) CaCO 3 precipitation along the vertical pathway sealed by a thin Mancos shale layer (Carmel Formation, Utah).(e) Calcite-cemented hanging-wall damage/mixed zone of the Sand Hill Fault, New Mexico (scale: >10 m).

Figure 2 .
Figure 2. Biofi lm growth and spreading models assume that the biofi lm is either represented as a continuum (A) or as discrete quantities (B and C) [Used by permission of IWA Publishing, modifi ed from IWA Task Group on Biofi lm Modeling (2006), Fig. 3.12, p. 106].

Figure 3 .
Figure 3. Schematics of the D2Q9 (left) and D3Q19 (right) models.The lattice velocities are indicated by arrows starting from the square/cube center.Note there is a zero lattice velocity e 0 = 0 in both the models.
32) which is a general formula that can describe all the three types of boundary conditions: the Dirichlet boundary condition, with b 1 = 0 and b 2 ≠ 0, the Neumann boundary condition, with b 2 = 0 and b 1 ≠ 0, and a mixed boundary condition, with b 1 ≠ 0 and b 2 ≠ 0.

Figure 4 .
Figure 4.A schematic illustration of a fl uid-solid interface node in the D2Q5 lattice model.F and S represent a fl uid and solid node, respectively.R represents the reactive interface node.

Figure 7 .
Figure 7. Developing biofi lm using the CA modeling approach on a fl at surface due to diffusion of a single limiting substrate from the top boundary.The system is transport limited [Used by permission of John Wiley and Sons, from Picioreanu C, van Loosdrecht MCM, Heijnen JJ (1998) Mathematical modeling of biofi lm structure with a hybrid differential-discrete cellular automaton approach.Biotechnology and Bioengineering, Vol.58, Fig. 6, p. 111].

Figure 10 .
Figure 10.Evolution of pore structure resulting from CaCO 3 precipitation along a cross-sectional slice perpendicular to the main fl ow direction (light gray for original grain and dark red for precipitated regions).3-D pore structure reconstructed from micro-CT images of a Berea sandstone was used [Used by permission of American Physical Society, from Jiang F, Tsuji T (2014) Changes in pore geometry and relative permeability caused by carbonate precipitation in porous media.Physical Review E, Vol.90, Fig. 8, p. 053306-6].
Transverse mixing-induced mineral precipitation experiments were performed by injecting solutions containing reactants (e.g., Ca 2+ , Mg 2+ , UO 2 separate inlets of a micromodel.The amount of mineral precipitation, snapshots of the morphology of precipitates, elemental analysis with Raman spectroscopy, and LSCM imaging and tracer tests

Figure 11 .
Figure 11.(a) A snapshot of calcium carbonate precipitation with 25mM Ca 2+ (upper inlet) and 25mM CO 3 2-(bottom inlet).Precipitation patterns at three different locations are highlighted [Used by permission of John Wiley and Sons, from Yoon H, Valocchi AJ, Werth CJ, Dewers T (2012) Pore-scale simulation of mixing-induced calcium carbonate precipitation and dissolution in a microfl uidic pore network.Water Resources Research, Vol.48, Fig. 1, p. W02524-3].(b) Tracer test result after pore blockage.A profi le of 3-D precipitates halfway down the length of the micromodel was obtained using confocal microscope, showing a complete pore blockage over a ~40 μm depth [Used by permission of John Wiley and Sons, from Fanizza MF, Yoon H, Zhang C, Oostrom M, Wietsma TW, Hess NJ, Bowden ME, Strathmann TJ, Finneran KT, Werth CJ (2013) Pore-scale evaluation of uranyl phosphate precipitation in a model groundwater system.Water Resources Research, Vol.49, Figures 10-11, p. 886].(c) Microscopy images of CaCO 3 polymorphs in the presence of 40 mM Mg 2+ with 10mM Ca 2+ .(upperinlet) and 10mM CO 3 2-(bottom inlet).The central mixing line in the micrometer is shown in an upper-left pore body [Used by permission of Elsevier, from Boyd T, Yoon H, Zhang C, Dehoff K, Fouke B, Valocchi AJ, Werth CJ (2014) Infl uence of Mg 2+ on CaCO 3 precipitation during subsurface reactive transport in a homogeneous siliconetched pore network.Geochimica et Cosmochimica Acta, Vol.135, Fig. 5, p. 329].

Figure 13 .
Figure 13.2-D slice from the 3-D model results ofPintelon et al. (2012).Biomass is shown in green along with the shear stress fi eld after 27 h of growth.Left panel is for an impermeable biofi lm, and the biofi lm permeability increases from left to right.Biomass was allowed to detach depending on adjacent shear stresses; the top row is a comparatively strong biofi lm, while the bottom row is a comparatively weak biofi lm.As the parameter X decreases, the biofi lm permeability increases.[Used by permission of John Wiley and Sons, from Fredrich JT, DiGiovanni AA, Noble DR (2006) Predicting macroscopic transport properties using microscopic image data.Journal of Geophysical Research, Vol.111, Fig. 7, B03201-10].

Figure 14 .
Figure 14.Comparison of the biofi lm spreading in the experiment (top) of Zhang et al. (2010b) with the model results (bottom) by Tang et al. (2013) after 11 days (left) and 24 days (right).[Top: Used by permission of American Chemical Society, from Zhang C, Kang Q, Wang X, Zilles JL, Mü ller RH, Werth CJ (2010b) Effects of pore-scale heterogeneity and transverse mixing on bacterial growth in porous media.Environmental Science & Technology, Vol.44, Fig. 3, p. 3088; Bottom: Used by permission of John Wiley and Sons, from Tang Y, Valocchi AJ, Werth CJ, Liu H (2013) An improved pore-scale biofi lm model and comparison with a microfl uidic fl ow cell experiment.Water Resources Research, Vol.111, Fig. 6, p. 8376].