The pore space in rocks, sediments, and soils can change significantly as a result of weathering (see Navarre-Sitchler et al. 2015, this volume), diagenetic, metamorphic, tectonic, and even anthropogenic processes. As sediments undergo compaction during burial, grains are rearranged leading to an overall reduction in porosity and pore size (Athy 1930; Hedberg 1936; Neuzil 1994; Dewhurst et al. 1999; Anovitz et al. 2013). In addition, geochemical reactions can induce the precipitation and dissolution of minerals, which can either enhance or reduce pore space (e.g., Navarre-Sitchler et al. 2009; Emmanuel et al. 2010; Stack et al. 2014; Anovitz et al. 2015). During metamorphism too, mineral assemblages can change, altering rock fabrics and porosity (Manning and Bird 1995; Manning and Ingebritsen 1999; Neuhoff et al. 1999; Anovitz et al. 2009; Wang et al. 2013). As the pore space in geological media strongly affects permeability, evolving textures can influence the migration of water, contaminants, gases, and hydrocarbons in the subsurface.
Although models—including the Kozeny–Carman equation (Kozeny 1927; Bear 1988)— exist to predict the relationship between porosity and permeability, they are often severely limited, in part because little is known about how pore size, pore geometry, and pore networks evolve in response to chemical and physical processes (Lukasiewicz and Reed 1988; Costa 2006; Xu and Yu 2008). In the case of geochemical reactions, calculating the change in total porosity due to the precipitation of a given mass of mineral is straightforward. However, predicting the way in which the precipitated minerals are distributed throughout the pores remains a non-trivial challenge (Fig. 1; Emmanuel and Ague 2009; Emmanuel et al. 2010, Hedges and Whitlam. 2013; Wang et al. 2013; Stack et al. 2014; Anovitz et al. 2015). Similarly, the impact of mechanical processes on the development of pore-size distributions and the anisotropy of pore networks has yet to be resolved to a satisfactory level (Day-Stirrat et al. 2008a,b, 2010). An important step towards understanding how pores evolve in different geological scenarios is to develop mathematical frameworks that can describe the effects of geochemical reactions and mechanical compaction on porous rock matrices. Here, we present three different phenomenological models that describe the effect of geochemical reactions and mechanical processes on the evolution of pore-size distributions in geological media. As both of these processes can alter pore size, we also develop a framework for coupling the geochemical and mechanical modes of evolution. In addition, we review some of the attempts to use these models to simulate actual field and laboratory data, and explore some of the open questions and directions for future research.
PORE-SIZE EVOLUTION DURING MINERAL PRECIPITATION AND DISSOLUTION
Mineral precipitation typically reduces the size of voids while dissolution increases them, and the use of powerful computational methods has enabled such processes to be simulated directly at the pore scale. Lattice-Boltzmann models have been used to map the evolution of the solid–fluid interface (e.g., Kang et al. 2002, 2005, 2007; Huber et al. 2014), while pore network models—which replace the complex geometry of the pore space with a network of interconnected pores and channels—have provided insight into the scale dependence of geochemical rate laws (Li et al. 2006; Meakin and Tartakovsky 2009). However, although pore-scale models can be powerful tools, they are unable to resolve processes that operate at scales smaller than the grid size. In many rocks, pore sizes are highly multiscalar, ranging from the nanometer scale up to hundreds of microns or larger; modeling processes that simultaneously occur at these different scales represents a major challenge. In addition, there is significant evidence (Rother et al. 2007; Anovitz et al. 2013b; Hedges and Whitelam 2013; Kolesnikov et al. 2014) that as pores become smaller the physical properties of the fluids contained within them are strongly altered, further complicating matters. In any given representative volume there may be many millions of individual voids of varying size and shape. Thus, if geochemical processes are size-dependent, it may be more practical to adopt an approach that attempts to simulate the evolution of pores across a wide range of spatial scales.
One approach that couples a statistical description of pore size with a mechanistic view of their evolution involves using a partial differential form of a continuity equation to track temporal changes to pores and pore-size distributions (Or et al. 2000; Leij et al. 2002). In its simplest geochemical formulation, the equation effectively tracks changes to the pore size frequency as a result of precipitation or dissolution:
where fr is a function describing the number of pores, characterized by a radius r, per unit volume of rock, sediment, or soil; the variable can be thought of as akin to a probability density function for pore size. In this framework the parameter v is the rate of change of pore size in units of distance per unit time and represents the rate of contraction or growth of the pores. This rate parameter is dependent on the mineral being precipitated or dissolved, as well as environmental conditions such as solution chemistry, temperature, and pressure, and the parameter can be dependent on pore size. This framework represents a much simplified version of other models describing pore size evolution (Or et al. 2000; Leij et al. 2002) in that it assumes that the total number of pores is conserved and that pores can only change as a result of geochemical processes. As a result, additional diffusional and sink terms are not included.
In the simplest scenario, v is a constant, and it can be related to a standard expression for the geochemical reaction rate by (Lasaga 1998)
where vm is the molar volume of the precipitating or dissolving mineral, k is a rate constant, Ω is the bulk degree of saturation (i.e., the ratio of the ion activity product to the solubility product), and a and b are positive exponents. In this convention if the sign of v is positive, dissolution occurs and pore size increases. Importantly, a constant value for v implies that the pores of all sizes respond to geochemical conditions in the same way so that all pores contract or expand at the same rate. Moreover, such behavior implies that precipitation will cause small pores to close first, although during diagenesis some studies have reported that this is not necessarily the case (Putnis and Mauthe 2001; Emmanuel et al. 2010; Anovitz et al. 2013a, 2015) Deviation from a uniform pore size reduction could be due to preferred flow paths and higher reactant fluxes though large pores. Alternatively, the effects of interfacial energy in very small pores could make v size-dependent (Emmanuel and Berkowitz 2007; Emmanuel and Ague 2009; Emmanuel et al. 2010; Hedges and Whitelam 2013; Rother et al. 2007; Anovitz et al. 2013b; Kolesnikov et al. 2014). Such effects will depend both on pore size and pore geometry. For pores with convex walls (see Fig. 2) as opposed to pores with concave walls (e.g., spherical pores) the expression could have the form (Emmanuel and Ague 2009):
where γ is the interfacial energy, θ is the dihedral angle between adjacent pore walls, R is the gas constant, and T is temperature. Crucially, this size-dependent modification effectively reflects the increased solubility of tiny crystals with high curvature. This means that grains with low curvature (large grains) can grow in big pores at the same time that smaller, high-curvature grains are prevented from growing in small pores. However, this effect is only expected to be important in nano-scale pores, or at very low levels of supersaturation. Elevated temperatures, too, may reduce size-dependent effects. For a much more detailed discussion of the solid–fluid interface at the nanometer-scale and its influence on fluid structure and fluid–solid interactions the reader is referred to other papers (e.g., Zhang et al. 2004, 2006; Vlcek et al. 2007; Mamontov et al. 2009; Wesolowski et al. 2012).
While Equation (1) possesses the relatively straightforward structure of an advection equation, the variable fr (i.e., the density of pores of a given size) is not a commonly reported parameter, in part due to the difficulty of assigning a representative size to the irregularly shaped pores typically found in natural porous media (Anovitz et al. 2009, 2013a, 2015). Descriptions of the various techniques used to measure pore-size distributions are given by Anovitz and Cole (2015, this volume) and Navarre-Sitchler et al, (2015, this volume). Such methods often report size distributions in terms of the volume of pores of a given size (or, more correctly, in a certain interval of sizes) per volume of sample. Using such data a cumulative porosity curve can be obtained so that the total porosity ϕ in the interval r1 to r2 is given by the following integral:
where ϕr is the probability density function for pore size multiplied by the total porosity. Practically speaking, ϕr can be calculated from pore size data by numerically differentiating the cumulative porosity curve with respect to pore size (e.g. Anovitz et al. 2015). If the pores of a given size possess a characteristic volume, Vr, then fr = ϕr /Vr, and Equation (1) becomes
For the simplest case in which the pores are spherical it follows that
and using the quotient rule it can be shown that the expression simplifies to
The extra term on the right hand side of the equation is a geometric term that effectively serves to increase overall porosity when pores expand, and reduce overall porosity when pores shrink. Similar expressions can be developed for different pore geometries; for example, in the case of cylindrical pores undergoing radial shrinkage it can be shown that the geometric term is 2 vϕr /r, while for cuboid pores with uniaxial shrinkage the term is vϕr /r. These equations can be readily solved for a range of initial distributions by numerically solving 1D partial differential equations.
In its basic form the model assumes a constant supply of reactants, although in principle it could be coupled to a reactive transport equation that conserves mass. While a fully coupled model would be more realistic, the main way in which pore-size distributions evolve can be shown by considering Equation (7) on its own. When v is constant, the model, as expected, predicts that as pores are filled, smaller pores are the first to be filled. To demonstrate the way pores of different sizes are predicted to change, Figure 3 shows a simulation of a bimodal distribution: the peak initially located at 30 μm is much higher than the peak at 60 μm; however, by the end of the simulation both peaks are approximately equal in height. While the v parameter determines the rapidity with which the curves evolve, the model predicts that under constant conditions pore-size distributions will evolve along fixed paths.
Although such models contain a number of simplifying assumptions, they can nevertheless provide important insights into the processes governing the evolution of pore-size distributions during mineralization. In a study examining varying levels of quartz cementation in sandstone from the Stø formation in the Barents Sea (Fig. 4; Emmanuel et al. 2010), reduction in pore size was found to be non-uniform across different sizes, so that pores with radii smaller than approximately 10 μm appeared to undergo very little change (Fig. 5a). Fluid flow in the system was thought to be minimal, and the study used a complex diffusion-reaction equation with moving boundary conditions to simulate the evolution of pore sizes. However, Equation (7) on its own can be used to model the evolution of the porous matrix from high porosity to low porosity in a much more straightforward way. Using the pore-size distribution in the least-cemented sample as an initial condition, the constant-v model performs adequately at simulating the shift of the peak initially located at 13.5 μm to ~12.3 μm (Fig. 5b). However, the model also predicts that the peak originally appearing at 8.5 μm should shift to around 7.4 μm and that the peak at 2 μm should disappear altogether. In contrast, in the low-porosity sample, both measured peaks remain unchanged. Significantly, the preferential filling of larger pores has been reported in other studies of sandstones (Putnis and Mauthe 2001; Anovitz et al. 2013a, 2015) and this type of behavior may be relatively common.
There are a number of mechanisms that could lead to the preferential filling of larger pores. Clay coatings in small pores may inhibit quartz cementation, and the timing of clay mineral formation, as well as the extent of coverage, may be critical factors. Alternatively, limited connectivity and subsequent limited solute flux in small pores may restrict the supply of dissolved silica available for quartz precipitation (Emmanuel et al. 2010). However, both simulations and measurements indicate that interfacial energy effects could also play an important role (Emmanuel and Berkowitz 2007; Emmanuel and Ague 2009; Emmanuel et al. 2010): when using the size-dependent expression for v in Equation (3), the model produces an excellent match with the measured data (Fig. 5c). Critically, in this simulation only two variables were changed to obtain a good fit: k and Ω. All other parameters were taken from the literature and are given in the caption for Figure 5. Importantly, this model suggests that the level of supersaturation during cementation was extremely low with Ω = 1.0004, a value that is similar to that calculated by Emmanuel et al. (2010) using a more complex reactive transport model. At higher levels of supersaturation, pores smaller than 10 μm are also filled, so that the model behavior evolves towards that of the constant-v model.
One of the main limitations of the pore-size distribution approach presented here is that although the sizes of the pores can be predicted, their spatial ordering and connectivity are not expressed. As pores get smaller, flow through throats can be restricted, and when throats < 10 nm the overlap of electrical double layers (EDL), particularly in clay-rich media, can further modify transport rates (see Tournassat and Steefel 2015, this volume). As a result, such models cannot, on their own, be used to simulate the development of preferential pathways for flow, and they are, therefore, more appropriate for regimes in which fluid flow is minimal. Despite this limitation, the simulations show that the integration of geochemical reaction kinetics into a continuity equation for pore-size evolution can produce realistic results, shedding light on the existence of size-dependent processes. While the model will not be able to describe all scenarios involving mineral reactions, it should be applicable to systems in which the rate of mineral reaction is limited by surface kinetics (i.e., the reaction rate is slow compared to the flux of reactant transport) as the concentrations of reactants in different pores are expected to be uniform. Such conditions are likely to prevail for the slow rates of quartz precipitation in many diagenetic scenarios.
PORE-SIZE EVOLUTION DURING MECHANICAL COMPACTION
While geochemical reactions can strongly influence the porosity and pore-size distributions in soils and rocks, mechanical compaction can also be an important factor in many scenarios. Sediment burial, soil tillage, and aquifer pumping can all lead to pore shrinkage and even pore collapse. To demonstrate the difficulties associated with developing models to describe pore-size evolution during compaction, we start by adopting a similar approach to that used in the previous section and consider the rate of change in the density of pores of a given size, r, as a result of an applied stress, σ. From the chain rule this can be expressed as:
where nr is a function describing the number of pores characterized by length scale r per unit mass of soil or rock. Note that this is slightly different from the fr function in Equation (1) which is defined as the number of pores per unit volume. This is because during compaction it is easier to track changes to a given mass than changes to a given volume. A further difference between Equation (1) and Equation (8) is that, for the case of mechanical compaction, the rate of change is expressed with respect to σ rather than time. As a result, this model is not temporally dynamic and considers only stress-induced changes at equilibrium. Such an approach is likely to be appropriate for systems in which the rate of change in stress is very slow, which is typically the case during burial and diagenesis. For simplicity, we adopt here a simplified approach which assumes that: (i) the pores are spherical; (ii) uniform compression is applied; and (iii) the material is homogenous with an effective bulk modulus, B. Using these assumptions, it can be shown that Equation (8) becomes
As the pores are spherical they have a characteristic volume, Vr, so that nr= ɛr/Vr, where ɛr is the pore space per unit mass associated with each pore size. Substitution of the definition for nr yields
which, for the case of spherical pores, simplifies to
In this formulation, the effective bulk modulus, B, is the inverse of the sediment compressibility. However, it should not be regarded as a true elastic modulus as it reflects both the elastic compression of the grains and their inelastic re-organization. As a result, B is not expected to be constant during the compaction process, although a stress-dependent B will affect only the scaling of the solution with respect to σ and not the evolution of the pore-size distributions at each stage of compaction.
Importantly, even though the rate of change of pore size is expressed in terms of stress, it can be seen that the structure of Equation (10) is similar to that of Equation (7) developed for the case of pore-size evolution due to mineral reactions. When the compaction is expressed in terms of changes due to stress, and the time scales of the processes are not of interest, this equation can be used to examine the compaction process. Although a slightly different equation was employed, such an approach was used to analyze the evolution of pore-size distributions during the compaction of Boston Blue Clay with applied stresses that ranged from 0.1 MPa to 10 MPa (Day-Stirrat et al. 2011; Emmanuel and Day-Stirrat 2012). The experimental data, obtained through mercury injection porosimetry, indicated that the main pore size peak at 110 nm shifted to approximately 65 nm (Fig. 6a); however this dramatic change was accompanied by a relatively small drop in porosity from 41.4% to 37.2%. Also puzzling, is that the 65 nm peak in the compacted sample is actually more pronounced than the primary peak in the unconsolidated sample, suggesting that shrinking pores in effect get stuck at that scale. Given these observations, it is perhaps unsurprising that the model (Eqn. 11) produces an extremely poor fit with the data (Fig. 6b), predicting that the shift in pore size should be accompanied by a much greater reduction in porosity. Emmanuel and Day-Stirrat (2012) suggested that the behavior of the system could be due to pores that become increasingly more resistant to compaction as they become smaller than 100 nm in radius, an effect that is most likely due to the difficulty of rearranging locked grains that have become tightly squeezed together. In addition, electrostatic repulsion between grains might also become significant at the sub-100-nm scale, and this could contribute to the resistance to grain compaction. The simplest way to adapt the model to account for such processes is to make B pore-size dependent, which would reflect different levels of effective compressibility in the material surrounding pores of different sizes. A similar approach was suggested by Emmanuel and Day-Stirrat (2012), and an increasingly stiff matrix at smaller scales was shown to be able to improve the fit between the measurements and simulations. However, further study is required to determine the most appropriate functional form required to describe the potential size dependence of the effective bulk modulus.
Although compaction is often considered to be dependent only on the applied stress, in many applications the temporal evolution of pore-size distributions may also be of interest. A time dependent equation can be obtained by multiplying Equation (10) by the temporal derivative of σ (i.e., the rate of applied stress), and defining the parameter α as α = ∂σ/∂t yields
Moreover, as ɛr= ϕr /ρb where ρb is the bulk density, the equation can be expressed in terms of ϕr such that
and it can be shown that
We note that both during compaction and geochemical reactions ρb will evolve over time; however, ρb is related to ϕr by
where ρs is the density of the solid phase. Thus, Equation (13) is actually an integro-differential equation with the form
While Equation (15) can be solved numerically, at present the main obstacle to testing the model’s ability to describe real systems is the absence of reliable datasets documenting the evolution over time of pore-size distributions during compaction.
CHEMO-MECHANICAL COUPLING AND PORE-SIZE EVOLUTION
In many situations geochemical and mechanical compaction do not act independently of each other, but rather act in concert, and the two approaches developed in the previous sections can be integrated to yield a single expression for the evolution of pore-size distributions. Summing the right hand sides of Equations (7) and (15) yields a single expression:
While solving this equation is a non-trivial task, it represents an important framework with which to tackle evolving pore-size distributions in geological scenarios in which chemo-mechanical coupling exists.
One scenario in which both mechanical compaction and mineral reactions can influence the evolution of porosity is during the intrusion of magma into sediments and concurrent contact metamorphism. In the country rock surrounding the intrusion, porosity is often filled as a result of reaction, recrystallization, and grain compaction (Baer 1991; Summer and Ayalon 1995; Young 2008). For example, scanning electron microscopy images of sandstone from the Jurassic Inmar formation, sampled near an intrusion in the Ramon Crater in Israel, reveal that metamorphosed sandstone has a much lower porosity than unaltered sandstone (Fig. 7a–b). In addition, backscattered electron imaging reveals the presence of minerals, such as kaolinite and Fe-oxides, in the quartzite that are not present in unaltered sandstone (Summer and Ayalon 1995). While the images clearly show a significant level of porosity reduction due to recrystallization and grain rearrangement, additional processes are also revealed from an analysis the pore-size distribution data obtained using the (U)SANS technique detailed by Anovitz and Cole (2015, this volume). Peaks at pore sizes of 1 nm and 10 nm present in the unmetamorphosed country rock are not present in the metamorphic quartzite (Fig. 7c–d), which could indicate the healing of micro-cracks or infilling of tiny pores during the metamorphic process (Brantley et al. 1990; Brantley 1992). Interestingly, this result is the opposite of that found during the diagenetically controlled mineralization in sandstones discussed earlier, which can leave small pores unfilled (Emmanuel et al. 2010), but similar to the changes was seen in metamorphosed marls (Wang et al. 2013). Such differences could be due to a higher level of supersaturation during metamorphic processes, preferential nucleation in small pores (Hedges et al. 2013), or mechanical compaction. In future studies, models such as that detailed in Equation (17) should help constrain the mechanisms and conditions influencing the evolution of pore-size distributions in complex diagenetic and metamorphic scenarios.
In this paper, we have outlined the equations governing the evolution of pore-size distributions during diagenesis and metamorphism. Three scenarios for changing porosity were explored: (i) mineral precipitation and dissolution; (ii) mechanical compaction; (iii) coupled chemo-mechanical effects. We showed that when compared with experimental and field data, such models can be used to identify important processes governing the evolution of voids at different spatial scales.
Although the equations detailed here provide a theoretical framework for the analysis of porosity evolution, a number of important challenges can be identified. Due to the lack of a diffusive term, the partial differential equations are relatively stiff, so that the numerical solutions can be unstable, particularly when the simulated processes are size dependent. Moreover, the equations in which the bulk density changes need to be accounted for (i.e., Eqns. 15 and 16) have not yet been solved. Software capable of solving integro-differential equations could provide a route to mapping the behavior of these more complicated models.
In addition to the practical difficulties associated with solving the partial differential equations, there are a number of more fundamental limitations associated with the approach explored here. One basic assumption is that all the pores undergo a gradual, predictable evolution as a result of mineral infilling or mechanical squeezing, and such an approach does not allow for stochastic processes that may lead to pores of any given size evolving in different ways. Moreover, during cementation existing voids may be filled by porous cement or newly formed minerals (Fig. 8), effectively subdividing large pores into smaller ones, so that relatively high levels of porosity are maintained even though a drastic reduction in pore size may occur. Such processes may be accounted for by introducing additional diffusive or sink/source terms in the equations, although this will require further parametrization of the model. On an applied level, even if models can reliably predict the evolution of pore-size distributions, the way such data can best be used to determine changes to permeability has yet to be resolved.
In many systems, in addition to the mechano-chemical processes, fluid flow can also add an extra layer of complexity. In addition to transporting dissolved species through the porous matrix, cement holding grains together may dissolve, allowing solid grains to be mobilized and transported as colloids through the porous matrix. Grain detachment has been observed in some studies of carbonate rocks (Emmanuel and Levenson 2014; Luquot et al. 2014), and the process may enhance porosity in regions in which grains are removed, while simultaneously reducing porosity in regions where grains are deposited. At present, determining how this mechanism influences spatial patterns of porosity and permeability requires new experimental data and more sophisticated reactive transport models.
One additional aspect crucial to the coupling between geochemical and mechanical process not addressed by the models presented here is the impact that porosity evolution has on the mechanical properties, such as the bulk modulus or fracture strength, of the geological media. Statistical techniques aimed at obtaining these data from descriptions of pore structures are summarized elsewhere in this volume (Anovitz and Cole 2015, this volume). Cementation typically increases rock strength and stiffness, and this is expected to have an effect on the evolution of pore size and permeability. There is an abundance of empirical data concerning the mechanical properties of rocks with different levels of porosity, but the effect that microstructural evolution has on rock strength is still not well understood. Thus, both experimental and theoretical studies aimed at developing a comprehensive framework linking geochemical reactions to the strength of porous geological materials represents an important direction for future research.
SE thanks the Israel Ministry of National Infrastructures, Water, and Energy. Effort by LMA was supported by research sponsored by the Division of Chemical Sciences, Geosciences, and Biosciences, Office of Basic Energy Sciences, U.S. Department of Energy. We acknowledge the support of the National Institute of Standards and Technology, Center for Neutron Research, U.S. Department of Commerce, and the High-Flux Isotope Reactor at the Oak Ridge National Laboratory in providing the research neutron facilities used in this work. This work utilized facilities supported in part by the National Science Foundation under agreement No. DMR-0944772. Certain commercial equipment, instruments, materials and software may be identified in this paper to foster understanding. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology, the Department of Energy, or the Oak Ridge National Laboratory, nor does it imply that the materials or equipment identified are necessarily the best available for the purpose. We are also grateful to Ben Gilbert and Carl Steefel for their helpful reviews.