## Abstract

In sedimentary basins and volcanic edifices, granular materials undergo densification that results in a decrease of porosity and permeability. Understanding the link between porosity and permeability is central to predicting fluid migration in the sedimentary crust and during volcanic outgassing. Sedimentary diagenesis and volcanic welding both involve the transition of an initially granular material to a non-granular (porous to dense) rock. Scaling laws for the prediction of fluid permeability during such granular densification remain contested. Here, based on collated literature data for a range of sedimentary and volcanic rocks for which the initial material state was granular, we test theoretical scaling laws. We provide a statistical tool for predicting the evolution of the internal surface area of a system of particles during isotropic diagenesis and welding, which in turn facilitates the universal scaling of the fluid permeability of these rocks. We find agreement across a large range of measured natural permeabilities. We propose that this result will prove useful for geologists involved in modeling porosity-permeability evolution in similar settings.

## INTRODUCTION

Initially granular materials deposited in sedimentary basins commonly densify by compaction and cementation (Fowler and Yang, 1998), and in volcanic systems by viscous processes of compaction and sintering (Quane and Russell, 2005; Wadsworth et al., 2014). These processes in turn significantly change the microstructure and the fluid permeability through the rocks (e.g., Bourbie and Zinszner, 1985; Heap et al., 2015). To date, no effective scaling of the permeability sufficient to describe comprehensively the entire process—from initially granular to finally dense and coherent—has been developed. This limits our ability to model the fluid flow in these systems pertinent to hydrodynamic fracturing efforts in hydrocarbon exploration (Norris et al., 2015), petroleum reservoir migration and storage (Honarpour et al., 1986), volcanic degassing time scales (Edmonds and Herd, 2007), and the propensity of magma to shear fracture (Mueller et al., 2008), amongst many other geological applications. While the trajectory of the porosity-permeability relationship during magma vesiculation is relatively well explored for many of Earth’s volcanic systems (Blower, 2001; Klug and Cashman, 1996; Mueller et al., 2005; Saar and Manga, 1999), the evolution of permeability during densification of an initially granular material is less well constrained (Blair et al., 1993; Bourbie and Zinszner, 1985; Martys et al., 1994). A deficit in the scaling efforts addressed herein has been that microstructural controls are rarely explored in detail when global approaches are attempted (e.g., Mueller et al., 2005).

## THE GRANULAR TO NON-GRANULAR TRANSITION IN MAGMAS AND SEDIMENTS

Magmas readily fragment to form granular materials (volcanic ash) upon rapid decompression (Alidibirov and Dingwell, 1996), or more generally when local strain time scales approach the structural relaxation time scale of the liquid (Dingwell and Webb, 1989), leading to the generation of in situ and proximal volcaniclastic deposits. Fragments are also deposited in clastic sedimentary processes. Such fragments, sedimentary or volcanic, are variably packed in granular beds which experience changes in fluid flow, pressure, and temperature, and undergo diagenesis or welding, densifying at the expense of pore space (Bourbie and Zinszner, 1985; Fowler and Yang, 1998). This densification also commonly occurs within volcanic interiors where fragmented magma is deposited in cracks within the still-ductile magma itself (Tuffen and Dingwell, 2005) or in subaerial volcanic basins (Branney and Kokelaar, 1992) where fragmental volcaniclastic deposits may remain hot enough to relax viscously (Dingwell and Webb, 1989) and to undergo densification by sintering (Wadsworth et al., 2014). Both of these processes—sediment compaction and viscous sintering—involve the time-dependent evolution of porosity, ϕ, and permeability, *k*, that accompanies the material transition from a packed granular state through a denser, coherent porous state to an impermeable, fully densified state (Fig. 1). This impermeable state is typically described in terms of the attainment of a percolation threshold, ϕ_{c}, at which the pores become isolated. For systems undergoing bubble nucleation, growth, and coalescence in a liquid, percolation thresholds of 0.2 ≲ ϕ_{c} ≲ 0.7 have been documented (Blower, 2001; Klug and Cashman, 1996; Mueller et al., 2005; Okumura et al., 2008; Westrich and Eichelberger, 1994). In contrast, during isotropic densification of particles, the percolation threshold can be much lower with ϕ_{c} ≈ 0.03 (Rintoul, 2000). This wide range in threshold values implies that densifying, initially granular systems are capable of maintaining permeable pathways to very low porosities.

Understanding how permeability evolves during the granular densification process will help us to better understand fluid pressurization that can lead to fracturing. In sedimentary basins, unexpectedly high pore fluid pressures can cause operational failure during exploration drilling, by drill-hole collapse (Fowler and Yang, 1998). Similarly, in volcanoes, it is the time scales and pressures associated with outgassing-welding processes that are thought to most readily influence eruption likelihood (Quane and Russell, 2005). In silicic crystal-poor volcanoes, the majority of shallow outgassing is thought to occur through cracks filled with densifying hot volcanic ash. This may be the primary degassing mechanism in obsidian-forming magmas (Castro et al., 2014). The efficacy of fluid pressure equilibration through permeable outgassing is a critical parameter in the assessment of whether or not explosive eruption or borehole collapse will occur. The initially granular nature of these in-conduit deposits renders the scaling laws proposed for vesiculation processes inappropriate for the processes of interest to us in this study.

While measurements of fluid permeability in granular sediments undergoing densification are reasonably common (Blair et al., 1993; Bourbie and Zinszner, 1985), those made on initially granular volcanic rocks have only recently become the subject of research (Heap et al., 2015; Okumura and Sasaki, 2014; Wright and Cashman, 2014). In both cases, data from the experimental densification of granular material and quantification of resultant permeability changes are somewhat sparse. Here, we collate a range of permeability data that involve this granular to non-granular transition in volcanic and sediment-derived rocks (Fig. 2) and find a universal scaling. This scaling, when coupled with existing kinetic laws for the evolution of porosity (Fowler and Yang, 1998; Quane and Russell, 2005; Wadsworth et al., 2014), yields the promise of providing predictive tools for the time scales of permeability loss in surficial environments.

## UNIVERSAL SCALING FOR FLUID PERMEABILITY

Scaling relationships used to define fluid permeability in vesiculating systems commonly require knowledge of the fluid pathway tortuosity or pore-throat length scale (Blower, 2001; Klug and Cashman, 1996; Le Pennec et al., 2001), whereas the length scale best suited to scaling the permeability in granular systems is the particle size. However, in densification scenarios where the system qualitatively evolves from granular toward fully dense (Fig. 1), a universal prediction of *k* is hampered by a lack of continuous description of the microstructural length scale involved. For models of compaction in sedimentary basins, the simplest power-law form is most commonly used: *k* = *k*_{r} (ϕ / ϕ_{r})^{b}, where ϕ is the porosity, *r* denotes reference (e.g., the initial) values, and *b* can range from 2 (Connolly et al., 2009; von Bargen and Waff, 1986) to ∼8 (Fowler and Yang, 1998). However, the reference *k*_{r}(ϕ_{r}) depends on parameterization of the particle sizes and initial packing, which are not functionally constrained and thus not universal. The most commonly used empirical scaling for the permeability when fluid flow is around spherical particles (granular system) is the Kozeny-Carman relation *k* = ϕ^{3}/2*s*^{2} (e.g., Scheidegger, 1963), where *s* is the specific surface—the ratio of pore surface area to sample volume—such that the particle length scale can be used to predict *s*. Adaptations of this scaling are frequently used in models of melt extraction from the asthenospheric mantle (Connolly et al., 2009; McKenzie, 1984). These approaches have been widely applied. Yet they have failed to prove universal, even for idealized geometries (Martys et al., 1994). An easily implemented universal model is the one proposed by Martys et al. (1994), who suggested that *k* is related to *s* via:

_{*}= 1 – (ϕ – ϕ

_{c}) and

*b*is 4.2 and may prove to be related to the initial particle geometry. We can calibrate the utility of Equation 1 using samples of sandstone and of sintered glass beads (Blair et al., 1993) for which all of the relevant parameters are measured quantities. This produces reasonable agreement over a large range of both

*k*and ϕ when ϕ

_{c}≈ 0.03 (Fig. 3, inset). The value of ϕ

_{c}≈ 0.03 is in excellent agreement with universal predictions for ideal spheres undergoing densification (Rintoul, 2000). For the majority of both sediments and variably welded volcanic rocks, the value of

*s*is unmeasured. Below, we tackle this deficit by estimating how

*s*evolves with ϕ in densifying granular systems.

## PREDICTING THE SPECIFIC SURFACE DURING DENSIFICATION

In order to use Equation 1 effectively, we need to know how *s* varies from an initial state as ϕ decreases during volcanic welding or sediment compaction. To do this, we use a statistical approach for heterogeneous granular materials developed by Torquato (2013) who provided solutions for (1) the size of pores between packings of spherical hard particles, and (2) the specific surface of packings of spheres that can overlap with one another randomly, termed overlapping, or that are “hard” and do not overlap, termed non-overlapping. For a given initial value of ϕ and an initial characteristic monodisperse particle radius *R* (a polydisperse solution is provided in the GSA Data Repository^{1}), we can compute the characteristic or “effective” pore radius *a* that is nestled between the particles, assuming all particles are initially non-overlapping, by starting from the cumulative probability density function *F*(*x*) of the pore size where *x* = *a*/*R* and β = 1 + *x*:

Here *y*_{0}, *y*_{1}, *y*_{2} and *y*_{3} are parameters that depend on ϕ (given in the Data Repository). The *n*^{th} moment of the probability density function of *x* (termed 〈*x*^{n}〉) is then related to the cumulative probability density function *F*(*x*) by integrating . Thus, the mean (i.e., *n* = 1) pore radius, which we take as the characteristic size, is given by *a* = *R*〈*x*〉.

Finally, we can use this value of *a* to compute the specific surface that would exist during densification assuming all of the pores interstitial to the densifying material remain overlapping:

The overlapping assumption for the pore space is valid as pores readily interact and are not hard spheres, unlike the initial particles. The exception to the use of Equation 3 is for material that is granular and below the maximum packing porosity of the particles ϕ_{m}, such as is the case for the high-ϕ experimental portion of our data sets in which particles have only undergone incipient sintering (Table 1). For these data, we compute the surface area directly from the monodisperse *R* by the simpler geometric relationship *s*(*R*) = 3(1 – ϕ)/*R*. Using this for the data from experiments using crushed volcanic ash from Volcán de Colima, Mexico (J. Kendrick, 2015, personal commun.), yields an estimation of ϕ_{m} ≈ 0.32 for these angular, polydisperse particles.

Using Equations 2 and 3 to compute *s* for use with Equation 1 provides us with tools to assess the scaling of the permeability for natural samples that were initially granular and for which we have constraint of the initial particle size *R*. The initial values of ϕ are taken as the highest measured value, and, in order to minimize the assumptions made, only data for which the initial characteristic *R* is quoted are used. The measured values of *R* and computed *a* and ϕ_{c} are quoted in Table 1. Figure 3 shows the results of this scaling coupled with the results for the calibration scaling (Fig. 3, inset). We find good agreement (*r*^{2} = 0.96) across many orders of magnitude of normalized permeability. In the Data Repository, we provide this scaling for other combinations of overlapping and non-overlapping particles and pores to show that these work less effectively than the method presented here (0.56 < *r*^{2} < 0.93).

## DISCUSSION

The scaling presented here can be used to predict a granular system permeability when ϕ and *s* are measured (Equation 1). Alternatively, knowledge of *R* of the initial granular system permits the use of these tools to estimate *s* (Equations 2 and 3) and thus *k* easily. The conceptual framework developed here is novel because the definition of the length scale used to estimate *s* changes from *R* to *a* as the system crosses below ϕ = ϕ_{m} upon deposition and densification, which is when the particles are sufficiently closely packed to be considered densifying in a bed (ϕ < ϕ_{m}) and not depositing from dispersion (ϕ > ϕ_{m}).

Despite the progress made using this method, caveats to our approach nonetheless remain and these are important to state explicitly. First and foremost, our proposed scaling does not address anisotropy of the permeability, which can develop during compaction of both sediments and volcanic ash (Wright and Cashman, 2014). Secondly, we treat the particles and the inter-particle effective pores as spherical and quasi-monodisperse (we provide polydisperse solutions for *s* in the Data Repository). Despite these implicit assumptions, from the apparent efficacy of the present scaling approach, we are led to infer that, for the wide range of samples tested (Table 1), particle angularity or polydispersivity do not have a first-order effect on the inter-particle fluid permeability. The dominant effect of angularity might be on the initial porosity (i.e., the maximum packing density of the particles at deposition), which will in turn impact the permeability and is captured by the bulk ϕ in Equation 1.

Our model carries the implication that it is the processes that affect *R* (and thus *s*; Equations 2 and 3) during particle formation, such as the fragmentation conditions (Fowler et al., 2010), bulk composition, and crystal content (Rust and Cashman, 2011), that play an important role in the permeability of the resultant deposits and their evolution during diagenesis or welding. This inference holds the potential of providing an explicit physical link between the bulk isotropic permeability and the mechanisms that control the development of the microstructural elements that are deposited. We conclude that a well-constrained state of the initial isotropic granular system appears to be sufficient to estimate the evolution of isotropic permeability during densification by compaction of a sediment bed or by welding in volcanic interiors or ignimbrites (Fig. 3).

## CONCLUSIONS

Initially granular materials undergoing densification evolve in terms of microstructure. These microstructural changes affect the bulk porosity and the permeability. We present a universal scaling between a normalized isotropic permeability and the porosity, using relatively few, easily obtainable, measured parameters. We calibrate this scaling using well-characterized materials before estimating its applicability to the literature data available. This scaling provides a tool to predict the fluid permeability across a wide range of porosities, invaluable for fluid flow modeling during volcanic degassing through granular filled cracks and in evolving sedimentary basins.

We acknowledge support from the European Union’s seventh program for research and technology under grant agreement VUELCO (282759) and European Research Council grants EVOKES (247076) and SLiM (306488). Comments from Bob Holdsworth, Heather Wright, and two anonymous reviewers greatly improved the clarity of our work.

^{1}GSA Data Repository item 2016065, alternative scaling arguments, is available online at www.geosociety.org/pubs/ft2016.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.