Dissolution and precipitation reactions are the primary mechanisms that occur when a rock (i.e., a mineral assemblage) is in contact with a fluid out of equilibrium. They play a critical role in natural processes (e.g., weathering, compaction, meteoric and marine diagenesis) and anthropogenic processes (e.g., reservoir acidizing, CO2 sequestration, acid mine drainage, contaminant mobility, bioremediation). Such fluid–rock interactions result in complex changes in pore structure and mineral composition, leading in turn to changes in flow, mechanical, and transport properties, such as permeability, dispersivity, strength, and pore accessibility. Consequently, geochemical disequilibrium can lead to complex modifications of hydrodynamic and transport properties in porous and fractured rocks.
Porous rocks are often characterized by complex textures and mineral compositions that are derived from their depositional and diagenetic environments. They typically have heterogeneous structures, the macroscopic physical properties of which depend on microscopic characteristics. Permeability, for example, is closely related to the microstructure, in particular the size and the spatial distribution of pore throats, pore roughness, and presence of fine clogging particles. The coupled hydrological, mechanical, and chemical (HMC) processes are highly non-linear and minor changes at the pore scale in one property can result in large modifications of the others properties. Prediction of system response to chemical conditions requires understanding how individual processes that occur at the microscopic scale contribute to the observed large-scale flow and transport distribution patterns.
Predictive modeling remains challenging for the time and spatial scales involved in geological processes and because of the lack of information about how the physical properties of the porous medium evolve as a result of chemical reactions. In particular, the role of microstructures and their possible effects on flow and transport have long been neglected. Consequently, upscaling the flow and transport properties remains poorly constrained by pore-scale observations despite a multitude of experiments, field observations, and modeling efforts.
Pore-scale processes have received increasing attention over the last decade, particularly in the context of subsurface CO2 injection and sequestration programs (Crawshaw and Boek 2013; Steefel et al. 2013). Advances in 3-D imaging techniques, such as X-ray microtomography (XMT) (Elliott and Dover 1982; Flannery et al. 1987; Sasov 1987a,b) have improved the capability for imaging porous media at the pore scale (e.g., Spanne et al. 1994; Coles et al. 1996; Lindquist et al. 2000; Noiriel et al. 2004; Bernard 2005; Blunt et al. 2013). The scale of observation from a few nanometers to centimeters can be used to link micro-scale properties of porous rocks to effective parameters at the aquifer or reservoir scale. The technique is particularly suited for 3-D rendering of structural details. In addition, the 3-D geometry derived in this way is then available for discrete calculations of porous media properties such as permeability, transport, electrical conductivity or fracture–matrix fluid transfer, and observations at the pore-scale can bring new information to improve predictability and reliability of numerical modeling (e.g., Fredrich et al. 2006; Flukiger and Bernard 2009; Noiriel et al. 2012; Bijeljic et al. 2013b; Blunt et al. 2013; Molins et al. 2014 and references below). The fact that XMT is a non-invasive, non-destructive technique allows for time-dependent observations and analyses, thus providing insights into HMC dynamic processes.
This article is concerned with the application of X-ray microtomography for the purpose of improving our understanding of dynamic geochemical processes at the pore scale. The techniques will be detailed and physical properties derived from imaging of porous rocks and fractures will be presented. The article will then examine fluid–mineral properties at the pore scale when changes of rock or material geometry are involved. Dissolution and precipitation reactions leading to particle displacement and migration are included in this analysis. Applications to dissolution and precipitation in porous and fractured rock will be presented to illustrate the capacity of XMT for investigating the effect of small-scale geometry modifications on flow and transport properties. This article is focused primarily on experimental investigations, and is not intended to cover the pore-scale modeling approaches and investigations, which are fully detailed in other articles (Liu et al. 2015; Mehmani and Balhoff 2015; Molins 2015; Steefel et al. 2015; Yoon et al. 2015, all in this volume).
Resolving pores using X-ray microtomography
The development of XMT relies on 3-D reconstruction of a sample from a series of 2-D radiographs. The technique has added a new dimension to the understanding of fluid–rock interactions in fractures and porous media through the direct visualization of pores, pore space, mineral distribution, fluid–rock and fluid–fluid interfaces, as well as fluid flow. Such pore-scale observations should be helpful in unraveling the processes that occur when a reactive fluid comes into contact with minerals.
One advantage of XMT over several other techniques is that it is non-invasive and non-destructive, and therefore can be used to characterize geometry changes during dynamic processes or experiments over time. In some cases, the imaging can be carried out in situ over the course of a dynamic experiment, in others it is conducted ex situ, before and/or after an experiment. The technique can be used as complementary to other techniques that provide additional physical, chemical, mineralogical and structural descriptions of materials, either while the experiment is running or post-mortem (e.g., measurements of fluid chemistry, permeability, acoustic velocity, electrical conductivity and resistivity, Hg-intrusion porosity, BET surface area and microscopy observations (e.g., scanning electron microscopy (SEM), focused ion beam coupled with SEM (FIB-SEM), confocal microscopy). When combined with other characterization techniques, direct comparisons can be made between the 3-D pore-structure and the different physical and chemical properties measured during experiments. The increase in computational capabilities has reduced the time needed for reconstruction, image processing and visualization of large data sets, greatly enhancing the size, resolution, and complexity of the porous materials that can be handled.
The application of XMT imaging to the field of geosciences and chemically reactive environments has recently increased thanks to the development of lab-based systems complementary to synchrotron-based systems, allowing for a micrometer resolution, and even nanometer resolution with the recent development of X-ray nanotomography (e.g., Brisard et al. 2012). XMT is based on the measurement of attenuation of an X-ray beam passing through a material. Most of the measurements are carried out using X-ray absorption properties (absorption contrast XMT, dual-energy XMT), but some specific setups also rely on coherent X-ray propagation (phase-contrast XMT) or secondary emission (X-ray fluorescence microtomography) properties. Technological evolution of instrumentation and detectors as well as the development of new reconstruction procedures continues to increase the capacities and utility of XMT. For example, the development of ultra-fast tomography allows for tracking of very rapid dynamic processes (e.g., Mokso et al. 2011; Armstrong et al. 2014a,b), as time needed to obtain a full 3-D data set can be less than 1 s. Techniques based on the energy dependence of the attenuation coefficient, such as dual-energy microtomography, can reveal the distribution of distinctive chemical compositions within a sample, as well as X-ray fluorescence microtomography (e.g., Van Geet et al. 2000; Egan et al. 2014). Detailed reviews of XMT development and applications are provided in Van Geet et al. (2001), Wildenschild et al. (2002), Carlson et al. (2003), Carlson (2006), Blunt et al. (2013), Cnudde and Boone (2013) and Wildenschild and Sheppard (2013).
Absorption contrast X-ray microtomography
When interacting with matter, X-ray incident intensity is progressively reduced. Through a homogeneous material, the attenuation is linear and depends on X-ray energy and sample atomic composition and density (in relation to its chemical composition, lattice parameter, and micro-porosity). In the energy range 1–100 keV, interactions between X-rays and matter include elastic scattering (Rayleigh scattering), inelastic scattering (Compton scattering), and the photoelectric effect. In the range 5–30 keV, within which most of the experiments are performed, the photoelectric effect is dominant.
The linear attenuation coefficient characterizes the attenuation of incident X-rays through a material, and it can be theoretically calculated for elements and mixtures, based on the NIST XCOM database (Berger et al. 2011; Fig. 1). For a single material, the total attenuation of a material over a distance x and at a given photon energy (E) follows the Beer–Lambert law:
with μTi(E) the linear attenuation coefficient of material i, and I0 and I the incident and transmitted beam intensities, respectively. Attenuation coefficients are additive over length, and a composite material crossed by an X-ray beam results, at a given photon energy (E), in a linear average attenuation coefficient, μ̄T(E):
with ϕ the length of the crossed composite material. Materials with high attenuation coefficients will absorb X-rays more efficiently, while materials with low attenuation coefficients allow more efficient transmission of X-rays.
The best experimental configuration depends on the X-ray source and beam specifications, which typically differ between a laboratory-based and a synchrotron-based setup. In a laboratory-based setup, X-rays are generated from low- or high-flux X-ray tubes resulting in a polychromatic cone beam. For low-flux sources, magnification is set by adjusting the sample position between the X-ray source and the detector, and a spatial resolution of ~1 μm can be reached. Alternatively, X-ray optics can be used with high-flux X-ray sources to achieve a higher magnification, leading to a spatial resolution as fine as 50 nm for recently developed systems.
In synchrotrons, a parallel beam of highly spatially coherent polychromatic X-rays is generated using insertion devices such as bending magnets (typical case for a X-ray imaging beamline), wigglers or undulators inserted in the storage ring. A monochromatic X-ray beam can be obtained using monochromators, which are pairs of crystals (Si(111) or multilayers) operating according to the Bragg Law of diffraction to select specific energies. Generally, the detector system operates through a scintillator converting the X-rays into visible light and an optical system (objective lenses and mirrors) ending with a detector recording the signal. Advantages in using a synchrotron source are that the data is of higher quality in terms of noise and artifact reduction (see below) and the ability to take advantage of the coherence of the X-ray beam when needed. It is also worth remarking that a pink beam (i.e., a polychromatic filtered white beam) can also be used, for example when an extremely high photon flux is needed (e.g., for ultrafast imaging). A synchrotron source may be preferred to a polychromatic cone beam, as the higher the signal-to-noise ratio, the more accurate and precise is data segmentation and quantification.
Image acquisition consists of recording a series of 2-D radiographic projections of the transmitted photons at different angular positions of the sample placed on a mechanical stage rotating over 180° (parallel beam) or 360° (cone beam) around an axis perpendicular to the beam. In some setups, the sample is fixed and the X-ray source describes a helical trajectory relative to the sample (spiral cone beam). In absorption contrast tomography, the sample is set immediately in the front of the scintillator to avoid phase contrast (see further discussion below). Between one and three thousand radiographs are usually recorded. The X-ray beam energy is adjusted to a value which allows for a high absorption contrast between the different materials forming the composite sample (e.g., mineral and air), while a reasonable transmission of X-rays (~20–30%) in the most absorbing zone of the sample is still available (to minimize noise). The exposure time can be increased to increase the transmission of X-rays, but should remain short enough to avoid saturation of the detectors when bright field images (i.e., radiographs without the sample in the field of view) are acquired. Phase-contrast imaging or holotomography can substitute for absorption modes when contrast between two phases of interest is not high enough. Phase contrast in the near-field region can also be used to highlight small features at or just below the resolution of the experimental setup, such as grain boundaries or cracks (e.g., Marinoni et al. 2009). For samples rich in heavy elements, higher energies are generally required, as increasing energy increases X-ray penetration. On a synchrotron source, the brightness (i.e., the photon flux) of a beam derived from an insertion device decreases at the higher energies, meaning that samples rich in heavy elements will generally require a longer acquisition time. This effect is in addition to the fact that the exposure time must also generally be increased for heavy elements, as X-ray linear attenuation is higher.
Radiographs (i.e., projections) of approximately 2000 × 2000 pixels are recorded. The intensity of transmitted photons is detected using a charge coupled device (CCD) or a flat panel detector coupled to scintillators converting X-rays to visible light. Before reconstruction (the procedure converting the set of radiographs in a 3-D volume), the radiographs must be corrected for background noise, beam structure, defects in the detector, and beam intensity fluctuation during acquisition. Fluctuations in beam intensity are usually corrected by normalizing the projections over a region of interest (ROI) where the sample is not present. Other corrections are carried out by subtracting dark field and flat field from original radiographs, using the following formula:
with I0cor and Icor the corrected intensities of the incident and transmitted beam, respectively, Idark the intensity of the dark field image, which is a projection collected without the beam (to record the background noise), and Iflat the intensity of the flat field, which is a radiograph without the sample in the field of view.
Image reconstruction and correction
3-D image reconstruction producing data sets of approximately 2000 × 2000 × 2000 voxels (i.e., volumetric pixels) can be achieved using either analytic (i.e., filtered back projection) or iterative (i.e., algebraic or simultaneous algebraic) reconstruction methods. The most popular is the filtered back-projection method based on the inverse Radon transform (Herman 1980). The set of projections over all the angular positions is first converted line by line in the frequency domain using the Radon transform. The set of projections for a specific line over all the angular positions in the frequency domain is called a sinogram. 3-D reconstruction, which corresponds to the Cartesian coordinate back transformation in the spatial domain, is then performed using a filtered back-projection algorithm. The reconstruction is performed slice by slice, and the 3-D volume is obtained by stacking the 2-D slices.
Several artifacts can result from XMT acquisition, depending on the beam characteristics, setup and sample composition, such as hot pixels, streak artifacts, ring artifacts, beam hardening and cone-beam effect. They appear in the 3-D reconstructed data sets as either a significant number of distinct geometrical features (e.g., lines, streaks, rings) or globally spread noise. A number of algorithms exist for improving the quality of 3-D reconstructed data sets (Fig. 2). Presently, most of the corrections are made routinely before, during or after reconstruction, and are included in reconstruction software (e.g., Dierick et al. 2004; Vlassenbroeck et al. 2007) without any end-user intervention.
Observations of geometry changes at different stages of a dynamic process require 3-D image registration of the data sets in the same coordinate system in order to map a 3-D volume to a reference one. Digital image correlation is often used to quantify local strain distribution when mechanical deformation of the sample has occurred (e.g., Lenoir et al. 2007; Hall et al. 2010). When there is no mechanical deformation, registration consists of a rigid transformation involving 3 rotations and 3 translations around and in the x-, y- and z-directions, combined with interpolation. Many optimization algorithms based on image intensity or feature extraction can be used to perform 3-D image registration in either spatial or frequency domains as, for example, after the calculation of the centre of gravity of intensities, the moment of inertia in the x-, y- or z- direction, or normalized cross-correlation (e.g., Maintz and Viergever 1998). However, when studying dynamic processes, automatic registration using iterative optimization algorithms often fails because geometry changes lead locally to modifications of features that change the image intensity distribution. In particular, it can be difficult to establish a correspondence between several distinct features. If this is the case, a manual registration remains the only alternative. Manual registration is based on the selection of several control points (or landmarks), which generally correspond to isolated pores or very absorbing minerals (e.g., iron oxides or sulfides). A minimum of four points are required to calculate the six unknown parameters of the rigid-transformation matrix, but a statistically representative collection of control points is preferred for error minimization. Different interpolation techniques can be then used to recalculate the grayscale of every pixel in the registered data sets, e.g., tri-linear interpolation (Gonzales and Woods 1992).
Phase contrast X-ray microtomography
The absorption contrast between materials of different densities (e.g., air and mineral) is generally large enough for them to be distinguished on 3-D data sets. However, XMT in the absorption contrast mode can fail to distinguish between different materials (or phases) when the absorption contrast is poor (e.g., between liquid and vapor), or when their densities are similar (e.g., water and organic material, calcite and dolomite, or quartz and feldspar). In these cases, propagation-based imaging, including phase contrast microtomography, can be used to enhance contrast at the interface of two phases (e.g., Snigirev et al. 1995).
Phase-contrast imaging relies on the fact that two isotropic materials of different density have a different refractive index, i.e., different X-ray propagation speeds. As a result of the difference in density, not only the intensity but also the phase of X-rays is modified when interacting with different materials, leading to a phase shift. Weak perturbation of the wave front generates interference patterns. The interference fringes are not directly linked to the phase itself, but rather to the second derivative of the phase, with the result that the method is more sensitive to abrupt changes in the decrement of the refractive index. Near-field propagation phase imaging enhances phase contrast in the regions of large changes in the refractive index between two phases, i.e., at their edges. For low-absorption materials or when absorption is similar between two phases, an alternative is to increase the sample–detector distance to allow for recording of the interference patterns (Cloetens et al. 1999). In this case, a high-resolution detector is required, although it limits spatial resolution. Retrieval of the phase distribution of the sample using specific algorithms can also be useful to obtain the 3-D distribution of the complex refraction index and thereby enhance the contrast between phases (e.g., Arzilli et al. 2015).
X-ray fluorescence microtomography (XFMT)
The photoelectric effect involves atom ionization (i.e., ejection of an electron) resulting from the interaction between X-rays and the inner shells of an atom. Filling the inner-shell vacancy in the atom is accompanied by either the emission of a second electron (Auger effect) or by emission of a secondary fluorescence X-ray photon (X-ray fluorescence). Fluorescence is emitted isotropically from all excited atoms, and the use of an energy-dispersive detector (e.g., a Si-drift diode detector) allows for recording integrated 2-D projections of the elemental distributions within the sample. The reconstruction of the 3-D elemental distributions is then possible from the 2-D projections. Combining XMT with XFMT means that both geometric and chemical mapping can be achieved. The technique is useful for distinguishing between minerals with similar densities (i.e., silicates, solid solutions), measuring element concentration ratios (Lemelle et al. 2004), and analyzing the distribution of radionuclides and metals in soils (Lind et al. 2013). However, the spatial resolution is limited (micrometer resolution) and acquisition time can take several days. Tracking of several ions (Cs, Ba, La) was also performed during diffusion experiments in porous media (Betson et al. 2005) and the technique appears promising for environmental studies, i.e., for quantifying trace elements in samples, although there are several technical challenges and analytical complexities (e.g., the absorption of the fluorescence X-ray photons by the sample itself) with its use.
Dual-energy X-ray microtomography
An XMT data set is a 3-D image of the X-ray attenuation of the sample. Dual-energy XMT is based on the dependence of the linear attenuation coefficient of an element or a material to the X-ray energy (E). The technique requires scanning the sample twice at two different energies. Van Geet et al. (2000) and Remeysen and Swennen (2010) used the technique with a Skyscan 1072 lab-source at 100 and 130 kV to recover densities and effective atomic numbers within carbonate samples. Nevertheless, due to the noise inherent in a polychromatic source, segmentation was required in order to calculate the repartitioning of the different rock-forming minerals.
Due to the ability to select specific energies, synchrotron source provides the best experimental setup for examining a specific element with dual-energy XMT, as a large difference in absorption exists in the absorption edge region of an element (see Fig. 1). In this case, the acquisition is performed at two energy settings, just above and below the absorption edge energy of the considered element, a range over which the attenuation coefficients of the other elements do not significantly change. The difference between the two XMT data sets therefore allows a 3-D chemical mapping of the considered element, which can in turn be used to obtain 3-D mineral distribution in materials such as rock specimens (Tsuchiyama et al. 2013). The technique is suitable for mapping elements with K- or L-absorption edges in the same energy range as that set up to image rock samples (e.g., Fe, As, Pb, Cs, U).
Image segmentation, in which voxels are classified into two (i.e., pore and solid) or more phases of interest (e.g., the different constituents of the solid phase, such as minerals), is an essential step for extracting geometrical properties of rocks. Digitization of the pore space and mineral matrix is also required for performing direct numerical simulations of flow and transport (e.g., Oren et al. 2007; Flukiger and Bernard 2009; Narsilio et al. 2009; Menon et al. 2010; Bijeljic et al. 2013a; Molins et al. 2014). A 3-D XMT data set of rock is assumed to consist of a fixed number of distinct mineral phases (solid phase) and fluid and/or gas phases (pore space). The dynamic range (i.e., the distribution of grayscale data) of several data sets can vary for different acquisitions of the same sample due to fluctuations in the X-ray beam or additional noise, even if similar acquisition and reconstruction parameters were chosen. As a consequence, normalization of every data set of a sample is necessary before segmentation.
Segmentation involves converting the grayscale data into several sets of voxels corresponding to the distinct phases of interest based on their different X-ray absorption properties. In most cases, the different phases cannot be clearly separated on the grayscale histograms as the data is noisy and can be complicated by partial volume effects (see further explanation below), with the result that applying a simple threshold method to the whole image would lead to misclassification. Different methods have been developed to improve segmentation accuracy (see Iassonov et al. 2009 and Schlüter et al. 2014 for reviews), although there is no clear test to evaluate their performance. For example, edge-preserving noise reduction filters, such as a 3-D median filter or an anisotropic diffusion filter, can be used before segmentation to improve the signal-to-noise ratio without altering information at the interface between two phases. Methods that incorporate local spatial information (i.e., based on the intensity of the neighboring voxels) are usually more appropriate over global threshold methods, e.g., indicator kriging (Oh and Lindquist 1999), watershed transform (Beucher 1992; Sheppard et al. 2004), or region growing (Pitas 2000). After segmentation, each phase is assigned a fixed value (e.g., 0, 1, 2, etc.). In some cases, segmentation can also produce image artifacts, even when applied to high quality XMT data sets, and filtering of the segmented data can be necessary.
In addition to the algorithm used, segmentation efficiency depends on the initial image quality (i.e., signal-to-noise ratio, which directly impacts the distribution of the grayscale values of the different phases on the histogram), and on whether the image resolution is sufficient to resolve the objects of interest. For example, a porosity difference of less than 0.5% was obtained by XMT consistent with chemical measurement (Noiriel et al. 2005). Uncertainties, however, are usually higher, particularly when the pore network is complex or there is a insufficient resolution in the image. A comparison between different segmentation methods by Andrä et al. (2013a) gave a maximum porosity difference of about 2% for Berea sandstone. In contrast, the difference was approximately 8% for Grosmont dolostone where the histogram of grayscale values was much more complex.
In some cases, contrast agents can be added to enhance contrast of the fluid phase or to study fluid partitioning in the pore space (e.g., Seright et al. 2002). Nevertheless, the presence of a contrast agent of intermediate grayscale between initial fluid and solid can affect how the other phases can be segmented (Prodanovic et al. 2006, 2007). It should be mentioned that segmentation is a critical step in image processing, as the quality of segmentation directly impacts the computed effective properties (Andrä et al. 2013b).
Porosity determination, pore geometry and pore-space distribution
Once segmentation has been completed, porosity, φ, can be calculated as simply the ratio of pore space to the total number of voxels, i.e., φ = Nvoid/(Nvoid + Nsolid), with Nvoid and Nsolid the number of voxels identified as pore and solid, respectively. As reactivity and other properties may change along the flow path, it may be important to calculate porosity as a function of sample depth (z-direction) or sub-sampled in every direction of space. When partial volume effects are too large to resolve pores, it may still be possible to estimate porosity without segmenting the data. In this case, either references of attenuation of the different materials forming the sample (e.g., μsolid the attenuation of the solid and μvoid the attenuation of air) or a porosity reference (e.g., determined from Hg-intrusion) are required (Fig. 3). Porosity of a voxel, in which attenuation is μ, is defined as:
Several parameters can also be extracted from 3-D images, for example to obtain an idealized representation of the complex pore space geometry and a description of topology (see Thovert et al. 1993; Spanne et al. 1994). Extraction of the pore-space structure is commonly based on medial-axis skeletonization, which can be used to partition the void space into pore-bodies and pore-throats. Lindquist et al. (1996, 2000) have developed algorithms to calculate pore connectivity, tortuosity, pore-body size, pore-throat size and surface area based on studies of the Berea and Fontainebleau sandstones. Al-Raoush and Willson (2005) calculated the spatial correlation between pore-body sizes for several unconsolidated porous media using semi-variograms and integral scale concepts. Bauer et al. (2012) developed a dual-pore skeletonization to obtain a separate description of macro-porosity and micro-porosity in limestone samples with bimodal pore-size distribution. Their model combines electrical transport properties of micro-porosity with those of the interconnected macro-pore network. Skeletonization is, however, sensitive to image resolution and noise, as small details or topological artifacts (e.g., edge-connected voxels of two non-connected objects) can lead to errors when estimating the medial axis. Increasing resolution does resolve this problem, but Plougonven and Bernard (2011) proposed a method based on the sequential detection of these topological artifacts associated with their removal to improve the reliability of skeletonization without adding resolution.
To overcome partial volume effects and solve for non-Fickian dispersive transport in porous limestone, Gouze et al. (2008) developed an alternative approach based on a combination of superimposed grayscale and segmented data sets. First, they combined a burning algorithm with a down-gradient search in order to obtain a dual-pore description based on the mobile–immobile fraction of connected porosity. Combined with the non-connected porous network (considering a voxel resolution of 5.01 μm) and the micro-porous matrix, three different domains can be identified: one is a mobile fraction composed of connected porosity where advection is dominant; the second is an immobile fraction made up of connected porosity in which flow is negligible and has both unconnected porosity and a micro-porous matrix, the latter of which is accessible for diffusion; and a third consisting of a solid fraction below which a critical porosity threshold for diffusion is effectively zero (Fig. 4).
Solid-phase distribution and quantification
The identification and characterization of the different minerals or materials forming porous media can be carried out in a similar manner to the characterization of pore space after segmentation. Smith et al. (2013b), for example, combined XMT with SEM analysis to quantify calcite and dolomite abundance in a vuggy limestone and marly dolostone based on grayscale intensities. They quantified the reactivity of both minerals during dissolution over the entire core samples as well as within the regions where wormholes developed later in the experiment.
Labeling of individual materials implies that they can be distinguished in the 3-D images, either based on their attenuation properties or on their morphology. Morphological criteria can be used to distinguish between materials, particularly when three-phase or more segmentation is not possible. However, the technique is often restricted to idealized materials with regular shapes, e.g., glass beads or mineral aggregates. Starting with binary images, the centroids of the different grains can be found by computing a distance map, the maxima regions of the distance maps corresponding to the centroids. Once every centroid has been identified, a fast watershed algorithm based on topographic distance can be applied, inverting the distance map in order to find out the separation lines that correspond to the boundaries between grains. However, it is worth mentioning that the efficiency of watershed separation depends on the degree of consolidation of the medium and the shape of individual grains, and partitioning of clusters into individual grains is often not trivial (Thompson et al. 2006). Once all the grains have been separated, they can be labeled. Morphological criteria, such as grain size or sphericity, can be used to classify different grain populations (e.g., Dann et al. 2011). Figure 5 illustrates the procedure, with a comparison with three-phase segmentation and labeling (Fig. 6).
Once the different grains have been labeled, it is possible to calculate individual and statistical parameters based on their morphology. Fonseca et al. (2012), for example, characterized the morphology of the grains forming Reigate sand, including the number of particles, particle size and surface area, orientations of the major, minor, and intermediate axes for each particle, and the distributions of particle aspect ratio, sphericity, and convexity. The authors observed that calculations of particle size and shape made by analysis of 3-D images differed appreciably from the values obtained from 2-D images.
Specific and reactive surface area measurements
Rates of fluid–mineral interactions are directly controlled by the mineral surface area available for sorption and dissolution/precipitation. The measurement of the surface area (Sr), which is often not simply related to reactive surface area, is essential for interpreting experimental results and for numerical modeling of reactive transport. Reactive surface area is generally determined from BET measurement (Brunauer et al. 1938) or a geometrical estimate of the specific surface area. However, imaged material suffers from resolution limitation compared to BET measurements, and different methods generally lead to different estimates of the geometric surface area (Seth and Morrow 2007). Estimation of geometric surface area is possible from grayscale XMT data sets using a marching cube algorithm (Lorensen and Cline 1987), but is mostly carried out using segmented data sets. At some point, improving resolution does not improve the estimate (Dalla et al. 2002; Porter and Wildenschild 2010). Although geometric surface area scales directly with the observation scale, XMT is a powerful tool for resolving surface area changes at the pore scale, particularly when calculations are performed on the same sample at different stages of experiment. For example, by superimposing XMT volumes (Fig. 7), Noiriel et al. (2005) were able to locally quantify and distinguish between reactive and non-reactive surface areas during dissolution of a porous limestone sample. The shifting (and non-shifting) of the fluid–mineral interface was directly linked to reactive (and non-reactive) mineral surfaces. XMT can also be combined with 2-D microscopic analyses for improved estimates of the contribution of each mineral to the total specific surface of the rock (Landrot et al. 2012; Smith et al. 2013a).
Geometrical estimates of the specific surface area changes were also carried out using XMT during dissolution (Noiriel et al. 2009; Gouze and Luquot 2011; Qajar et al. 2013; Molins et al. 2014) and precipitation (Noiriel et al. 2012) experiments. Although XMT makes it possible to determine specific surface area, a poor correlation between reactive and specific surface areas is often reported (see, for example, Anbeek 1992, who reports a larger specific surface for weathered minerals compared to unweathered minerals). The weathered surfaces actually correspond to low-reactivity surfaces, similar to what is observed when different dissolution rates between minerals contribute to the increase of specific surface area (e.g., Gouze et al. 2003). To complicate matters further, dissolution–precipitation reactions may result in changes of reactive surface area in addition to the effects on specific surface area (Noiriel et al. 2009; Gouze and Luquot 2011). Unfortunately, the changes are not necessarily correlated. For example, Noiriel et al. (2009) measured the changes in reactive surface (from chemistry) and specific surface (from XMT images) areas, and observed different trends during a dissolution experiment in a porous limesto ne containing two different calcites (micrite pellets and sparite cement), with different chemical signatures, particularly in Sr and Ba. However, their geometric model of spherical pore growth and micrite sphere reduction failed to reproduce observations of fluid chemistry. A semi-empirical sugar-lump model based on observations carried out at a finer scale using SEM was more successful in describing the increase of sparite crystal boundaries exposed to the fluid in this example.
Flow in single fractures is closely related to aperture distribution, fracture wall roughness, tortuosity and asperities. Peyton et al. (1992), Jones et al. (1993), Keller (1998), and Vandersteen et al. (2003) used different estimates based on grayscale levels to define and measure fracture aperture, am, from XMT images with partial volume effects. Gouze et al. (2003) used segmented data to obtain the void structure, then extracted fracture walls and aperture. However, their method has some limitations, particularly when altered fracture walls exhibit some overlap or when secondary fracture branching occurs. Characterization of the void structure (Fig. 8) permits computation of aperture distribution and related statistics, e.g., fracture volume, tortuosity, and contact areas (e.g., Gouze et al. 2003; Karpyn et al. 2007; Noiriel et al. 2013). Methods of describing fracture wall roughness can also be applied, e.g., determination of the roughness factor (Patir and Cheng 1978), roughness coefficient (Myers 1962), or fractal dimension (Brown and Scholz 1985). Either the 3-D void structure or the 2-D maps of aperture or fracture wall elevation can be used as an input for flow modeling (e.g., Noiriel et al. 2007a, 2013; Crandall et al. 2010a,b). Crandall et al. (2010a) generated different meshes from the same data set on a fractured Berea sandstone core, by changing the fracture wall roughness properties through different refinement procedures. Given the different results for tortuosity and transmissivity, Crandall et al. (2010a) established a relationship between the method by which a fracture is meshed and the flow properties calculated by numerical models.
Rocks typically have complex structures over a wide range of length scales. XMT is unable to resolve pores at a scale smaller than the spatial resolution, with the result that there is insufficient resolution to capture the complete geometry of features, leading to the so-called partial volume effects. When this occurs, voxel intensity balances intensities of several phases (e.g., water and solid matrix at the fluid–rock interface), and this can result to misleading representations of the pore size and connectivity. In contrast, a higher resolution is better for describing the fluid–rock interface, pore shape and size, but it can fail to provide a reasonable estimate of the size of the representative elementary volume (REV), which is the smallest volume which allows for a continuum description of a representative property, for example permeability. In addition, the resolution achievable in XMT is not only limited by the beamline and detector properties, but also by the sample size, as the sample must usually remain within the field of view during data acquisition (although local XMT is possible, see below). Although mosaic scanning can overcome part of the problem, the sample size can remain a limiting factor for investigating porous media, particularly when the pore structure is heterogeneous and a REV is sought for upscaling the flow and transport properties. In other words, although scales much smaller than the resolution used for imagery can play an important role in flow and transport, these features cannot be seen, leading to erroneous parameter estimation or poor process modeling in reactive transport. The most sensitive parameter is certainly the mineral surface area, the measurement of which depends on scale, and thus on the ability to resolve micro-porosity. Permeability can also be poorly estimated, especially when its value is constrained primarily by small pore throats or narrow channels. To overcome this problem, Prodanovic et al. (2015) derived a two-scale network from their XMT data sets to estimate flow and transport properties: a macro-network that mapped the inter-granular porosity that was clearly resolved with XMT; and a micro-network that mapped the micro-porous regions unresolved with XMT, the properties of which could be derived from SEM observations at a higher resolution. Other studies include XMT investigations coupled with SEM or FIB-SEM analyses (e.g., Sok et al. 2010; Landrot et al. 2012; Beckingham et al. 2013).
Multi-resolution XMT, which consists of data acquisition at various resolutions, can be used to obtain a detailed description over several different scales. It sometimes involves local XMT, which is a technique used when the sample is larger than the field of view and requires a specific reconstruction method of the region of interest to compensate for the incomplete projection data. Several studies have examined the multi-scale characterization of rock samples (Bera et al. 2011; Dann et al. 2011; Peng et al. 2012, 2014; Hébert et al. 2015). Hébert et al. (2015) have investigated the intrinsic variability and hierarchy of the connected pore space of limestone and dolostone samples at different XMT resolutions ranging from 0.42 to 190 μm. The distribution of porosity values between 0.42, 1.06, and 2.12 μm resolutions highlights partial volume effects (i.e., more porosity details are visible when resolution is increased, Fig. 9), and challenges the ability of XMT to provide reliable representations of pore networks and permeability estimates for low-porosity and highly heterogeneous samples. However, Peng et al. (2014) observed that the contribution of small pores to the permeability was minor in Berea Sandstone, where the pore network appears to be well connected at a resolution higher than 5.92 μm. However, the micro-porous phase that was resolved at a resolution of 1.85 μm (but not at 5.92 μm) was shown to increase porosity, surface area, and pore network connectivity estimates (Fig. 10).
COMBINING EXPERIMENTS, 3-D IMAGING AND NUMERICAL MODELING
The first investigations of 3-D fractures and porous media at the micrometer scale combined with experimentation in the geosciences were focused on single-phase flow experiments and modeling of macroscopic flow and transport properties (Coles et al. 1996; Arns et al. 2001, 2005; Enzmann et al. 2004a,b; Fredrich et al. 2006; Bijeljic et al. 2013b; Kang et al. 2014). Colloidal deposits and their effects on permeability reduction were also investigated (Gaillard et al. 2007; Chen et al. 2008, 2010). Studies of multi-phase flow and gas trapping during imbibition or drainage experiments for the purpose of enhancing oil recovery or for better understanding of groundwater contamination by non-aqueous phase liquid (NAPL) or CO2-trapping processes have also benefited from XMT (e.g., Coles et al. 1996; Seright et al. 2002; Wildenschild et al. 2002; Arns et al. 2003; Culligan et al. 2004; Al-Raoush and Willson 2005; Karpyn et al. 2007; Karpyn and Piri 2007; Prodanovic et al. 2007, 2010; Al-Raoush 2009, 2014; Iglauer et al. 2011; Silin et al. 2011; Armstrong et al. 2012; Blunt et al. 2013; Ghosh and Tick 2013; Kneafsey et al. 2013; Smith et al. 2013a,b; Andrew et al. 2014; Celauro et al. 2014; Landry et al. 2014; see also Wildenschild and Sheppard 2013 for a review). Several studies have also examined multi-phase flow and transport in soil aggregates (Carminati et al. 2008; Koestel and Larsbo 2014), including evaporation and salt crystallization in the form of discrete efflorescence or desiccation cracks (Shokri et al. 2009; Rad et al. 2013; DeCarlo and Shokri 2014).
X-ray microtomography for monitoring reactive transport
Reactive transport studies include dissolution or precipitation experiments, biofilm growth, or weathering processes. Most of the dynamic experiments performed in the lab are based on flow-through reactors (e.g., Noiriel et al. 2012; Kaszuba et al. 2013). The apparatus is generally composed of cm-scale columns packed with solid material (e.g., glass beads, crushed crystals, or rock aggregates), or core holders containing rock samples. In some cases, the experimental setup can be mounted directly on the X-ray beamline. However, the core holder must be made of material transparent or weakly attenuating to X-rays, or the samples must be removed from the core holder before imaging. Some studies have used XMT to visualize changes in pore structure, while other parameters can be recorded to better clarify the evolution of flow and transport properties and processes, e.g., permeability measurement, chemistry analyses, pH, or tracer experiments. Further measurements can also be made at different stages of the experiment (e.g., Vialle et al. 2014), including geophysical monitoring (e.g., P-wave velocity, electrical conductivity), BET surface area measurement, or determination of the porosity distribution by Hg-intrusion.
XMT images can also be used directly or indirectly to estimate solute transport properties in low-porosity materials or to follow the propagation of fronts in reactive transport experiments. Polak et al. (2003), Altman et al. (2005), and Agbogun et al. (2013) have used X-ray absorbent tracers (NaI, CsCl, or KI) to directly visualize their displacement within the porous network. Nakashima and Nakano (2012) combined 3-D image analysis of porosity and surface-to-volume ratio with a tracer experiment (KI) to determine tortuosity. Mason et al. (2014) developed a method that relies on the linear attenuation coefficients to estimate the volume percentage of different materials (i.e., carbonates, amorphous components, and porosity) and the Ca-concentration profiles in cement during alteration by a CO2-rich brine. Burlion et al. (2006) followed the alteration front displacement in a mortar during an accelerated leaching process by an ammonium nitrate solution (Fig. 11). They observed an exponential decrease in the rate of advance of the propagation front through the altered cement paste over time, while aggregates remained unaltered. The experimental results can be compared with a theoretical or numerical solution of the diffusion or reaction–diffusion equation. Noiriel et al. (2007b) also observed that the rate of propagation of a calcite dissolution front across fracture walls in an argillaceous limestone decreased exponentially due to modification of transport mechanisms. Initially, the transport mechanisms were advection-dominant within the fracture. Then, they evolved to a combination of advection in the fracture and diffusion in the newly formed micro-jporous clay coating which grows over time.
The propagation rate of a reaction front in diffusion-dominated systems does not necessarily result in a parabolic (t1/2) dependence on time, particularly when feedback between porosity and the effective diffusion coefficient is involved. Navarre-Sitchler et al. (2009) quantified porosity changes associated with weathering processes in a basalt clast from Costa Rica as a result of mineral dissolution and precipitation of secondary phases. They combined 3-D image analysis with laboratory and numerical diffusion experiments to examine changes in total and effective porosity and effective diffusion coefficient across the weathering interface. Due to an increase in the porosity as a result of chemical weathering, the diffusive transport of aqueous weathering products away from the core/rind interface was enhanced (Navarre-Sitchler et al. 2015). This occurs at a critical porosity of about 9%, beyond which the number of connected pathways for diffusive transport increases dramatically. Reactive transport modeling further demonstrates that the rate of advance of the weathering front can be constant over time, (i.e., linear) even if the dominant transport mechanism is diffusion, in the case where an increase of the effective diffusion coefficient occurs due to porosity enhancement (Navarre-Sitchler et al. 2011).
A number of studies of porosity and permeability evolution associated with injection or sequestration of CO2-rich fluids in reservoir rocks, caprocks, or cements have been carried out (e.g., Gouze et al. 2003; Noiriel et al. 2004, 2005; Luquot and Gouze 2009; Gouze and Luquot 2011; Luquot et al. 2012, 2013, 2014a,b; Jobard et al. 2013; Smith et al. 2013a; Ellis et al. 2011, 2013; Abdoulghafour et al. 2013; Deng et al. 2013; Jung et al. 2013; Luhmann et al. 2013; Mason et al. 2014; Walsh et al. 2014a,b). Complementary to the experiments, numerical modeling has also been used to evaluate mass transfer at the pore scale (e.g., Flukiger and Bernard 2009; Molins et al. 2014). Direct simulation on 3-D segmented images appears to be the modeling approach of choice for single-phase flow and transport, since the complexity of the pore space is preserved with this approach (Blunt et al. 2013).
Most experiments involving reactive transport have focused primarily on mineral dissolution, but mineral precipitation experiments are also reported in the literature, under either abiotic or biotic conditions. Precipitation experiments involving biofilm growth (Davit et al. 2010; Iltis et al. 2011) or biomineralization (Armstrong and Ajo-Franklin 2011; Wu et al. 2011) were conducted to improve in situ degradation of organic pollutants or sequestration of radionuclide or trace metals into solid phases. Noiriel et al. (2012) evaluated upscaling of precipitation rates from an integrated experiment and modeling approach to the study of calcite precipitation in columns packed with glass beads and calcite spar crystals. Although upscaling was possible using kinetic data determined from well-stirred reactor experiments, the study highlights that nonlinear, time dependence of reaction rates, as related to evolving surface area and/or reactivity, can be difficult to assess in natural contexts.
In other cases, a more complex suit of reactions are involved, including dissolution of primary minerals by the reactive fluid that leads to precipitation of secondary phases. Cai et al. (2009) and Crandell et al. (2012) quantified the changes in porosity and pore size characteristics in Hanford sediments exposed to simulated caustic tank waste. The authors found that secondary precipitation of sodalite and cancrinite, following the dissolution of quartz and alumino-silicates, resulted in a decrease in both the abundance and size of large pore throats.
X-ray microtomography for monitoring geomechanical evolution
Evaluation of mechanical or mechano-chemical properties of rocks, concrete, and aggregates have also benefited from XMT observations. Dautriat et al. (2009) measured the dependence of porosity and permeability on the stress path in elastic, brittle, and compaction regimes in Estaillades limestone samples. The permeability drop was linked to pore collapse associated with the propagation of micro-cracks through the dense calcite (i.e., bioclasts and cement) and filling of the pores by some fine calcite fragments. In this case, the less dense calcite (i.e., micro-porous red algal debris and altered micro-sparitic cement) appeared to accommodate strain in a more diffusive manner. By comparing the digital images of the specimen in the reference and the deformed states, Higo et al. (2013) were able to apply digital image correlation to obtain the full-field surface displacements in Toyoura sand during triaxial tests. They were able to observe the rotation of sand particles and their displacement in shear bands as a result of the expansion of the voids associated with dilation. Cilona et al. (2014) investigated the effects of rock heterogeneity on the localization of compaction in porous carbonates. They performed a micro-structural characterization of deformation bands under different strain and stress conditions and determined the crack distribution and density. Compaction bands were found to be localized in the laminae where porosity was initially higher.
Renard et al. (2004) performed a series of uniaxial stress-driven dissolution experiments on halite aggregates for the purpose of studying their deformation as a result of pressure solution. They measured directly the vertical shortening of the sample on the radiographs, determining that the reduction in permeability from 2.1 to 0.15 Darcy (after 18.2% compaction) was linked to grain indentation and pore connectivity reduction by precipitation on the free surface of pore throats. Peysson et al. (2011) investigated permeability alteration due to salt precipitation during drying of brine, and reported the accumulation of salt near the sample, which led to a reduction in permeability. In this study, complete pore sealing due to the precipitation of salt resulted in very low permeability. In contrast, a drying experiment conducted by Noiriel et al. (2010), in which intense fracturing was induced by halite precipitation, involved an increase in permeability.
Rougelot et al. (2009) were able to link micro-crack formation in a cementitious material containing 35% glass beads to the leaching of calcium from the cement paste by ammonium nitrate on the basis of the observation of a higher density of micro-cracks around the glass beads. Using numerical simulations complementary to their observations, the authors showed that cracking was inherent to the initial pre-stressing of the composite, inducing tensile stress to develop around glass beads as the mechanical properties of the leached cement paste deteriorated. Dewanckele et al. (2012) quantified porosity changes in a building rock sample during weathering processes by SO2. Despite a reduction in pore size, the authors observed an increase in porosity as a result of the formation of micro-cracks in the rock caused by efflorescence.
Effect of mineral reaction kinetics on evolution of the physical pore space
The dissolution and precipitation kinetics of minerals, combined with the effect of transport of species to and from the fluid–mineral interface, exerts an important control on the processes of mass transfer between fluids and minerals along a flow path. The relative importance of surface versus transport control during dissolution/precipitation has been extensively discussed (e.g., Plummer et al. 1978; Berner 1980; Rickard and Sjöberg 1983; Brantley 2008; Steefel 2008). The overall reaction rate integrates both the effects of geochemical reactions and transport of species, which is difficult to access in transport-controlled conditions, as transport processes can be highly heterogeneous in rocks.
The transition state theory (TST) is generally considered appropriate to describe geochemical reaction rates (Lasaga and Kirkpatrick 1981). A general equation to describe dissolution/precipitation reaction rates of a mineral by TST is given by:
with kr and kr′ the reaction and intrinsic kinetic rate constants, respectively (mol·m−2·s−1), Sr the reactive surface area, ai the activity of the inhibitor and catalyst species i of the reaction, and Ω the saturation index (Ω = IAP/Ksp), with Ksp the solubility and IAP the ion activity product. j, n and n′ are semi-empirical constants the values of which depend on the kinetic behavior involved in the chemical reaction. Rates of mineral–water interaction have been conventionally determined from bulk measurement during well-stirred flow-through experiments, after determination of the mineral surface area using the BET method; this assumes that the surface area SBET scales linearly with the reactive surface area Sr. However, the generally observed increase in geometric surface area that occurs without a corresponding change of reactivity (Gautier et al. 2001), along with the large discrepancies in observed roughness factors between fresh and weathered minerals (Anbeek 1992), suggests that (i) the density of reactive sites is poorly correlated with geometric surface area and (ii) the reactive surface may change during reactions as observed in numerous studies (Arvidson et al. 2003; Lüttge et al. 2003; Hodson 2006; Noiriel et al. 2009).
More recently, the determination of bulk rates of dissolution/precipitation have been augmented by studies of surface processes at the microscopic scale using various techniques such as atomic force microscopy (AFM) (e.g., Hillner et al. 1992; Stipp et al. 1994; Dove and Platt 1996; Shiraki et al. 2000; Teng et al. 2000; Teng 2004) and optical interferometry (e.g., Sjöberg and Rickard 1985; Lüttge et al. 1999; Arvidson et al. 2003; Fischer and Luttge 2007; Colombani 2008; Cama et al. 2010; King et al. 2014). With these characterization techniques, experiments can be performed under continuous flow conditions and the dynamics of the mineral surface can be resolved with a lateral and vertical resolution of ~50 nm and ~0.1 nm, respectively. Although they have also been applied to the polished surfaces of fine-grained rocks (e.g., Emmanuel 2014; Levenson and Emmanuel 2014), these methods are mostly limited to the study of micrometer-scale oriented crystal faces due to the limited depth of field possible. XMT can be used in complement to these methods to provide full 3-D measurements of precipitation or dissolution rates. In addition, coupling flow-through experiments with XMT observation gives access to overall reactions rates that integrate both chemical reactions and transport effects over time and space.
Rates of dissolution/precipitation reactions in porous rocks
As a result of mineral dissolution/precipitation associated with reactive fluid injection or mixing, the fluid–rock interface can evolve over time, thus modifying the flow and transport properties of rocks. The velocity of the moving interface is generally non-uniform through space and time as a result of the difference in reaction rates between minerals and the modification of transport processes at the pore scale. XMT, after segmentation and registration of 3-D data sets, offers an alternative method for calculating rates in such dynamic systems at the pore scale. Rates integrate the effects of both chemical reactions and local transport conditions; they can be measured through time-lapse positioning of the fluid–solid interface. Mapping of the displacement velocity of the fluid–mineral (or fluid–solid) interface is possible from image subtraction. The surface-normal velocity (m·s−1) is defined by:
with Ifs the vector position of the fluid–solid interface and n the normal to the surface. By combining a distance transform (e.g., Pitas 2000) of the original data set with the image difference, it possible to evaluate the displacement of the interface over time and thus calculate the velocity. Integration of the velocity over a surface allows for calculation of the local dissolution/precipitation rate, r (mol·s−1):
with S the surface area (m2), V the local volume dissolved/precipitated (m3) and V the mineral molar volume (m3·mol−1). The kinetic rate kr (mol·m−2·s−1) can also be calcul ated, assuming that the reactive surface area is known, i.e., kr = r/Sr.
An example of a pore-scale rate calculation is shown in Figure 12, which shows the calcite precipitation rates determined in a flow-through experiment in a column initially packed with glass beads and aragonite grains and injected with a supersaturated solution. Image subtraction after 3-D registration of the data sets makes it possible to localize the newly formed crystals (Fig. 12a) that precipitated primarily around the aragonite grains (a few crystals also developed locally on glass beads). The sizes of the new crystals indicate that precipitation rates were heterogeneous over the duration of the experiment. Different rates are linked to chemical conditions at the fluid–rock interface at the micrometer or smaller scale (e.g., the local saturation index Ωloc), growth rate according to the different spar crystal faces, growth competition between crystals, and nucleation stage. Note that there is no indication that crystals have nucleated and grown on glass beads from the very beginning of the experiment. Combining a distance transform of the original material with the newly formed crystals along with the determination of their local maxima and propagation of these maxima allows for determination of the growth velocity of every new crystal, from Equation (6) (Fig. 12b). The growth velocity is shown to be between 0 (i.e., no precipitation) and 4.5 μm·d−1.
Rates of dissolution/precipitation reactions in fractures
Application of Equation (7) to dissolution within a planar fracture reduces to:
An example of a calcite dissolution rate calculation from a 95-hour flow-through experiment in an artificially fractured limestone is shown in Figure 13. The experiment consisted of injecting deionized water equilibrated with a partial pressure of CO2 of 1 bar at a constant flow rate Q of 10 cm3·h−1. Solution pH was determined to an average of 3.9 and 5.5 at the inlet and outlet, respectively. A trapezoid fracture of 15 mm long and 7.5 mm wide was obtained after stitching two flat (polished) fracture walls together. To induce a heterogeneous flow field over the fracture width, the fracture aperture was varied from 20 μm (shorter base) to 45 μm (longer base). As the rock used was a very fine-grained and pure limestone, the reactive surface area can be assumed to be almost constant over space and time. A map of fracture aperture at the end of the experiment shows the development of a principal flowpath within which dissolution has been enhanced, and this corresponds to the portion of the fracture with a larger initial aperture.
Assuming that the fracture can be represented as a continuum, dissolution is steady-state, and reactive surface area scales with the XMT geometric surface area (see further discussion), it is possible to recalculate the apparent kinetic rate (kr′) and even a theoretical pH based on kinetic law. Taking S = Sr, for example, and the kinetic law obtained by Plummer et al. (1978) (kr′ = k1aH+ + k2aCO2 + k3) yields:
Results presented in Figure 13 demonstrate that kr′ measured locally in the fracture void integrates transport effects, the highest values being obtained within the main flowpath, where residence time is shorter and the average flow velocity higher. Formation of preferential flow pathways is often associated with a reduction in the flux of dissolved species at the outlet (Noiriel et al. 2005; Luquot and Gouze 2009). This could be interpreted as the result of modification of the reactive surface, although it appears more likely that it is actually the result of a reduction in the “transport efficiency” in areas where the flow rate was reduced. The result is a situation in which transport progressively shifts from advective–dispersive-dominant to diffusion-dominant, as a result of the decrease in fluid velocity.
Comparison with the experimental results shows that the average pH is overestimated by the calculation at the inlet (pHcal ~ 4.1) and underestimated at the outlet (pHcal ~ 4.4). These differences can arise from the kinetic formulation used, the estimation of the surface area, the assumption of a continuum (i.e., constant pH values across the fracture walls) and steady-state dissolution. For the calculation, the dissolution rate was taken to be constant everywhere within the fracture over the course of the experiment, although flow-field calculations indicated that the flow velocity decreased by one order of magnitude. However, only a direct or inverse fully coupled modeling approach of reactive transport, taking into account feedback between chemistry, flow, and transport at the micro-scale, associated with upscaling from the pore scale to the sample-scale, could allow for a proper interpretation of the XMT observations. As shown in Figure 13, the initial heterogeneous flow field within fractured (or porous) media can quickly lead to localization of the flow through the primary flow paths (see also Fig. 16) as a result of heterogeneous transport along the different flow paths. However, incomplete transverse mixing across the fracture walls may also result in heterogeneous rates of reaction, particularly when there are differences between areas of higher and smaller fracture aperture, and for which velocity profiles vary. Li et al. (2008) showed that development of concentration gradients within single pores or fractures can lead to scale-dependent reaction rates. The highest discrepancies between pore-scale and continuum-scale models were observed for a mixed kinetic control (i.e., similar characteristic times for transport and surface-reaction) due to comparable rates of chemical reaction and advective transport. Molins et al. (2014) also compared pore-scale with continuum-scale simulations of calcite dissolution inside a capillary. They pointed out the effect of larger diffusive boundary layers formed around grains and in slow-flowing pore spaces that exchanged mass by diffusion with fast flow paths. Overall, the assumption of well-mixed conditions—the conceptual basis of the continuum scale model—did not apply in the capillary tube because of the local importance of transport limitations to the bulk (or effective) rates.
Porosity and permeability development in porous media and fractures
Porosity and permeability development in porous media
The relationship between porosity and permeability in rocks is complex and varies according to the porous and mineral networks and the geological processes involved (Bourbié and Zinszner 1985). Several models linking porosity to permeability have been proposed (e.g., Carman 1937; Bear 1972), but their predictive capacity is usually limited due to the lack of in situ observations of the parameters that control the dynamic of the porosity–permeability relationship. Combining XMT observation and determination of porosity, along with the measurement of permeability, allows us to examine more closely the relationship between these two major reservoir properties. Noiriel et al. (2004), for example, linked two distinct porosity–permeability power law relationships (Fig. 14a) to different processes that occurred successively during a flow-through dissolution experiment. First, a permeability increase of an order of magnitude for a porosity change of only 0.3% was linked to micro-crystalline particle migration within the rock (Fig. 7a). Only a small fraction of the pore space was involved at that stage of the experiment. Permeability increased later in the experiment as both the pore wall roughness decreased and the pore network connectivity increased, as shown in Figure 14.
Porosity–permeability evolution in reactive rocks results from the physical evolution of the pore space, which in turn results from interplay between reaction kinetics and advective and diffusive transport. Possible feedback between the flow regime and geochemical alteration can lead to instabilities and localization or divergence of the flow, depending on whether dissolution or precipitation is involved. The two parameters used for characterizing and predicting these phenomena are the Péclet (Pe) and Damköhler (Da) numbers (e.g., Steefel and Lasaga 1990; Hoefner and Fogler 1998), defined locally as: Pe = vL*/Dm and Da = krL*2/Dm, where v is the fluid velocity (m·s−1), Dm is molecular diffusion (m2·s−1), kr is a first order kinetic constant (s−1), and L* is a characteristic length (m), e.g., the fracture aperture or pore size.
Starting from a very homogenous rock in terms of structure and mineralogy, the changes in the rock geometry resulting from dissolution generally follow a pattern determined by the couple of the Damköhler and Péclet numbers. When reaction kinetics are slow compared with transport (low Da), dissolution is rather uniform and ramified reaction fronts result (Hoefner and Fogler 1998; Golfier et al. 2002), whereas compact and channelized dissolution can be observed at high Da, with a dissolution front advancing from the point of injection. Between these values, conical, dominant, or ramified wormholes can form. The patterns are directly linked to dissolution instabilities: as the reactive fluid infiltrates areas of higher permeability, a positive feedback between transport and chemical reactions develops, and leads to the growth of the wormholes (Ortoleva et al. 1987; Steefel and Lasaga 1990; Daccord et al. 1993).
Several experiments have been carried out to investigate the effect of dissolution regime on porosity development in the context of CO2 sequestration. Different values of Da were explored by either adjusting the fluid velocity (Luquot and Gouze 2009; Luhmann et al. 2014) or reactivity (Carroll et al. 2013; Smith et al. 2013a). Luquot and Gouze (2009) observed homogeneous dissolution over the length of a core (named D3) at the highest flow rate (lowest Da), whereas gradients in reaction developed at lower flow rates (higher Da), resulting in wormholes associated with highly heterogeneous porosity development within two other cores (named D1 and D2, Fig. 15). Smith et al. (2013a) investigated the effects of pore-space heterogeneity on the development of dissolution fronts in a vuggy limestone and marly dolostone. They observed that a homogeneous pore-space distribution (90% of pore sizes differing by only one order of magnitude) resulted in stable, uniformly advancing dissolution fronts associated with porosity increase and only minor changes in permeability. Conversely, heterogeneous pore space distributions (90% of pore sizes spanning at least 3.5 orders of magnitude) resulted in greater variability in local fluid velocity and mass transfer rates, leading to the formation of unstable dissolution associated with dramatic permeability increases of several orders of magnitude.
Although the Péclet and Damköhler numbers can provide some useful information about the evolution of a particular system, these non-dimensional parameters must be used with care. First, the existence of feedback that enhances permeability implies that the fluid velocity is not constant with time. Second, the dissolution kinetics of a reactive fluid through a rock is not constant and can vary over one or two orders of magnitude between far-from-equilibrium and close-to-equilibrium states. As a result, both Da and Pe numbers evolve with space and time. Finally, heterogeneities present initially in rocks can exert a first-degree influence on the type of dissolution front and resulting relationship between porosity and permeability (Smith et al. 2013a), thus competing with dissolution instabilities linked purely to the couple (Pe, Da).
Porosity and permeability development in fractures
The development of wormholes can also be observed in fractures, even in very flat artificial ones. The following example illustrates how dissolution kinetics can lead to different dissolution patterns in fracture, even though Da and Pe values were similar at the inlet of the samples at the beginning of each experiment (i.e., at t = 0 and z = 0). The effects of reaction kinetic law have been explored through two flow-through dissolution experiments. Two different inlet fluids at pH 3.9, either deionized water + HCl or deionized water equilibrated with PCO2 = 1 bar, were injected at a flow rate Q = 100 cm3·h−1 in artificial, almost flat, limestone fractures of initial aperture ~25 μm. While the fracture injected with HCl shows the formation of dominant wormholes (Fig. 16a), the fracture injected with CO2 does not exhibit any particularly localized wormholes (Fig. 16b), despite the flow beginning to localize due to slight heterogeneities in the initial aperture field. The difference arises from different kinetic paths taken by the two different inlet fluids during calcite dissolution. The dissolution rate of calcite is influenced by pH, PCO2 and surface reaction (Plummer et al. 1978). As the fluid reacts with calcite, pH increases and dissolution rates decrease. However, the decrease is higher (~2 orders of magnitude) and more abrupt for HCl than for CO2 (~1 order of magnitude) (Fig. 17), as a result of the buffering effect of the carbonated species. This example shows clearly that very different dissolution patterns can develop despite similar initial Da and Pe values, due to the kinetic reaction path followed by the reactive fluid.
Effects of texture and mineralogy on complex porosity–permeability relationships and transport
Even as 3-D mineralogical mapping of rocks remains a challenge, advances in XMT have illuminated the role of rock mineral composition and spatial distribution during geochemical processes. Primary rock texture and mineralogy have an important influence on the evolution of flow and transport properties, as pointed out in recent studies at the pore scale. For example, opposite trends of the porosity and permeability evolution have been observed, despite the fact that net dissolution occurred in the system (e.g., Noiriel et al. 2007a; Ellis et al. 2013; Luquot et al. 2014a). Ellis et al. (2013) even reported opposite outcomes with respect to fracture permeability evolution for two very similar experiments performed on nearly identical rock samples. In the experiment of Noiriel et al. (2007a), the heterogeneity of the dissolution rates (differing by about one order of magnitude for calcite and dolomite (Chou et al. 1989) and ten orders of magnitude for clays and carbonates (Köhler et al. 2003)) of the minerals comprising the rock matrix led to the development of dissolution heterogeneities at the fluid–rock interface. Increase of tortuosity and fracture roughness led in turn to a reduction in permeability, despite net dissolution. Roughness and tortuosity increases associated with large differences in the dissolution rate of minerals were also reported by Gouze et al. (2003), Noiriel et al. (2007a, 2013), and Ellis et al. (2011) (Fig. 18; see also Fig. 8).
Substantial reductions in permeability were also observed or implied in several experiments after the matrix-forming minerals had lost their cohesion and individual grains were being removed and transported by the fluid, thus contributing to pore clogging (Noiriel et al. 2007b; Ellis et al. 2013; Mangane et al. 2013; Qajar et al. 2013; Sell et al. 2013).
An increase in the roughness of the fluid–rock interface due to spatially and mineralogically heterogeneous reaction rates affects the transport of reactants and dissolved species to and from the fluid–mineral interface. For example, Noiriel et al. (2007b) observed an exponential decrease with time of the flux of dissolved species in a fractured argillaceous limestone during dissolution. By increasing diffusion compared to advection close to a mineral surface, dissolution of calcite grains combined with the development of a micro-porous clay coating nearby the fluid–mineral interface resulted in a decrease of transport “efficiency” (Fig. 19).
The scale dependence and the effects of mineral spatial distribution on reaction rates is also an important topic, as demonstrated in various experimental and numerical modeling studies (Li et al. 2006, 2007, 2013; Kim and Lindquist 2013; Salehikhoo et al. 2013). Such effects originate from the differences in the rates of mass transport between reactive and non-reactive pores and also depend on the spatial distributions of reactive minerals (Li et al. 2007). Geometrical changes accompanying dissolution of anorthite in a sandstone also resulted in continuously changing upscaled reaction rates, as shown by Kim and Lindquist (2013) using a pore-network model that took into account changes in pore volume, reactive surface area, and topological changes.
Effects of pore-scale heterogeneity on permeability reduction
As discussed above, permeability is closely linked to pore-scale heterogeneity. In the case of mineral precipitation, the location of the newly formed crystals is constrained by the induction time for nucleation (i.e., the time required before growth can initiate) and the degree of epitaxy between pre-existing mineral surfaces and newly formed crystals. In addition, the growth rate of newly formed crystals depends on their orientation (linked to the orientation of pre-existing mineral faces when growth is epitaxial), the different facets exhibiting different growth rates (e.g., Hilgers and Urai 2002; Nollet et al. 2006). Growth competition between neighboring crystals also involves variations of growth rate over time. Geometry changes resulting from growth at the fluid–mineral interface can have a significant effect on the hydrodynamics and transport fluid within a porous network. An example is shown here for different nucleation-and-growth scenarios derived from an experimental study of precipitation in columns packed with glass beads and calcite spar crystals (Noiriel et al. 2012). The growth of newly formed calcite crystals was epitaxial on calcite spar, leading to a smooth surface; conversely sparse rhombohedra of calcite developed on glass beads, leading to a rougher surface (Fig. 6). The experiment illustrates, that under similar geochemical and flow rate conditions, the variability of location, shape, and size of the new crystals depend on the pre-existing mineral substrate (i.e., either calcite spar or glass beads).
To examine the effects of evolving surface roughness on pore hydrodynamics, different growth scenarios were generated and combined with flow modeling. Starting from an initial sub-volume of 300 × 300 × 300 voxels, an algorithm was used to nucleate and grow cubic crystals at the calcite spar surface to mimic experimental observations (Fig. 20; see also Fig. 6) so as to compare the evolution of permeability between uniform (case 1) and heterogeneous (case 2) case of precipitation. Permeability was calculated by solving Stokes’ equations based on a volume averaging technique (Batchelor 1967; Quintard and Whitaker 1996). The pores with heterogeneous precipitation (case 2) exhibited a higher reduction in permeability in the simulations compared with the uniform case (case 1, Fig. 21a). Reduction in permeability is linked to the decrease in porosity that induces pore-throat size reduction. However, this example illustrates than a similar reduction in porosity can lead to different degrees of permeability reduction, as it also depends significantly on the fluid–mineral interface roughness, and so the location and shape of the newly formed crystals. Quantification of the fluid–solid surface roughness factor illustrates that trend (Fig. 21b).
The effect of roughness development on the transport of reacting species and precipitation rate near the interface is far more complex to assess, as the dynamics of growth are much more difficult to unravel in systems where flow velocity and species concentration fields evolve continuously along with the surface area of the crystal facets. Nollet et al. (2006) noted changes of growth rates of crystal faces through time, which they attributed to growth of individual facets, growth competition, flow direction, and changes of flow velocity. Local modifications of hydrodynamics and transport near the interface are likely involved in this case.
Rocks and other subsurface materials exhibit both a structural and behavioral complexity during hydro-mechano-chemical processes, which makes predicting flow and reactive transport in them challenging. Advances in XMT have begun to provide a better mechanistic understanding of the HMC coupling through direct rendering of the flow, transport, and chemical processes at the pore scale. XMT used together with physico-chemical measurements during experiments and finer-scale observations offers the opportunity to better access and understand reaction-induced changes of pore structure, permeability and surface reactivity, as a result of the interplay between geochemical reactions, changing the rock geometry, and the hydrodynamic and transport properties. In particular, the ability offered by 4D XMT to follow the movement of the fluid–rock interface and to measure local rates of dissolution/precipitation through time gives new insights about rock and mineral reactivity in open hydrological systems.
The potential of combining XMT experiments with numerical simulations to study pore-scale processes in rocks is far from mature, but shows encouraging developments. Upscaling experimental laboratory data to reservoir scale is one of the most challenging issues for predicting the long-term evolution of reservoirs in the presence of reactive fluids. Refining of macro-scale modeling will be possible once pore-scale modeling has reached maturity. Even if direct simulation and network modeling of flow and solute transport are now used routinely to calculate averaged properties at different scales, the coupling between chemical reaction, rock geometry modifications, and hydrodynamic and transport properties remains insufficiently constrained, particularly in physically and mineralogically heterogeneous rocks. Indeed, experimental observations have established some of the relationships between micro-scale features and macro-scale properties and this can lead in certain cases to potentially unexpected (emergent) effects. Our knowledge of fluid–rock interaction in chemically-reactive environments will have improved significantly once it is possible to integrate 4D XMT data with fully coupled modeling of reactive transport that include multicomponent chemistry, fluid velocity distribution and transport of solute species in fluid phase.
I am grateful for the constructive reviews of the manuscript and helpful comments provided by Maša Prodanovič, Marco Voltolini and Carl Steefel.