Micro-Continuum Approaches for Modeling Pore-Scale Geochemical Processes

The recent profusion of microscopic characterization methods applicable to Earth Science materials, many of which are described in this volume (Anovitz and Cole 2015, this volume; Noiriel 2015, this volume), suggests that we now have an unprecedented new ability to consider geochemical processes at the pore scale. These new capabilities offer the potential for a paradigm shift in the Earth Sciences that will allow us to understand and ultimately quantify such enigmas as the apparent discrepancy between laboratory and field rates (White and Brantley 2003) and the impact of geochemical reactions on the transport properties of subsurface materials (Steefel and Lasaga 1990, 1994; Steefel and Lichtner 1994; Xie et al. 2015). It has only gradually become apparent that many geochemical investigations of Earth materials have suffered (perhaps inadvertently) from the assumption of bulk or continuum behavior, leading to volume averaging of properties and processes that really need to be considered at the individual grain or pore scale. For example, a relationship between reaction-induced porosity and permeability change can perhaps be developed based on bulk samples, but ultimately a mechanistic understanding and robust predictive capability of the associated geochemical and physical processes will require a pore-scale view.

The question still arises: Do we need pore-scale characterization and models in geochemistry and mineralogy? The laboratory–field rate discrepancy (or enigma) is a good example of where a pore-scale understanding may provide insights not easily achievable with bulk characterization and models. If the reasons for this apparent discrepancy between laboratory and field rates cannot be explained, then it appears unlikely that scientifically defensible and quantitative models for a number of important Earth Science applications ranging from chemical weathering and its effects on atmospheric CO2, to subsurface carbon sequestration, to nuclear waste storage, to contaminant remediation and transport, …


INTRODUCTION
The recent profusion of microscopic characterization methods applicable to Earth Science materials, many of which are described in this volume (Anovitz and Cole 2015, this volume;Noiriel 2015, this volume), suggests that we now have an unprecedented new ability to consider geochemical processes at the pore scale.These new capabilities offer the potential for a paradigm shift in the Earth Sciences that will allow us to understand and ultimately quantify such enigmas as the apparent discrepancy between laboratory and fi eld rates (White and Brantley 2003) and the impact of geochemical reactions on the transport properties of subsurface materials (Steefel andLasaga 1990, 1994;Steefel and Lichtner 1994;Xie et al. 2015).It has only gradually become apparent that many geochemical investigations of Earth materials have suffered (perhaps inadvertently) from the assumption of bulk or continuum behavior, leading to volume averaging of properties and processes that really need to be considered at the individual grain or pore scale.For example, a relationship between reactioninduced porosity and permeability change can perhaps be developed based on bulk samples, but ultimately a mechanistic understanding and robust predictive capability of the associated geochemical and physical processes will require a pore-scale view.
The question still arises: Do we need pore-scale characterization and models in geochemistry and mineralogy?The laboratory-fi eld rate discrepancy (or enigma) is a good example of where a pore-scale understanding may provide insights not easily achievable with bulk characterization and models.If the reasons for this apparent discrepancy between laboratory and fi eld rates cannot be explained, then it appears unlikely that scientifi cally defensible and quantitative models for a number of important Earth Science applications ranging from chemical weathering and its effects on atmospheric CO 2 , to subsurface carbon sequestration, to nuclear waste storage, to contaminant remediation and transport, can be fully developed and applied.The reasons for the discrepancy (apparent or real) have been widely discussed, and over time the number of possibilities for explaining it have narrowed.One potentially important effect that contributes to this apparent laboratory-fi eld rate discrepancy is geochemical in origin and has to do with the fact that most laboratory studies do not consider mineral dissolution as regulated by the precipitation of a secondary phase, that is, as an incongruent reaction.As proposed by Zhu and co-workers and as investigated further by Maher and co-workers, the slow precipitation of secondary clay minerals as a result of primary silicate system conditions (Navarre-Sitchler et al. 2011).Or even where the pores are connected, the reactive phases within the pore are coated with other (typically secondary) phases.The fl uids in the connected pores then might "see" only the secondary phases (e.g., clay or iron hydroxide) rather than the reactive feldspar.In this respect, both the impacts of the porosity and the reactive mineral surface area cannot be completely separated from the transport properties of the medium under consideration.Similarly, the physical surface area of reactive minerals is another parameter that may be relatively straightforward to determine on bulk samples, whether using a grain size geometric analysis or using gas adsorption isotherms, e.g., BET methods (Brunauer et al. 1938).But the physical surface area of the minerals may not translate to a unique reactivity without consideration of the reactive site density, which can vary signifi cantly between natural samples (Beig and Lüttge 2006;Bracco et al. 2013).
Two-dimensional and three-dimensional images can be analyzed to extract reactive transport model parameters.This includes the sample specifi c parameters of porosity, mineral volume fractions and surface areas, diffusivity and even permeability.Details regarding the collection and underlying principles of various imaging methods, including 2-D Scanning Electron Microscopy (SEM), and 3-D Focused Ion Beam-Scanning Electron Microscopy (FIB-SEM), X-ray Computed Micro-Tomography (X-ray microCT), optical petrology and Small Angle Neutron Scattering (SANS) of nanoscale porosity have been discussed in other articles (Anovitz and Cole 2015, this volume;Noiriel 2015, this volume) and thus will only be mentioned briefl y here.
In this article, we will discuss the data processing approaches that can be used to extract model parameters from two-and three-dimensional SEM, FIB-SEM, and X-ray microCT images.In addition, some of the outstanding issues, such as image segmentation and resolution, will be discussed in the context of their effect on parameter estimation.Image segmentation refers to the partitioning of the histogram of pixel intensities (SEM imaging) or voxel attenuation values (X-ray CT imaging) into one or more categories such as pores and grains.Comparison of parameters from bulk samples (e.g., reactivity) and from image-derived micro-continuum samples, however, should provide insight into the scaling relationships in reactive porous media.Finally, it should be noted that several of these parameters can be measured on bulk samples in the laboratory.Parameter estimation from images allows for more discrete parameter evaluation, including the ability to map parameters at multiple scales or associated model grid cell sizes.

Porosity
Sample porosity can be easily determined from imaging techniques.Given that porosity is an intensive property, it can be computed from either 2-D or 3-D images.Perhaps the simplest approach is determining porosity from 2-D SEM images of a polished section.Polished sections, including thin sections, can be easily prepared by impregnating samples with epoxy, curing, and then cutting and polishing the samples to the desired thickness.Numerous commercial companies also offer inexpensive sectioning services, removing the limitations of experience or facilities.Further guidance and extensive details on sample preparation for SEM imaging can be found in several existing texts (Echlin 2011).
Scanning electron microscopes, now widely available, use an electron beam to capture electron-sample interactions.Most SEM instruments are also equipped with a backscattered electron (BSE) detector.In this mode, the degree of backscatter is proportional to the mean atomic number, producing an image with varying grayscale intensities (Krinsley et al. 2005).Before imaging, polished samples are typically coated with a thin layer of conductive material such as carbon or gold using an ion beam sputterer to prevent surface charging when imaging.This is unnecessary if the instrument is operated in environmental mode or under low vacuum.Pores are easily distinguishable in SEM BSE images of most geologic samples in which the pore space shows a signifi cant contrast with the minerals.

Steefel, Beckingham & Landrot
Benchtop or synchrotron CT images can provide a three-dimensional depiction of a coherent geologic sample.In this method, a series of radiographs or projections are collected over a range of angles (Cnudde and Boone 2013;Wildenschild and Sheppard 2013) and a 3-D image.Several variations in image acquisition as well as complexities in image reconstruction exist and are covered in more detail in existing reviews (Wildenschild and Sheppard 2013;Noiriel 2015, this volume).In X-ray CT imaging, voxel attenuation is proportional to the energy of incident X-rays, material density, and atomic number (Cnudde and Boone 2013;Wildenschild and Sheppard 2013).The contrast in attenuation allows pores to be distinguished from grains, although varying degrees of phase retrieval may be required to make the approach fully quantitative (Wildenschild and Sheppard 2013).
Micron-scale resolution is possible with both SEM and microCT imaging, but this resolution is not suffi cient to capture the abundant sub-micron-scale pores in carbonates and shales.FIB-SEM imaging, however, can be used to characterize nanoscale porosity (Curtis et al. 2012;Landrot et al. 2012).Using FIB-SEM, high-resolution (nm) 3-D images are developed from a series of 2-D SEM images reconstructed to a 3-D volume.Before imaging, a trench is milled in front of the area of interest using the ion beam.A SEM image is then captured before the ion beam is used to mill away a small layer from the sample.This sequence of image capturing and milling is cycled, generating a series of 2-D images.Given the close spacing of the slices, these 2-D images can be reconstructed into a 3-D volume (Curtis et al. 2012;Landrot et al. 2012).Small angle neutron scattering or SANS can also be used to investigate nanoscale pore distributions and processes, but it is a statistical rather than mapping technique and is not discussed further here, although interested readers can fi nd discussion in other articles (Anovitz and Cole 2015, this volume;Navarre-Sitchler et al. 2015, this volume).
The porosity in SEM, X-ray CT, and FIB-SEM images can be determined by computing the ratio of pore pixels or voxels to total pixels or voxels in the 2-D or 3-D image.This fi rst requires segmentation of pore and grain pixels/voxels.There are a variety of existing segmentation techniques that have been used with varying success in the literature.In general, the pore-grain threshold occurs at a minimum in the histogram of grayscale intensities between the individual intensity distributions for pores and grains (Peters 2009).The choice of thresholding technique should be carefully made so as to be optimal for the sample of interest, as further discussed below.In addition, some samples will require extensive pre-segmentation fi ltering and even manual correction to remove image artifacts.Once the segmented image is produced, pore and grain pixels can be easily summed using commercially available image processing software or using individually developed computer programs.
While porosities can be determined from either 2-D or 3-D images, a suffi cient number of images are required in order to ensure that the volume used is representative of the sample, thus obtaining reliable porosity, mineral volume fractions, and mineral surface areas.The representativeness of the area or volume can be determined by computing the porosity on smaller volumes subsampled from the original image.As the sampled volume is increased, the computed porosity should approach a uniform value as a representative elementary volume (REV) is reached.

Mineral volumes
Despite recent efforts based on 3-D imaging (Mutina and Koroteev 2012), mineral volume fractions in mineralogically complex systems can be reliably determined only from 2-D SEM imaging.Where only a single mineral is present, as in the study by Noiriel et al. (2012), or where there is a signifi cant contrast in density between the minerals present, 3-D X-ray synchrotron mapping may be able to provide quantitative determinations of mineral volumes.Mineral volume determination at the microscopic scale is possible on SEMs equipped with Energy-Dispersive X-ray Spectroscopy (EDS or EDX) capabilities.The 2-D BSE imaging approach described above makes it possible to distinguish between quartz, clay, and other reactive minerals (Peters 2009) with a few hundreds of nm per pixel length resolution.EDX imaging further allows for mineral identifi cation by determining the elements present at the microscopic scale.Two-D elemental maps can be captured and processed to determine the mineral volume percentages within each pixel (Landrot et al. 2012).Processing of EDX signals can be carried out with commercial software, such as QEMSCAN, or though customized thresholding and processing codes such as those presented in Landrot et al. (2012).These methods couple information on BSE intensities with elemental intensities so as to identify individual minerals.In either method, knowledge of the bulk mineralogy is needed to aid in mineral characterization, as many minerals contain similar elemental compositions and BSE intensities.This can be obtained from X-ray Diffraction (XRD) or X-ray Fluorescence (XRF) on bulk samples.Identifi cation of minerals can be challenging in highly heterogeneous samples even when commercial software like QEMSCAN is used.In addition, grain edge artifacts at grain-epoxy boundaries can alter BSE intensities (Dilks and Graham 1985), potentially causing misidentifi cation of minerals at grain boundaries.
Following mineral segmentation, computation of mineral volume fractions can be carried out by counting the number of pixels corresponding to each mineral and dividing by the total grain pixels.An additional advantage of determining mineral volumes from imaging is the ability to assess explicitly the accessibility of minerals to reactive fl uids if the porosity is imaged as well, and if the imaging is at suffi ciently high resolution that the connectivity of the pore network can be quantifi ed.As observed in recent studies (Peters 2009;Landrot et al. 2012), mineral abundance alone may not always accurately refl ect mineral accessibility and pore connectivity.At least in continuum models, it is typically more accurate to discretize minerals based on pore accessibility fractions rather than based on total volume fractions.Mineral accessibilities can be computed by fi rst identifying the grain-pore boundary and then by counting the number of associated pixels for each mineral that are present at the interface, or adjacent to the pore space.Recent studies have also considered the connectivity of the pore space as well, since this is a primary control on the transport of ions to and from reactive surfaces (Navarre-Sitchler et al. 2009;Landrot et al. 2012).

Mineral surface area
There are several complexities and intricacies involved in obtaining and interpreting mineral surface area, a full review of which is beyond the scope of this chapter.Instead, we will briefl y give a few examples of ways to interpret surface area from 2-D and 3-D images.It should be noted that it is uncertain what the appropriate mineral surface areas are for reactive transport modeling of porous media.A range of surface area estimates have been used without evaluation of their impact or success.This includes liberally interchanging geometric surface area (GSA), specifi c surface area (SSA), and reactive surface area (RSA).Mineral GSA typically refers to a surface area computed from an average grain size and assuming a particular geometric grain shape, typically perfectly smooth spheres (White et al. 2005;Alemu et al. 2011).Mineral SSA refers to the total or rough surface area per gram mineral, often measured via the Brunauer, Emmett, Teller (BET) analysis (Brunauer et al. 1938).RSA accounts (or attempts to account) for the distribution of reactive sites on a mineral surface and is usually estimated by applying a scaling factor of one to three orders of magnitude to SSA or GSA (Zerai et al. 2006), although the basis for applying this factor is not clear.Geometric, specifi c, and reactive mineral surface areas can be approximated from 2-D and 3-D images using a variety of approaches.Geometric surface areas can be approximated from the grain sizes observed in 2-D or 3-D images assuming typical grain geometry and either an average grain size or range of grain sizes and corresponding surface areas.An alternative approach is to assume mineral-specifi c geometries and use image-observed grain dimensions to compute surface areas based on the ideal geometries for each mineral (Bolourinejad et al. 2014).This approach is somewhat limited as real minerals often deviate from ideal geometries or do not have perfectly smooth surfaces (Bolourinejad et al. 2014).It should also be noted that these 2-D approaches additionally require bias correction (Weibel 1979;Crandell et al. 2012).Specifi c surface areas that refl ect surface roughness can be estimated from GSA determinations by applying a roughness factor, although the choice of a factor is diffi cult to justify and only an average at best (Zerai et al. 2006).Alternatively, the BET surface area can be measured in the laboratory and distributed among minerals as identifi ed in 2-D images.In 3-D images, specifi c surface areas can be successfully computed by creating a triangulated or polygonized surface mesh, such as through the marching cubes algorithm (Lorensen and Cline 1987;Landrot et al. 2012), then summing the area of the polygons.This is necessary as it has been shown that simply counting surface voxel faces in 3-D images overestimates surface area, as is explained in the review by Wildenschild and Sheppard (2013).This approach can be used to compute surface areas from both FIB-SEM images and X-ray CT images.

Diffusivity
Diffusivity from bench-top experiments.Diffusivity within homogeneous porous samples may be determined experimentally with a range of techniques, some of which do not require pore-scale imaging.The most common experimental techniques are based on a diffusion cell in which one end of the cell is held at a constant concentration (effectively a Dirichlet or fi xed boundary condition), while a reservoir on the other side of the diffusion cell is monitored for solute breakthrough, either under transient (e.g., where the reservoir is stagnant) or steady-state conditions (where the reservoir is subject to fl ow).For example, Figure 1 shows a schematic for an experimental setup that has been used to determine ion diffusivity in bentonite clay (Tachi and Yotsuji 2014).
Another possible experimental approach for porous medium samples relies on chemical mapping of the diffusion profi le.For example, Navarre-Sitchler et al. ( 2009) used micro-X-ray synchrotron mapping of bromide fl uorescence (XRF) to determine the diffusivity of samples of unweathered and weathered basaltic andesite.An effective, upscaled diffusion coeffi cient was determined for the porous andesite by fi tting the profi le with a numerical solution of the diffusion equation where * i D is the diffusion coeffi cient for species i in porous media that incorporates the effects of tortuosity, C i is the tracer concentration,  is the divergence operator,  is the porosity, and t is time (Steefel et al. 2015).For example, a XRF map of a sample of weathered basalt from a Costa Rican chronosequence shows the distribution of bromide tracer (blue) after 7 days resulting from diffusion from left to right, as shown in Figure 2  Diffusivity from numerical experiments.Experiments and/or diffusion profi le mapping work well for determining diffusivity where the properties of the porous material are largely homogeneous, since in this case a single parameter value can be fi tted to the data.In some cases, it may also be possible to estimate diffusivities for heterogeneous samples if the distribution of properties (e.g., grain size) is known, although the likelihood of a non-unique fi t increases with the number of different properties in a sample.Alternatively, if the distribution of properties is unknown, it is possible to estimate diffusivity even in highly heterogeneous samples using numerical modeling based on 2-D or 3-D pore structure characterization.Navarre-Sitchler et al. ( 2009) made use of X-ray synchrotron microtomography to map the pore structure of samples of weathered basalt similar to that shown in Figure 2. Using a simple implementation of thresholding to map basalt versus pores, they were able to delineate chemical weathering related macroporous zones (> 4.4 m voxel resolution) that were connected in 3-D.Navarre-Sitchler et al. (2009) then carried out numerical tracer diffusion experiments in 3-D cubes of weathered basalt by assuming a low diffusivity of 1.75 × 10 -14 m 2 s -1 for the largely unconnected pore structure of unaltered basalt and a free ion diffusivity of 10 -9 m 2 s -1 (corresponding to a tortuosity of 1.0) for connected pores that can be fully resolved with the 4.4 m discretization.Implemented in this way, the numerical tracer diffusion simulations with the code CrunchFlow are similar to what could be done with a pore network model.A 2-D slice through the skeletonized pore structure of one of the weathered zones is shown in Figure 3A   used and they are based on the segmented 3-D porosity map of the weathered sample.The upscaled effective diffusion coeffi cient, in contrast, is determined by fi tting a 1-D profi le to the 3-D simulation results.The diffusion of the bromide tracer was assumed to follow a modifi ed Archie's Law model that incorporates a critical porosity threshold value of 9%, below which the porosity is considered to be largely unconnected: w here D p is the diffusivity in the unweathered parent basalt, D 0 is the diffusion coeffi cient in pure water, and D e is the effective diffusion coeffi cient in porous media and where the effective porosity,  e , is defi ned as: The parameter a in Equation ( 3) is taken as 1.3 while the parameter  is assumed to be 1.0 for 3-D volumes measuring 220 m on a side (Navarre-Sitchler et al. 2009), while  T refers to the total porosity as mapped with the X-ray synchrotron microtomography with a 4.4 m voxel resolution, and  c is the critical porosity estimated as 9% based on the same data.
A summary of the 3-D tracer diffusion simulation results for the Costa Rican basalts are given in Table 1.Samples with less than 10% porosity are essentially unweathered.Note that the results are to some extent dependent on sample size, and also on the porosity of the volume considered as expected from Equation (2).The results of the tracer diffusion simulations agree broadly with the laboratory tracer diffusion experiments using bromide (Fig. 4).
Using numerical modeling based on pore-scale imaging for the purposes of estimating diffusivity offers practical advantages when the material of interest within a sample is too small to easily investigate with laboratory tracer experiments.As an example, consider the chloritefi lled pores within the Lower Tuscaloosa Formation that hosts the Cranfi eld CO   0.17 3.7 × 10 -15 2 0.17 1.7 × 10 -12 10 0.17 6.1 × 10 -12 14 0.17 1.9 × 10 -11 21 0.17 4.8 × 10 -11 30 0.17 the nano-crystalline chlorite as reconstructed in 3-D with the techniques described above.
Carrying out a laboratory tracer experiment on such a small zone would be quite diffi cult, but determining an approximate upscaled diffusion coeffi cient is relatively simple using a true porescale model (Molins 2015, this volume), or the approach as described in Figures 3 and 4. As in the Costa Rican basalts described above, the Tuscaloosa Formation sandstone is divided into pores and chlorite mineral and these are represented directly in the pore-scale model using the same grid resolution as the image resolution (14 nm).A non-reactive tracer with a concentration of 0.01 M is released from the left hand side of the cube (Dirichlet boundary) and allowed to diffuse to the right.All other boundaries are treated as no-fl ux.Based on the average diffusion profi le, one can estimate an upscaled, effective diffusion coeffi cient of 7 × 10 -12 m 2 /s (Fig. 6).
The upscaled diffusion coeffi cient can then be used to describe diffusivity in the chlorite-rich zones within micro-continuum representations, as described more fully below.

Permeability
Imaging methods have been increasingly used as the basis for the prediction of permeability (Caubit et al. 2009;Algive et al. 2012;Beckingham et al. 2013).In addition to empirical relationships such as the Kozeny-Carmen equations that compute permeability from experimental or image-computed porosity (Kozeny et al. 1927;Carman 1939), permeability can also be estimated from pore networks extracted from 2-D and 3-D images.Pore and porethroat size distributions as well as connectivities are needed to recreate a representative pore network.A suite of different approaches have been used to characterize pore and pore-throat size distributions from 2-D and 3-D images.This includes, for example, multiple point statistics (Okabe and Blunt 2004), image erosion-dilation (Crandell et al. 2012), maximum inscribed spheres (Baldwin et al. 1996), and watershed segmentation (Beucher and Lantuéjoul 1979;Silin and Patzek 2003).It should be noted that some of these 2-D methods require bias correction (Crandell et al. 2012) and may not be able to determine pore connectivities without relying on information from 3-D images (Beckingham et al. 2013), or determined by some other means.
From these statistical distributions, simple pore network models can be created and used to compute continuum-scale permeability.In these models, a series of pores are defi ned on a regular cubic lattice and connected by pore throats.These statistical distributions of pore sizes also provide information on pore connectivity for the models.Fixed fl uid pressures are then applied at the inlet and outlet (Beckingham et al. 2013).Using Poiseuille's law to describe pore throat conductance, and assuming an incompressible fl uid, the pressure in each node can be determined (Li et al. 2006).At the continuum scale, Darcy's law can then be applied to solve for the permeability (Li et al. 2006).While this method has been shown to successfully predict permeability (Beckingham et al. 2013), image resolution and segmentation effects can affect predictions, as discussed further below.
Estimating permeability using pore-scale numerical simulation.As an alternative to these statistical methods, pore network structures have also been directly extracted and simulated from 3-D images.Direct 3-D network extraction from 3-D images has been carried out using a variety of techniques, such as skeletonization based on the medial axis transform as in seminal work by Lindquist and co-workers (Lindquist et al. 1996).Direct modeling using Lattice Boltzmann based on the image-derived pore structure is often used to estimate fl ow and network properties (Blunt et al. 2013).The 3-D imaging could be X-ray synchrotron microtomography, or for potentially higher resolution, FIB-SEM techniques.For example, Oostrom and coworkers (Oostrom et al. 2014) determined permeability for micro-models based on the pressure drop across the model for a range of imposed fl ow rates using Darcy's Law.While the micromodels investigated by Oostrom et al. (2014) were based on idealized geometries, the approach can be generalized to more complex natural pore structures, for example a capillary tube fi lled with crushed calcite grains (Molins et al. 2014), or in 3-D to fractured shales imaged with high resolution FIB-SEM techniques (Trebotich and Graves 2015).In Molins et al. (2014) and in Trebotich and Graves (2015), the pore structure is fully resolved and the Navier-Stokes equation is solved for the domain, in the case of the fractured shale (Fig. 7) at a resolution of 48 nm.The approach in order to determine an upscaled permeability in the case of multidimensional heterogeneous samples requires an averaging technique for the pressure, since there is no single pressure value at the downstream side of the volume.

Imaging issues impacting parameter estimation
Image segmentation.To date, a variety of segmentation and interpretation methods have been applied to 2-D and 3-D images with varying success (Sezgin and Sankur 2004;Dong and Blunt 2009;Iassonov et al. 2009;Peters 2009;Porter and Wildenschild 2010;Bhattad et al. 2011;Wildenschild and Sheppard 2013).Despite the range of automated thresholding methods, there remains a lack of consistency in results even on images of the same sample (Wildenschild and Sheppard 2013).Typically, operator input remains necessary (Cnudde and Boone 2013).
Several image artifacts and issues continue to make image segmentation challenging.One complexity that is particularly common in 3-D images is partial volume or edge effects The surface is colored according to velocity magnitude in the pore, with blue representing lower velocities, shades of green, yellow, and red representing higher velocities.Resolution 48.4 nm with 1920 × 1600 × 640 total grid cells.Original pore structure obtained from FIB-SEM mapping, which indicates a bulk porosity of 18%.[Reproduced from Trebotich D, Graves DT (2015) An adaptive fi nite volume method for the incompressible Navier-Stokes equations in complex geometries Communications in Applied Mathematics and Computational Science, Vol.10(1), p. 43-82, Fig. 13, with permission from Mathematical Sciences Publishers.] that occur, for example, when a single voxel contains multiple substances or straddles a pore-mineral boundary.In these cases, the corresponding composite voxel is assigned an intermediate X-ray attenuation value that refl ects some combinations of the neighboring voxel compositions (Ketcham and Carlson 2001).Additional diffi culties result from image artifacts, such as noise, surface charging, and ring artifacts.Also, variations in intensities across images are common and occur in SEM images when surfaces are not completely fl at, or as a result of beam hardening in X-ray CT imaging (Wildenschild and Sheppard 2013).Filtering of some of these artifacts can be achieved through relatively simple approaches involving noise removal by pixel fl ipping of isolated misidentifi ed phases Peters (2009).Others artifacts are not easily correctable and require additional more complex fi ltering processes (Blunt et al. 2013;Cnudde and Boone 2013;Wildenschild and Sheppard 2013) or even manual image correction (Crandell et al. 2012) before or after thresholding to correctly segment images.Given the range and extent of segmentation methods and pre-processing procedures, a detailed review is not included here, but the interested reader is referred to other reviews on the subject (Cnudde and Boone 2013;Wildenschild and Sheppard 2013).
Imprecision and errors in image segmentation can have several impacts on parameter estimation.Small-scale features may be misinterpreted as noise and thus erroneously removed pre-or post-thresholding, altering the total porosity, mineral volume fractions, and mineral surface areas.This may also impact the estimation of permeability if these features are removed in pore throats and result in erroneously increasing pore throat size(s).A recent review in Iassonov et al. (2009) evaluated porosities determined from a range of segmentation methods applied to 3-D X-ray CT images of porous media.They found large discrepancies of one to over two orders of magnitude in porosity resulting from different segmentation methods, even on simple samples such as glass beads (Iassonov et al. 2009).Given that this initial segmentation step is needed to defi ne the grain-pore boundaries, segmentation error directly impacts analysis of pore and pore-throat sizes, and thus the permeability that is predicted.Similar pore-grain segmentation methods have been used to process BSE and EDS images to identify minerals (Peters 2009;Landrot et al. 2012).It is thus likely that similar potentially signifi cant errors may result in the case of mineral volumes and surface areas as well.
Care should be taken to reduce segmentation error whenever possible.This can be achieved by evaluating the appropriate segmentation method for each sample.In some cases, different images from the same sample may require different segmentation procedures (Wildenschild and Sheppard 2013).Additionally, image determined parameters should be compared with measured parameters whenever possible.Successful segmentation should produce porosity estimates that agree with laboratory-measured values, taking into account the fraction of porosity that is accounted for in each method.Similarly, image determined mineral volume fractions should agree with laboratory-measured mineralogy from XRD or XRF, keeping in mind that the sensitivity of the instruments may not be suffi cient to characterize the minor phases that can be captured in high-resolution imaging.

Image resolution
High resolution images typically require small sample sizes.There are several methods that have been used with some success to increase the sample size of high-resolution images.In one approach that is often used with 2-D SEM images, a series of overlapping highresolution images are captured and stitched together into a larger, high-resolution, composite image (Crandell et al. 2012).Similar mosaic techniques have been used with 3-D images as well (Mokso et al. 2012).In addition, alternative reconstruction and transform approaches have been successful in some cases in increasing the sample size of high resolution images (Defrise et al. 2006;Cnudde and Boone 2013).The use of multi-scale imaging that relies on a range of imaging techniques or resolutions has recently been gaining interest as well.These include, for example, combining X-ray CT images at different resolutions or fusing X-ray CT images with FIB-SEM and SEM BSE images (Sok et al. 2010;Landrot et al. 2012).These studies used high-resolution FIB-SEM imaging to evaluate the porosity and surface area of clay minerals that could not be characterized with the lower resolution pore-scale imaging approaches (Sok et al. 2010;Landrot et al. 2012).
Regardless of the approach, features smaller than the voxel or pixel size cannot be distinguished with the approaches discussed here.These features, however, may be present via the partial volume effect (Cnudde and Boone 2013).In CT imaging for example, this occurs where a single voxel contains more than one phase, with the result that an intermediate attenuation value is assigned to the voxel.This results in image imprecision, particularly at interfaces between materials (Cnudde and Boone 2013), and can introduce potentially large errors when small features and pores are characterized.Imprecision in small-scale features may have little impact on overall porosity and mineral volume fractions, but large impacts on other image-estimated parameters.For example, there could be a large effect on the determination of mineral surface areas given that small-scale features often account for a large portion of the surface area.This is in addition to the impact that lower image resolutions have on correctly capturing surface roughness and thus surface area.Misclassifying smallscale features in either 2-D or 3-D images also introduces error in permeability and diffusivity predictions when using image-informed network models (Caubit et al. 2009;Beckingham et al. 2013).Accurate permeability predictions may still be possible from models informed by lower resolution images in the case where smaller pores (below the resolution of the technique) do not have a large impact on permeability and fl ow (e.g., Blunt et al. 2013).

MICRO-CONTINUUM MODELING APPROACHES
If the imaging methods are suffi ciently high resolution that individual grains and pores and most importantly their interfaces can be resolved, then the 2-D and 3-D maps can be used directly in high resolution pore-scale modeling (e.g., Molins et al. 2014).However, this requires access to pore-scale reactive transport modeling software (still a specialized fi eld) and to high performance computer hardware.Thus it may not be practical for all researchers approach, even if the full pore-scale approach is arguably the most rigorous.If the domain size is too large and/or spatial resolution too low, or if a true pore-scale code is not available, then micro-continuum modeling is another possibility.The data can be used at the spatial resolution at which is collected, or volume averaged to a coarser discretization so as to handle a larger domain, or simply to make the simulations computationally feasible.Scalar quantities like the porosity and mineral volume fractions can be volume averaged directly, vectorial quantities like permeability or diffusivity are more likely scale-dependent and require special treatment (see discussion below).

Volume averaging of porosity and mineral volume fractions
Given an initial map of porosity and mineral volume fractions at some resolution, it is a relatively simple procedure to calculate equivalent quantities at coarser resolutions by using volume averaging.However, these scalar quantities are infl uenced by the connectivity and thus the transport properties of the medium.If the connectivity and the accessibility of reactive surface area are scale-dependent, then volume averaging as an upscaling procedure may introduce errors into the simulations.How accurate the volume-averaged quantities are at various scales is a future research area, one that has been neglected to date perhaps because researchers interested in upscaling in porous media have focused primarily on the physical properties, or because the routine use of high-resolution chemical and mineralogical mapping is still in its infancy.
If pores and individual mineral grains are fully resolved (i.e., no voxel includes more than a single mineral, or a mixture of pore space and minerals), then producing coarser representations of the medium consists of adding up the various image voxels corresponding to porosity and the individual minerals.The percentage of porosity per unit volume porous medium, for example, is just the number of voxels made up completely of pores divided by the total number of voxels in the volume of interest.The volume procedure, however, is equally straightforward where individual voxels consist of a mix of either different minerals or minerals and porosity, since the porosity (connected or unconnected) will be just the weighted average of the porosity in the individual voxels.An example of volume averaging is given by choosing an image of the porosity and mineral distribution from the Lower Tuscaloosa Formation at the Cranfi eld CO 2 sequestration site with a resolution of 331 nm (Fig. 8), which is adequate to resolve all but the nano-crystalline chlorite-fi lled pores (Landrot et al. 2012).With suffi cient computational resources (software and hardware), it should be possible to simulate the porescale geochemical processes at the original resolution of 4 m, that is, use the information on porosity (accessible and inaccessible) and mineralogy shown in Figure 8 directly.In order to make the reactive transport simulations tractable for our purposes, however, we assume a 256 by 256 2-D section with 16-m grid resolution that produces a section measuring 4.1 mm by 4.1 mm.As an example, volume averaging of porosity and mineral abundance produces data like that shown in Table 2, which represents a small portion of the Lower Tuscaloosa Formation sandstone sample investigated in Landrot et al. (2012).The volume averaging to 16 m 2 produces a porosity map (Fig. 9A) and mineral abundances for quartz, chlorite (chamosite), and illite as shown in Figures 9B, 9C, and 9D respectively.These maps are used as initial conditions in the reactive transport modeling of the 4.1 mm by 4.1 mm 2-D section discussed below.

Micro-continuum reactive transport simulations of fractured tuff.
The fi rst microcontinuum reactive transport modeling described in the literature that we are aware of was presented by Glassley and co-workers (Glassley et al. 2002).This mysteriously unrecognized and largely uncited publication was the fi rst to make use of mineralogical and chemical data collected with modern synchrotron techniques at the micro-scale and then used as initial conditions for high resolution reactive transport simulations.This study developed spatially distributed representations of porosity and mineralogy based on a combination of optical mineralogy and -XRF mapping at a resolution of approximately 1 m in 2-D.The rock samples consisted of fractured tuffaceous rock from Yucca Mountain, Nevada.A sample area of 10 6 m 2 was mapped in detail and the resulting element and porosity maps were digitized, thus creating a domain decomposed into 12,208 grid cells that were 8.77 m on a side (Fig. 10).A bulk porosity of about 6% was estimated based on averaging of the entire sample.Simulations were conducted in which a dilute fl uid enters the discretized porous medium at two different fl ow rates of 0.1 and 1.0 m 3 m -2 yr -1 , assumed to be uniform across the domain.The fl uid is reacted with the rock at 90 ºC.Simulations involving the slower fl ow rate, in which the fl uid residence times is approximately 3.65 days, provide fl uid composition results at the downstream end that are very similar to those obtained from homogeneous mineral distribution representations.At the higher fl ow rate of 1.0 m 3 m -2 yr -1 (residence time of approximately 8.76 hours), however, the fl uid composition differs between the heterogeneous and homogeneous cases along the entire length of the fl ow path.The authors concluded that the simulation results demonstrate that the fl uid composition characteristics in the homogeneous and discrete mineral representations will be similar only when the bulk average contact times for the individual mineral phases along the fl ow paths are approximately equivalent (within a few percent).

Micro-continuum reactive transport simulations in Lower Tuscaloosa Formation (Cranfi eld) sandstone.
The Cranfi eld Oil Field in Mississippi has been used as a subsurface CO 2 injection pilot site, with super-critical CO 2 injected into the lower Tuscaloosa Formation at about 300 m depth.The Tuscaloosa Formation is a 15-m-thick heterogeneous fl uvial sandstone that was the subject of an experimental study by Lu and co-workers (Lu et al. 2012), who reported low reactivity for the sandstone in contact with CO 2 -infused brine.The question arises as to whether the low reactivity is due primarily to the limited availability of reactive surface area?This can be evaluated more quantitatively with the micro-continuum approach.The bulk reactivity can be estimated by carrying out 2-D diffusion-reaction simulations using the volume-averaged porosity and mineral distributions presented in Figure 9, with the left hand boundary set as a fi xed or Dirichlet boundary condition.This effectively makes the CO 2 reservoir (5 bars, 25 ºC) infi nite.The simulations are carried out over the 4 mm by 4 mm domain for a period of 365 days, which The white regions map fine-grained groundmass, which was treated as inter-grown cristobalite+alkali feldspar+Ca-smectite+hematite.The gray sinuous band running from top to bottom of the figure is a silicafilled fracture.Area is 1000 × 1000 m.[Reproduced from Glassley WE, Simmons AM, Kercher JR (2002) Mineralogical heterogeneity in fractured, porous media and its representation in reactive transport models.Applied Geochemistry, Vol.17, p. 699-708, Fig. 2, with permission from Elsevier.] is a suffi cient amount of time for the system to come to a steady state.An effective diffusion coeffi cient of 10 -10 m 2 s -1 is assumed for the domain, with the exception of zones containing 50% or more chlorite where the effective diffusion coeffi cient is assumed to be 7 × 10 -12 in agreement with 1-D simulation shown in Figure 6B.No correction is made for the 3-D connectivity in the 2-D simulations, so access to some reactive phases is likely underestimated.How to incorporate 3-D tortuosity into 2-D modeling domains is a future research topic.In the case where the local porosity is zero, the diffusive fl ux will be zero because of the use of harmonic means to calculate properties at grid cell interfaces.Thus, unconnected porosity is automatically accounted for within the resolution of the discretization.Sub-grid connectivity, however, is not accounted for, although this should be possible with a more rigorous treatment of the data.
The spatial distribution of chlorite dissolution rates in the 2-D Tuscaloosa Sandstone section is shown in Figure 11.The sparse distribution of accessible chlorite certainly contributes to the low bulk reactivity of the material, so perhaps the observations Lu et al. (2012) of low reactivity are understandable.Bulk rates calculated from the micro-continuum simulations are given in Table 3.The rates are normalized to a cubic centimeter of Tuscaloosa Sandstone, which is close to the size of the rock sample used in the Lu et al. (2012) experiments.The low reaction rates of the lower Tuscaloosa Formation sample in contact with the CO 2 -infused brine is thus a consequence of both the low volume fractions of the reactive phases (quartz dominates) and the poor accessibility of some of the phases.

Multi-continuum approaches
There are many examples of soils, sediments, and rocks that are characterized by multiple length scales, each with its own set of physical and/or chemical properties.Probably the best known example is that of fractured rock.Here fl ow in the fractures is described by meter or larger length scales, while transport in fi ner-grained, unfractured material within the same volume is dominated by diffusion length scales (mm-cm).Other examples of hierarchical porous media can be mentioned, as for example where pore or smaller scale parameters and processes affect larger macroscale behavior.The models used to describe these systems are typically referred to as multi-continuum models.
In some cases, these scales are separated in space and can be represented as discrete spatial zones within a fl ow or reactive transport model-an example might be a discrete fracture located in what is otherwise an unfractured, low permeability material.In other cases, however, the domains may overlap within the same representative elementary volume (REV).In these cases, the domains are typically represented as two or more distinct continua with their own set of mass balance equations and physical-chemical properties (Barenblatt et al. 1960;Pruess and Narisimhan 1985).Functions describing exchange between the various continua may be included as well.Multi-continuum models had their origin primarily in fractured rock systems where it has been diffi cult or impossible to represent all of the fractures discretely with a single representative elementary volume.The general approach was fi rst introduced by Barenblatt andco-workers in 1960 (Barenblatt et al. 1960) and it has since been implemented in various forms, including: (1) the equivalent continuum model, or ECM (Wu 2000), ( 2) the dual permeability model (DPM), dual or multi-porosity model (Warren and Root 1963), and (3) multiple interacting continua or MINC approach (Pruess and Narisimhan 1985;Aradóttir et al. 2013).Among these three commonly used approaches, the dual-continuum approach has been most extensively applied in different subsurface environments (Arora et al. 2011), perhaps because it is relatively simple compared to the other approaches and because it is capable of describing many natural subsurface materials.The dual-continuum approach considers two interacting regions, one associated with the less permeable soil or rock matrix, and the other characterized by fl ow and/or diffusion in macropores or discrete fractures.In this approach, the representative elementary volume (REV) is partitioned into sub-volumes of each domain such that: where V REV refers to the total volume of the porous medium, V macro and V micro refer to the volume of the macro and micro (or equivalently, fracture and matrix) domains, respectively.Fe-hydroxide 2.8 × 10 -4 1.0 -2.9 × 10 -15 Calcite 0 0.04 0 Magnesite 0 0.04 0 Amorphous Silica 0 0.024 0 Table 3. Bulk rates calculated in 2-D section of Lower Tuscaloosa Formation sandstone using porosity and mineral abundance distributions shown in Figure 9. Intrinsic rates per unit surface area mineral used are standard literature values.Bulk porosity = 0.13.
The fraction, e, of volume occupied by the macropores and micropores, then, can be described respectively as (Lichtner 2000) These relations can be easily extended to include multiple interacting domains, as in the MINC approach (Pruess and Narisimhan 1985).The dual or multi-domain conceptualizations can be distinguished by their different formulations of the governing equations of fl ow in the fracture domain and/or by their different approaches to establish exchange between the two overlapping continua (or multiple domains).Various reviews of the different multi-continuum approaches, including the governing equations and exchange functions, are available elsewhere (Berkowitz and Balberg 1992;Lichtner 2000;Šimůnek et al. 2003;MacQuarrie and Mayer 2005;Aradóttir et al. 2013).
Multi-continuum models applied to micro-continuum systems.While the multicontinuum approach has been applied to various problems in which two or more domains with contrasting permeability and/or diffusivity are identifi able, the approach has less commonly been used to capture micro-continuum scale effects, for example, interactions between grain and pore scale processes (e.g., nm to mm scale) and the larger domain within which fl ow and reactive transport occurs (e.g., m scale).Perhaps the earliest contributions that considered interactions between microscopic and macroscopic domains via diffusion were those of Ortoleva and co-workers (Dewers and Ortoleva 1990;Sonnenthal and Ortoleva 1994).The case of diffusive exchange between a macroscopic melt and microscopic discrete crystals was considered by Wang (1993).
Wanner and Sonnenthal presented a three region model for kinetic Cr isotopic exchange (shown schematically in Fig. 12) that considered a mobile region within which advective fl ow occurs, an immobile region within which transport is only via diffusion, and a "mineral region" in which all of the reactions take place (Wanner and Sonnenthal 2013).The advantage of the MINC approach is that it is able to handle diffusive fl uxes, J D , between the minerals and  (Wanner and Sonnenthal 2013).By treating the mineral where the reduction occurs as a discrete domain, mixed diffusion-surface reaction rates can be considered.J D1 and J D2 refer to the fl uxes of Cr(VI) from the mobile and immobile continua to the mineral surfaces, respectively.both the other domains within which transport occurs, thus allowing for an explicit treatment of diffusive limitations to the rate (Xu 2008).Providing a surface reaction-controlled rate at the mineral surface in combination with the MINC approach allows one to consider a mixed diffusion-surface reaction control on the rate, as in the discrete model presented by Noiriel et al. (2012).
In a study by Aradóttir et al. (2013), the method of 'multiple interacting continua' (MINC) was applied to include microscopic rate-limiting processes operating at the grain scale within continuum (cm to m) scale reactive transport models of basaltic glass dissolution.In contrast with the nanometer-scale resolution model for glass dissolution discussed below, the approach taken by Aradóttir et al. (2013) allows for the use of a coarse numerical grid while capturing the interaction with the microscopic grains via the multi-continuum approach.The MINC method involves dividing the system up to ambient fluid and grains, using a specific surface area to describe the interface between the two (Fig. 13A).The various grains and regions within grains are then described by dividing them into continua separated by dividing surfaces.Millions of grains can thus be considered within the method without the need to explicitly discretize them individually.Four continua were used for describing a dissolving basaltic glass grain; the first one describes the ambient fluid around the grain, while the second, third and fourth continuum refer to a diffusive leached layer, the dissolving part of the grain and the inert part of the grain, respectively (Fig. 13B).

Resolution of nanoscale reaction fronts
In addition to their incorporation in MINC, micro-continuum models also have important application to the simulation of nanoscale reaction fronts, particularly where there has been interest in the long term performance of the engineered or natural materials (Grambow 2006;Gin et al. 2013b).Unlike the multi-continuum approach in which a relatively coarse numerical discretization is used (that approach relies on the multiple interacting continua to capture microscopic behavior), very high resolution gridding is used in this section to capture microscopic effects.There has been increasing interest in this topic in recent years as various characterization methods have dramatically improved the spatial resolution of reaction fronts that can be achieved.In particular, the higher resolution of the newer chemical profi ling techniques, which include atom probe tomography (APT), scanning transmission electron microscopy (STEM), energy fi ltered transmission microscopy (EFTEM) , and time-of-fl ight secondary ion mass spectrometry (ToF-SIMS), has called into question the long-standing model of glass and mineral dissolution in which diffusion and hydration lead to the selective release of cations from the surface-altered zone (Geisler et al. 2010;Hellmann et al. 2012Hellmann et al. , 2015;;Gin et al. 2013aGin et al. , 2015)).The higher resolution techniques demonstrate convincingly that the broad sigmoidal profi les interpreted as inter-diffusion cation profi les are largely an effect of the low resolution techniques that average elemental concentrations across a broad region.The broad inter-diffusion profi les in glass and minerals that investigators thought they saw in the past had led to models for dissolution that involved selective leaching of elements from the glass or mineral structure.In contrast, Hellmann et al. (2015) report nm to sub-nm scale reaction front widths in altered borosilicate glass for all ions except for H + , the measurement of which they suspect to be subject to too much error for high resolution mapping.The sharp nm-scale fronts indicated by higher resolution profi ling, along with isotopic studies targeting the gel layers formed from glass corrosion (Geisler et al. 2010), have led to reinterpretation of these as dissolution-precipitation rather than inter-diffusion fronts (Hellmann et al. 2012(Hellmann et al. , 2015)).It should be pointed out that some high resolution studies like that on the nuclear glass altered for 25 years (reported by Gin et al. 2013a, see Fig. 14A) still indicate ~20-nm fronts for H + and Li + , suggesting that at longer time scales, it may still be possible for the inter-diffusion fronts to develop, even if they are much narrower than previously thought.It is noteworthy that even in the case of the 25-year glass investigated by Gin et al (2013a), however, the reaction front for boron and sodium are narrower than the fronts for Li + and H + , arguing that a dissolutionprecipitation mechanism c ontrols the release of B and Si (Fig. 14B).

Analytical and numerical models for reaction fronts.
A number of models for diffusion and reaction have been presented over the years, with noteworthy contributions by Thompson, Korzhinskii, and Weare (Korzhinskii 1959;Thompson 1959;Weare et al. 1976).A special class of analytical models have been developed for the case of inter-diffusion of cations applied primarily to the problem of nuclear glass corrosion (Doremus 1975;Hellmann 1997).Lichtner et al. (1986) presented perhaps the fi rst numerical reaction-diffusion model that could accommodate kinetic models for mineral (or glass) reaction (Lichtner et al. 1986).The fi rst comprehensive numerical study of a geological diffusion-kinetic reaction system may have been that presented by Steefel and Lichtner in 1994, a study that highlights some of the advantages of this approach over those relying on analytical solutions to the diffusion or diffusion-reaction system (Steefel and Lichtner 1994).The advantages of the numerical versus analytical models is their ability to couple multicomponent diffusion and kinetically controlled mineral reaction while considering aqueous and surface complexation.In addition, the grid-based numerical formulation allows one to consider changes in porosity, diffusivity, and permeability resulting from the chemical reactions.
Numerical models for glass alteration-the GRAAL model.While numerical diffusionreaction models are now becoming more common, in part because their advantages in coupling multicomponent transport and reaction as discussed above (see, for example, Marty et al. 2015), applications to the micro-or nanoscale are more rare.The applications to date have primarily been to glass corrosion, although the approaches should be applicable to mineral dissolution as well.A noteworthy fi rst effort in this regard are the series of papers by Frugier and co-workers that describe the basis for their glass corrosion model (referred to as the GRAAL model) and its application to the problem of glass corrosion (Frugier et al. 2008).The key features of glass corrosion implemented in the GRAAL (glass reactivity with allowance for the alteration layer) model include (Fig. 15): • Relatively rapid exchange and hydrolysis reactions involving the mobile glass constituents (alkalis, boron, etc.); • Slower hydrolysis involving silicon, which results in the formation of an amorphous silica-rich layer at the glass/solution interface; • The amorphous layer itself dissolves as long as the external solution is undersaturated with respect to silica; • The amorphous layer becomes a barrier to diffusion (referred to as a Passive Reactive Interface, or PRI), which at steady-state becomes the rate-limiting mechanism.The PRI is assumed to form at the interface between the gel layer and the inter-diffusion zone and it represents a densifi cation of a portion of the gel layer; • Crystallization of secondary phases may occur in the broader gel layer.
The principal simplying assumptions in the GRAAL model are: • The rate-limiting step for glass corrosion is water diffusion within the PRI; • Water diffusion in the glass and proton/alaklin ion exchange are ignored, and; • The reactivity of the PRI is described by its thermodynamic state relative to the leaching solution.
The model equations that are solved, then, include: 1.An equation describing the kinetics of dissolution of the PRI as a result of undersaturation with respect to the solution; 2. An equation describing the kinetics of formation of the PRI as a result of water diffusion; 3. An equation describing the kinetics of formation of secondary phases in the gel; 4. Mass balance for silicon, and: 5. Mass balance for boron.
Numerical models for glass alteration-the KC model.As an alternative to the GRAAL model in which diffusion through the Passive Reactive Interface (PRI) is assumed to be rate-limiting, it is possible to develop a more general model that makes no a priori assumptions about the rate-limiting step in the overall glass alteration process.The Kinetic Micro-Continuum (KC) model that has been developed is based on the reactive transport software CrunchFlow (Steefel et al. 2015) and includes the following processes: • Diffusion of water through the pristine glass and its alteration products; • Ion exchange between water and the cations in the glass; • Kinetically controlled hydrolysis reactions resulting in breaking of glass network bonds (Si, B, Al, etc.).The rate may be described by either a linear or a nonlinear transition state theory (TST) law with an affi nity control supplied by a specifi c phase (e.g., amorphous silica), or with an irreversible rate law with no affi nity control.In either case, far from equilibrium dependencies of the rate on other dissolved (e.g., pH, Al, silica) or sorbed species can be included; • Multicomponent diffusion of ions through the glass corrosion products; • Precipitation reactions for amorphous and/or crystalline phases of variable composition that are kinetically and thermodynamically controlled; • Kinetically controlled ripening and/or densifi cation reactions that can reduce the porosity and/or pore connectivity (and thus the diffusivity) of the corrosion products; • Kinetically and thermodynamically controlled formation of new crystalline phases (e.g., smectite, zeolite), with possible consequences for the transport properties of the corrosion layer; • Flow and diffusion in the aqueous phase adjacent to the glass surface.
The KC model incorporates the possibility (unlike the requirement in the GRAAL model) of diffusion-limited glass corrosion by considering explicitly the kinetically-controlled densifi cation of either (1) a residual silica-rich glass network in which other important components (e.g., the cations and network former boron) have been leached, or (2) of a newly precipitated silica-rich gel layer.Whether a passivating layer (i.e., defi ned as the Passivating Reactive Interface (PRI) by Frugier et al. 2008) forms in the model depends on the relative rates of (1) silica recrystallization and densifi cation, (2) leaching of the glass constituents, and (3) dissolution and/or recrystallization of the corrosion products.
Application to the 25-year glass alteration test.As an application of the KC model described above, the 25-year glass alteration experiment as described by Guittoneau and coworkers (Guittonneau et al. 2011;Gin et al. 2013a) is simulated.The principal objective of the modeling is to capture the width of the various reaction zones close to the pristine glass surface and their relative positions as recorded by Atom Probe Tomography (APT) rather than to match the long term corrosion rate.We assume a pure diffusion-controlled regime and a constant grid spacing of 1 nm.The assumption that a continuum model applies at this spatial scale is a severe approximation given that pore sizes are close to this value, but the approach allows us to compare results with elemental profi les and avoids the additional requirement of a full computationally expensive and chemically simplifi ed atomistic treatment (Bourg and Steefel 2012).
As a boundary condition at one end of the reactor-glass specimen system, we consider a Dirichlet or fi xed concentration condition corresponding to the mineral water used to replenish periodically the experimental reactor (Guittonneau et al. 2011).The fi xed concentration boundary condition in 1-D is probably a good approximation to a fl ow-through or continuously fl ushed system.At the other end of the 1-D system, we assume a no-fl ux condition, which is reasonable as long as the corrosion front does not fully penetrate the glass specimen.Within the fi rst 50 nm of the reaction, the system is characterized by a porosity of 0.41 (as in the experimental system reported by Guittonneau et al. 2011) and a mixture of quartz sand and granite upon which the borosilicate glass coupon rests.A diffusivity of 10-12 m 2 s -1 for all ions was assumed for the sand-granite mixture.The alloy specimens included in the experiments were not considered in the modeling.From 50 nm out to 200 nm, the system was assumed to consist of a borosilicate glass with a porosity of 1%.The diffusivity of the borosilicate glass was assumed to follow a threshold type of model (Navarre-Sitchler et al. 2009), with a value of 5 × 10 -24 m 2 s -1 (in approximate agreement with the value proposed by Gin et al. 2013a) for values of the porosity below 50% and a value of 10-12 m 2 s -1 for porosity values above 50%.Modeling carried out on weathered basalts (Navarre-Sitchler et al. 2011) indicate that a simple porosity dependence (as in an unmodifi ed Archie's Law formulation) cannot replicate the observed concentration profi les, since the reaction front continuously widens due to the simulated porosity and diffusivity enhancement.Some form of a threshold model, based on the idea that dissolution and porosity enhancement increase the rate of diffusivity by increasing connectivity (Navarre-Sitchler et al. 2009, 2011), appears to be required.
In the modeling, the dissolution of the glass is assumed to follow a Transition State Theory (TST) rate law with a dependence on the saturation state (departure from equilibrium, or affi nity) with respect to amorphous silica (Grambow 2006).Since a linear TST dependence on the departure from equilibrium does not capture the sharp B front and places the B (and Na) dissolution fronts too close to the Li-H inter-diffusion front, the dissolution of the glass is assumed to have a cubic dependence on the departure from equilibrium with respect to amorphous silica.The rate of glass corrosion is also assumed to depend on the hydration state of the glass: before the H 2 O diffusion front has penetrated the pristine glass and hydrates it, the rate is effectively zero.A higher-order dependence on the concentration of hydrated sites in the borosilicate glass is also required so as to locate the boron release front further from the Li-H inter-diffusion front.The rate law used is therefore where k is the rate constant and a H-hydrated is the concentration of hydrated sites in the glass, and Q am-silica and K am-silica are the ion activity products and equilibrium constants with respect to amorphous silica.This formulation could be reconciled with a model in which the number of hydrated sites needs to reach some (high) threshold value before the dissolution of the glass accelerates appreciably.
Here we present a semi-quantitative comparison of the simulation results from the KC model with nanometer discretization to the data from the Gin et al. (2013) study.The focus is on the relative position of the fronts, and in general, the width of the fronts as they evolve over time rather than on the total extent of alteration or even the rate of alteration.Schematically, the geometry that we wish to capture in the modeling is given in Figure 14 above (Gin et al. 2013a).
The model results for the 1-D run after three years are shown in Figure 16A.The simulations predict that the Li-H inter-diffusion front maintains a relatively constant width of about 20 nm over three years (time evolution not shown).The B release (dissolution) front is even sharper and is located further from the pristine glass interface than is the Li-H interdiffusion front, a result that can be justifi ed based on the assumption that the dissolution of the glass occurs rapidly when hydration of the glass is nearly complete.Thus, the key component of the model is the coupling of glass dissolution to the extent of hydration, which is driven by diffusion into the pristine glass (not the PRI discussed by Frugier et al. 2008).According to the simulations, the rate-limiting step for the overall glass alteration process is the diffusion into the pristine glass-once hydrated, the borosilicate glass dissolves quickly, as indicated by its sharp front.The simulations also predict an early time period when the steady-state 20-nm inter-diffusion zone is not fully developed (results not shown), which might help to reconcile the observations by Hellmann et al. (2015)   In addition, the simulations predict a linear rate of front propagation over time once the initial period (less than 1 year) is passed (Fig. 16B).This can be explained by a constant-width zone over which diffusion is limiting, in agreement with earlier results on weathering of basalt (Navarre-Sitchler et al. 2011).In the case of the nuclear glass altered for 25 years, the constant width is a result of the development of an approximately 20-nm-wide inter-diffusion zone at the edge of the pristine glass.Between this zone and the dissolving gel at the outer boundary, the increase in porosity is suffi cient that diffusion is not limiting (even if the zone were to grow with time).If diffusion through a continuously growing silica-rich gel layer was limiting the rate, the dependence on time should be parabolic.A constant width PRI, as discussed by Frugier et al (2008), could also result in a linear front advance rate.

SUMMARY AND PATH FORWARD
While true pore-scale models are arguably the most rigorous way to treat geochemical processes operating at the pore-scale, micro-continuum modeling approaches offer some advantages in terms of their relative ease of use, ability to apply well-tested software (e.g., Steefel et al. 2015), and computational effi ciency.This approach offers the additional advantage that micron to even nm-scale mineralogical, chemical, and physical heterogeneities can be incorporated into the simulations.The disadvantage of the approach is that one still faces many of the standard limitations of continuum representations of the pore scale, namely the need to average geochemical, mineralogical, and physical properties and the inability to explicitly resolve interfaces between solids, gases, and fl uids.Many important parameters and processes still operate at the sub-grid scale (e.g., nanopore connectivity) and these must be accounted for in order to achieve a realistic simulation of pore-scale geochemical processes.The challenge of dealing with reactive surface area, for example, persists in any continuum treatment of the pore-scale, in contrast to the more rigorous geometric methods in true pore-scale models where the fl uid-mineral interface is resolved explicitly (Molins 2015, this volume).
Certainly the advent of new microscopic characterization techniques, including increasingly higher resolution X-ray microtomography, BSE-SEM, and FIB-SEM, are motivating the search for novel and complementary modeling methods.The longer time and space scales that achievable with the micro-continuum models is another reason why they will not soon be replaced by either molecular dynamics (MD) modeling approaches or even true pore-scale models.The true pore-scale and MD approaches have an important role to play here, however, since they can be used to provide upscaled parameters for the micro-continuum models.Eventually, one expects the development of a new class of hybrid models that link the molecular, true pore-, and micro-continuum scales within a single dynamic, multi-scale framework.
Much remains to be done in the fi eld of micro-continuum modeling of pore-scale geochemical processes.In fact, the fi eld is still in its infancy.Since the mineralogical mapping is predominantly carried out in 1-D or 2-D, we require an improved treatment of how 3-D effects and parameters are incorporated.Tortuosity and permeability are two of the most obvious examples.We also need improved representations of the correlation between mineralogy and physical transport properties like diffusivity and permeability at the nanoand micro-scale.Consistent upscaling strategies are required so that it is possible to change model resolution without undue loss of process fi delity.Ultimately, it is clear that the microcontinuum approaches will need to incorporate a more formal multi-scale framework, particularly where there is interest in capturing nanoscale features like pore connectivity and reaction fronts within larger scale domains.

Figure 1 .
Figure 1.A) Example of a through diffusion cell setup: (a) inlet reservoir, (b) peristaltic pump, (c) throughdiffusion cell, and (d) outlet reservoir.Arrow heads indicate the circulation of water from the reservoir to the fi lter in order to homogenize the inlet and outlet solutions compositions.B) Plot of concentrations in inlet and outlet reservoirs versus time.[Reproduced from Tachi Y, Yotsuji K (2014) Diffusion and sorption of Cs + , Na + , I − and HTO in compacted sodium montmorillonite as a function of porewater salinity: Integrated sorption and diffusion model.Geochimica et Cosmochimica Acta, Vol.132, p. 75-93, Figs. 1 and 3, with permission from Elsevier.] (Navarre-Sitchler et al. 2009).Navarre-Sitchler et al. (2009) took the additional step of comparing these results to estimates of the diffusivity from pore-network modeling (see discussion below).
(basalt in blue, pores in red).Results of the numerical tracer experiment are shown in Figure 3B, with diffusion of the tracer from the bottom of the Figure 3B towards the top.Note that in these 3-D tracer diffusion simulations, only two distinct tortuosities (or diffusion coeffi cients) are

Figure 2 .
Figure 2. Bromide tracer (blue) diffusion profi le after 7 days based on XRF mapping at the Advanced Light Source, Lawrence Berkeley National Laboratory.Diffusion is from left to right.See Navarre-Sitchler et al. (2009) for discussion of simulations.
2 sequestration site investigated by Landrot et al. 2012.One such pore of approximately 2 m diameter imaged with scanning electron microscopy (SEM) is shown in Figure 5A.The Figure shows a quartz grain on the right side partly milled with focused ion beam techniques, while Figure 5B shows

Figure 3 .
Figure 3. A) Segmented X-ray synchrotron microtomographic data collected at the Advanced Light Source at Lawrence Berkeley National Laboratory with a voxel resolution of 4.4 m.Macropores developed as a result of chemical weathering in the basalt and are connected primarily in the third dimension (into the page), with red indicating pores, blue indicating basalt.B) Tracer diffusion simulation results using the pore structure shown in Figure 3A.Results are shown after 7 days of diffusion of the tracer from bottom to top.Simulation assumed a Dirichlet boundary condition at the bottom with a fi xed tracer concentration of 0.01.See Navarre-Sitchler et al. (2009) for description of experiments.

Figure 4 .
Figure 4. Comparison of results from laboratory and numerical (simulated) diffusion experiments.A) The XRF image of Br concentrations measured in basalt samples after 7 days (high concentration in white, low concentration in darker shades).B) Contour plot of simulated Br (tracer) distribution after 7 days in the same sample based on pore structure determined with X-ray synchrotron microtomography.C) 1-D effective diffusion coeffi cient fi t to the 3-D data.[Reproduced from Navarre-Sitchler A, Steefel CI, Yang L, Tomutsa L, Brantley SL (2009) The evolution of dissolution patterns: Permeability change due to coupled fl ow and reaction.In: Chemical Modeling of Aqueous Systems II.Vol 416.Melchior D, Bassett RL (eds)American Chemical Society, Washington, p 212-225, Fig. 9 with permission from the American Geophysical Union.]

Figure 6 .
Figure 6.A) Simulated 2-D slice of tracer concentration from a 3-D cube of chlorite and pores from the Lower Tuscaloosa Formation at the Cranfi eld CO 2 sequestration site, which was characterized with FIB-SEM (Landrot et al. 2012).The code CrunchFlow (Steefel et al. 2015) is used to carry out a 3-D tracer diffusion simulation, with release of the tracer at the left boundary at a concentration of 0.01 M. B) Fit of 1-D diffusion profi le assuming a single homogeneous diffusivity of 7 × 10 -12 m 2 /s (line) versus concentration from 3-D simulation.

Figure 7 .
Figure 7. Steady-state fl ow through fractured shale based on Navier-Stokes equation.The surface plot shows the interface between pores and solid (shale).The surface is colored according to velocity magnitude in the pore, with blue representing lower velocities, shades of green, yellow, and red representing higher velocities.Resolution 48.4 nm with 1920 × 1600 × 640 total grid cells.Original pore structure obtained from FIB-SEM mapping, which indicates a bulk porosity of 18%.[Reproduced from Trebotich D, Graves DT (2015) An adaptive fi nite volume method for the incompressible Navier-Stokes equations in complex geometries Communications in Applied Mathematics and Computational Science, Vol.10(1), p. 43-82, Fig. 13, with permission from Mathematical Sciences Publishers.]

Figure 8 .Figure
Figure 8. Connected porosity mapping in a 29 mm 2 region of the Lower Tuscaloosa Formation sandstone (Cranfi eld), considering the nanoporosity within chlorite as 100% connected.The porosity is mapped with a 331 nm/pixel resolution.The white color represents the connected pore network that starts from the edges of the image and propagates inward through the map, and the turquoiseblue color represents the chlorite fraction that is linked to the connected pore network measured in the 29 mm 2 region.[Reproduced from Landrot G, Ajo-Franklin J, Yang L, Cabrini S, Steefel CI (2012) Measurement of accessible reactive surface area in a sandstone, with application to CO 2 mineralization.Chemical Geology, Fig. 8, with  permission from Elsevier.]

Figure 10 .
Figure10.Discretized mineralogy of tuff sample generated from digitized element maps.Gray regions represent pure silica (treated as cristobalite in the simulations), and the black regions are K-feldspar.The white regions map fine-grained groundmass, which was treated as inter-grown cristobalite+alkali feldspar+Ca-smectite+hematite.The gray sinuous band running from top to bottom of the figure is a silicafilled fracture.Area is 1000 × 1000 m.[Reproduced fromGlassley WE, Simmons AM, Kercher JR (2002) Mineralogical heterogeneity in fractured, porous media and its representation in reactive transport models.Applied Geochemistry, Vol.17, p. 699-708, Fig.2, with permission fromElsevier.]

Figure 11 .
Figure 11.Spatial distribution of chlorite dissolution rates in units of mol L -1 fl uid s -1 after 365 days of diffusion-reaction (no fl ow) simulation.Given these results, the low reactivity of the Tuscaloosa Sandstone samples investigated by Lu et al. (2012) is perhaps understandable.

Figure 12 .
Figure 12.Schematic representation of multiple-interacting continua (MINC) representation of kinetic Cr isotopic fractionation(Wanner and Sonnenthal 2013).By treating the mineral where the reduction occurs as a discrete domain, mixed diffusion-surface reaction rates can be considered.J D1 and J D2 refer to the fl uxes of Cr(VI) from the mobile and immobile continua to the mineral surfaces, respectively.

Figure 13 .
Figure 13.A) A four-dimensional MINC interpretation of basaltic glass dissolution.The left figure shows a zoom-in of real grains in the simulated column, which is packed with basaltic glass grains of size 125-250 m, yielding a porosity of 0.45.The middle figure shows a blow up of several grain clusters within the column and their interpretation as four interacting continua within the MINC approach.Each grain cluster consists of approximately 25,000 individual grains.B) A schematic illustration of elements and connections in the four-dimensional MINC setup, with each column representing a different continuum.[Reproduced from Aradóttir ESP, Sigfússon B, Sonnenthal EL, Björnsson G, Jónsson H (2013) Dynamics of basaltic glass dissolution-Capturing microscopic effects in continuum scale models.Geochimica et Cosmochimica Acta, Vol.121, p. 311-327, Figs. 2 and 5 with permission from Elsevier.]

Figure 14 .
Figure 14.A) High resolution atom probe tomography (APT) profi le across glass alteration front.Note sharper front for B and Na as compared to Li and H. B) Schematic representation of distribution of fronts for 25-year altered glass shown on left.[Modifi ed slightly after Gin S, Ryan JV, Schreiber DK, Neeway J, Cabié M (2013a) Contribution of atom-probe tomography to a better understanding of glass alteration mechanisms: Application to a nuclear glass specimen altered 25 years in a granitic environment.Chemical Geology, Vol.349, p. 99-109, Figs. 3 and 6, reproduced with permission of Elsevier.]

Figure 15 .
Figure 15.Simplifi ed diagram of predominant mechanisms accounted for in the GRAAL model for glass corrosion.The Passive Reactive Interface (or PRI) is interpreted to form at the interface between gel and the ion inter-diffusion zone (not shown) and is assumed to be the rate-limiting step in the overall glass corrosion process.[Reproduced from Frugier P, Gin S, Minet Y, Chave T, Bonin B, Godon N, Lartigue J-E, Jollivet P, Ayral A, De Windt L (2008) SON68 nuclear glass dissolution kinetics: Current state of knowledge and basis of the new GRAAL model.Journal of Nuclear Materials, Vol.380, p. 8-21, Fig. 13, with permission from Elsevier.] of a nanometer to even sub-nanometer Li + front in their shorter term experiments (recall that the Gin et al. 2013 study was based on experiments conducted over nearly 25 years, as described in Guittoneau et al. 2011).

Figure 16 .
Figure16.A) Simulation results using the KC model, with alteration proceeding from left to right.The original glass wafer edge is located at 50 nm.Note the position of the boron release front further from the pristine glass than the Li-H inter-diffusion front, in qualitative agreement with the observations shown in Figure14.The simulations predict that the Li-H inter-diffusion front maintains a relatively constant width of about 20 nm over three years (time evolution not shown).B) Position of borosilicate glass corrosion front as a function of time.The linear advance rate is only compatible with a model in which diffusivity is enhanced within the altered zone, resulting in a nearly constant width reaction zone over time.

Table 1 .
Results of fi tting of tracer distribution from 3-D simulations with a 1-D diffusion model.

Table 2 .
(Landrot et al. 2012)raging of porosity and mineral percentages from Lower Tuscaloosa Formation (Cranfi eld) sandstone.Volume averaging to 16 m is based on an original image resolution of 4 m(Landrot et al. 2012).Each pixel in the table corresponds to a 2-D section measuring 16 m 2 .