This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.

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 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, 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–field 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 mineral dissolution (e.g., feldspar) can result in an approach to thermodynamic equilibrium with respect to the dissolving primary phase, and thus a slowing of the rate of reaction (Zhu et al. 2004; Maher et al. 2006, 2009). A second potentially important effect is related to the heterogeneous nature of natural porous media, which can result in bypassing of reactive zones by groundwater and surface water flow (Malmström et al. 2000). At the pore scale, this effect can occur when pores and/or micro-environments are not connected to the macropores within which most of the flow occurs. The result can be that minerals lining these pores contribute little or nothing to the overall reactivity of the formation (Peters 2009; Landrot et al. 2012). In addition, diffusion boundary layers can form around reactive grains, further reducing the rate of reaction relative to the experimental surface reaction rates determined in the absence of these transport effects (Li et al. 2008; Noiriel et al. 2012; Molins et al. 2014).

One can argue that the ability to demonstrate a predictive capability for geochemical processes, including those operating at the pore scale, is the ultimate test of understanding (Steefel et al. 2005). The development of new characterization techniques for the pore scale has important implications for the models that can be applied to these systems. One possibility is to import the microscopic characterization or mapping as initial or final conditions into true pore-scale models (Molins et al. 2012, 2014, 2015, this volume; Steefel et al. 2013; Yoon et al. 2015, this volume). Another option that typically allows for larger spatial domains to be considered is the use of the characterization data in pore network models (Mehmani and Balhoff 2015, this volume). A third possible approach that is summarized in this chapter is to make use of “micro-continuum” models informed by high-resolution geochemical, mineralogical, and physical data to describe geochemical pore-scale processes. Micro-continuum geochemical models are typically coarser than either the true pore-scale or pore-network models and thus cannot resolve pore-scale interfaces between mineral, liquid, and gas. The approach suffers from most of the same limitations that apply to larger scale continuum descriptions of porous media, namely the inability to resolve pore-scale solid–liquid–gas interfaces and the requirement that many parameters and properties (e.g., permeability or reactive surface area) need to be averaged or upscaled in some fashion. However, the approach is capable of improving on coarsely resolved (meter-scale) models by assigning differing mineralogical/geochemical and physical properties (porosity, permeability, and diffusivity) values to the domain, thus making it possible to calculate larger scale (bulk) reaction rates and transport properties.

MAPPING OF MODEL PARAMETERS FROM IMAGE ANALYSIS

An important first step in developing micro-continuum pore-scale geochemical models is the collection and interpretation of data on the mineralogical, geochemical, and transport properties at a fine (< mm) scale. Detailed reviews are provided elsewhere on the range of characterization techniques available to describe pore-scale geochemical processes (Anovitz and Cole 2015, this volume; Navarre-Sitchler et al. 2015, this volume; Noiriel 2015, this volume). Here we focus on the approaches that are specifically suited for the development of pore-scale parameter distributions for micro-continuum modeling, although many of the techniques could also be used for direct pore-scale or pore-network modeling where the resolution is sufficiently high.

The pore-scale parameters of interest for micro-continuum modeling include the porosity, mineral volume fractions, and mineral reactive surface area, along with the more challenging transport-related parameters of permeability and diffusivity. While porosity is typically considered as a scalar quantity and therefore relatively easy to quantify with a variety of mapping/characterization techniques, the more important quantity for the purposes of reactive transport modeling is the connected porosity (Navarre-Sitchler et al. 2009; Peters 2009; Landrot et al. 2012). As demonstrated with modeling of basalt weathering over hundreds of thousands of years, the connected porosity is the parameter that controls reactivity under open 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 fluids 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 significantly 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 specific 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 briefly 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 significant contrast with the minerals.

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 sufficient 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 find 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 first 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 filtering 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 sufficient 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 significant 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 identification 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. Identification 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 misidentification 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 fluids if the porosity is imaged as well, and if the imaging is at sufficiently high resolution that the connectivity of the pore network can be quantified. As observed in recent studies (Peters 2009; Landrot et al. 2012), mineral abundance alone may not always accurately reflect 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 first 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 briefly 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), specific 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, specific, 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-specific 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).

Specific surface areas that reflect surface roughness can be estimated from GSA determinations by applying a roughness factor, although the choice of a factor is difficult 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 identified in 2-D images. In 3-D images, specific 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 fixed 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 flow). 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 profile. For example, Navarre-Sitchler et al. (2009) used micro-X-ray synchrotron mapping of bromide fluorescence (μXRF) to determine the diffusivity of samples of unweathered and weathered basaltic andesite. An effective, upscaled diffusion coefficient was determined for the porous andesite by fitting the profile with a numerical solution of the diffusion equation

 

Cit=-·[φDi*Ci]
(1)

where Di* is the diffusion coefficient for species i in porous media that incorporates the effects of tortuosity, Ci 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 (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).

Diffusivity from numerical experiments

Experiments and/or diffusion profile 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 fitted 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 fit 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 m2 s−1 for the largely unconnected pore structure of unaltered basalt and a free ion diffusivity of 10−9 m2 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 (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 coefficients) are used and they are based on the segmented 3-D porosity map of the weathered sample. The upscaled effective diffusion coefficient, in contrast, is determined by fitting a 1-D profile to the 3-D simulation results. The diffusion of the bromide tracer was assumed to follow a modified Archie’s Law model that incorporates a critical porosity threshold value of 9%, below which the porosity is considered to be largely unconnected:

 

De=Dp+D0(φe)2,
(2)

w here Dp is the diffusivity in the unweathered parent basalt, D0 is the diffusion coefficient in pure water, and De is the effective diffusion coefficient in porous media and where the effective porosity, φe, is defined as:

 

φe=a(φT-φc)β;φT>φcφe=0;φT<φc.
(3)

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 chlorite-filled pores within the Lower Tuscaloosa Formation that hosts the Cranfield CO2 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 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 difficult, but determining an approximate upscaled diffusion coefficient is relatively simple using a true pore-scale 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-flux. Based on the average diffusion profile, one can estimate an upscaled, effective diffusion coefficient of 7 × 10−12 m2/s (Fig. 6). The upscaled diffusion coefficient 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 pore-throat 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 defined 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 fluid 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 fluid, 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 flow 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 flow 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 filled 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 multi-dimensional 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 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 reflects some combinations of the neighboring voxel compositions (Ketcham and Carlson 2001). Additional difficulties 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 flat, 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 flipping of isolated misidentified phases Peters (2009). Others artifacts are not easily correctable and require additional more complex filtering 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 define 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 significant 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 sufficient 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 high-resolution 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 small-scale 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 flow (e.g., Blunt et al. 2013).

MICRO-CONTINUUM MODELING APPROACHES

If the imaging methods are sufficiently 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 field) 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 influenced 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 Cranfield CO2 sequestration site with a resolution of 331 nm (Fig. 8), which is adequate to resolve all but the nano-crystalline chlorite-filled pores (Landrot et al. 2012). With sufficient computational resources (software and hardware), it should be possible to simulate the pore-scale 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 μm2 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 first micro-continuum 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 first 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 106 μm2 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 fluid enters the discretized porous medium at two different flow rates of 0.1 and 1.0 m3 m−2 yr−1, assumed to be uniform across the domain. The fluid is reacted with the rock at 90 ºC. Simulations involving the slower flow rate, in which the fluid residence times is approximately 3.65 days, provide fluid composition results at the downstream end that are very similar to those obtained from homogeneous mineral distribution representations. At the higher flow rate of 1.0 m3 m−2 yr−1 (residence time of approximately 8.76 hours), however, the fluid composition differs between the heterogeneous and homogeneous cases along the entire length of the flow path. The authors concluded that the simulation results demonstrate that the fluid 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 flow paths are approximately equivalent (within a few percent).

Micro-continuum reactive transport simulations in Lower Tuscaloosa Formation (Cranfield) sandstone

The Cranfield Oil Field in Mississippi has been used as a subsurface CO2 injection pilot site, with super-critical CO2 injected into the lower Tuscaloosa Formation at about 300 m depth. The Tuscaloosa Formation is a 15-m-thick heterogeneous fluvial 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 CO2-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 fixed or Dirichlet boundary condition. This effectively makes the CO2 reservoir (5 bars, 25 ºC) infinite. The simulations are carried out over the 4 mm by 4 mm domain for a period of 365 days, which is a sufficient amount of time for the system to come to a steady state. An effective diffusion coefficient of 10−10 m2 s−1 is assumed for the domain, with the exception of zones containing 50% or more chlorite where the effective diffusion coefficient 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 flux 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 CO2-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 flow in the fractures is described by meter or larger length scales, while transport in finer-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 flow 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 difficult or impossible to represent all of the fractures discretely with a single representative elementary volume. The general approach was first introduced by Barenblatt and co-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 flow 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:

 

VREV=Vmacro+Vmicro,
(4)

where VREV refers to the total volume of the porous medium, Vmacro and Vmicro refer to the volume of the macro and micro (or equivalently, fracture and matrix) domains, respectively. The fraction, e, of volume occupied by the macropores and micropores, then, can be described respectively as (Lichtner 2000)

 

ɛmacro=VmacroVREV,ɛmicro=VmicroVREV.
(5)

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 flow 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 multi-continuum approach has been applied to various problems in which two or more domains with contrasting permeability and/or diffusivity are identifiable, 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 flow 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 flow 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 fluxes, JD, between the minerals and 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 profiling techniques, which include atom probe tomography (APT), scanning transmission electron microscopy (STEM), energy filtered transmission microscopy (EFTEM), and time-of-flight 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. 2012, 2015; Gin et al. 2013a, 2015). The higher resolution techniques demonstrate convincingly that the broad sigmoidal profiles interpreted as inter-diffusion cation profiles are largely an effect of the low resolution techniques that average elemental concentrations across a broad region. The broad inter-diffusion profiles 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 profiling, 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, 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 dissolution– precipitation mechanism controls 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 first numerical reaction–diffusion model that could accommodate kinetic models for mineral (or glass) reaction (Lichtner et al. 1986). The first 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 diffusion-reaction 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 first 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 densification 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 affinity control supplied by a specific phase (e.g., amorphous silica), or with an irreversible rate law with no affinity 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 densification 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 densification 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., defined 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 densification, (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 profiles and avoids the additional requirement of a full computationally expensive and chemically simplified atomistic treatment (Bourg and Steefel 2012).

As a boundary condition at one end of the reactor-glass specimen system, we consider a Dirichlet or fixed concentration condition corresponding to the mineral water used to replenish periodically the experimental reactor (Guittonneau et al. 2011). The fixed concentration boundary condition in 1-D is probably a good approximation to a flow-through or continuously flushed system. At the other end of the 1-D system, we assume a no-flux condition, which is reasonable as long as the corrosion front does not fully penetrate the glass specimen. Within the first 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 m2 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 m2 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 m2 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 unmodified Archie’s Law formulation) cannot replicate the observed concentration profiles, 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 affinity) 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 H2O 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

 

Rcorr=kaH-hydrated10(1-Qam-silicaKam-silica)3
(6)

where k is the rate constant and aH-hydrated is the concentration of hydrated sites in the glass, and Qam-silica and Kam-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 inter-diffusion front, a result that can be justified 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) 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).

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 sufficient 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 efficiency. 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 fluids. 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 fluid–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 field of micro-continuum modeling of pore-scale geochemical processes. In fact, the field 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 fidelity. Ultimately, it is clear that the micro-continuum 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.

ACKNOWLEDGMENTS

This work was supported as part of the Center for Nanoscale Control of Geologic CO2, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Contract No. DE-AC02-05CH11231 to Lawrence Berkeley National Laboratory. This work was also supported in part by the Director, Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, of the U.S. Department of Energy under the same contract to Lawrence Berkeley National Laboratory. We are grateful to Qingyun Li, Jennifer Druhan, and Bhavna Arora for their careful reviews of the chapter.