An analytical form of the pore-emanating crack model in 3-D with implications for the uniaxial compressive strength of sandstone

SUMMARY

The uniaxial compressive strength |${\sigma _c}$| of rocks is a key material property in a wide range of applications. Models for |${\sigma _c}$| typically either require numerical solutions, restricting their wide utility, or are empirical and therefore confined to a specific case. Here, we study the theoretical pore-emanated crack model and provide an analytical emulator function that matches the 2-D and 3-D solutions to a high degree of accuracy over all porosities, |$\phi $|⁠. A key input to both the full solution and to our emulator functions is the pore radius, assumed in the model to be circular or spherical, in a porous rock. In most porous lithologies, including sandstone, the notion of a pore radius is poorly defined since they are built from compacted or lithified grains. And so here we explore statistical methods to find a characteristic pore length scale, |${l_2}$|⁠, from an initial particle radius; this method is provided as an easy-to-use supplementary tool. We advocate for the use of our 3-D function |${\sigma _c} \approx 1.57{K_{{\bf{Ic}}}}/( {{\phi ^{0.43}}\sqrt {\pi {l_2}} } )$|⁠, where |${K_{{\bf{Ic}}}}$| is the fracture toughness of the solid matrix. A compilation of |${K_{{\bf{Ic}}}}$| values for minerals and rocks allows us to explore the effect of this parameter and to make recommendations for appropriate values in the model. We compare our simple emulator function for |${\sigma _c}$| with existing data sets across a wide range of sandstones to demonstrate the utility of this law for applied cases. We find that our function performs particularly well for relatively low porosity sandstones (⁠|$\phi \mathbin{\lower.3ex\hbox{$\buildrel<\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 0.15$|⁠) representative of mature basin systems from a diagenetic point of view; we discuss alternative models that are more appropriate for higher porosity sandstones.

1 INTRODUCTION

Rock properties are key inputs to all models for rock deformation, crustal mechanics, sedimentary basin evolution, compaction and diagenesis (Churcher et al. 1991; Bjørlykke 2014; Zhang 2016; Heap et al. 2020). For this reason, substantial effort has been made to constrain the properties of rocks as a function of relatively simple rock characteristics such as porosity, grain size, pore geometry, cement and clay content and elastic constants (Bell 1978; Sammis & Ashby 1986; Cargill & Shakoor 1990; Zhang et al. 1990, 2020; Bell & Culshaw 1998; Palchik 1999; Kahraman 2001; Cuss et al. 2003a,b; Demarco et al. 2007; Zorlu et al. 2008; Zhu et al. 2010, 2011; Wong & Baud 2012; Duda & Renner 2013; Baud et al. 2014; Wasantha et al. 2015; Kim & Changani 2016; Chen et al. 2018; Tang et al. 2018; Heap et al. 2019; Carbillet et al. 2021, 2022, 2023; Noël et al. 2021; Hill et al. 2022; Qi et al. 2022). Among all rock properties, the strength of rocks is of particular interest because it is central to fault initiation and propagation (Ikari et al. 2011; Kilgore et al. 2012), localization of strain in the crust (Baud et al. 2006; Eichhubl et al. 2010; Wong & Baud 2012), drilling efforts (Onyia 1988; Shi et al. 2015), geological engineering including construction (Sonmez et al. 2003; Ramamurthy 2004; Perras & Diederichs 2014) and geothermal or oil reservoir stimulations (Zimmermann & Reinicke 2010; Davies et al. 2012), among other applications. Despite a very wide range of applications in both natural and industrial processes, rock strength models remain difficult to use, to implement or to validate against data. A key missing step is the pursuit of theoretically grounded analytical models that can be used easily to compare with data at a range of scales.

Models for the strength of rocks can be grouped into three main categories on the basis of model geometry [other types of micromechanical model are reviewed in Kemeny & Cook (1991)]. First, there are so-called ‘grain crushing models’ which assume a granular or particulate geometry, and are based on Hertzian contacts between grains (Zhang et al. 1990; Ergenzinger et al. 2011). Because most sedimentary rocks, or granular volcanic rocks such as tuff, start out with this granular packing geometry before being compacted or lithified, these grain crushing models are thought to be particularly suitable for high porosities of initial granular packings in sedimentary basins (Zhang et al. 1990; Wong & Baud 2012; Carbillet et al. 2023) or volcanic terrains (Zhu et al. 2011; Heap & Violay 2021). Secondly, there are so-called ‘pore-emanating crack models’ which assume that there are circular (2-D) or spherical (3-D) pores embedded in a solid matrix of rock (Sammis & Ashby 1986; Klein & Reuschlé 2003; Wong & Baud 2012; Baud et al. 2014; Vasseur et al. 2017; Heap & Violay 2021; Heap et al. 2021). The idealized geometry assumed in these pore-based models makes them most suitable for more mature (i.e. less porous) sedimentary settings in which diagenetic processes render the rock microstructure close to being approximated by spherical pores in a solid matrix. This is therefore a case typical of lower relative porosities compared with the grain crushing models (Carbillet et al. 2023) or, if in the case of volcanic lavas, the rock was never granular, then the pore-emanating crack model may be valid across most porosities (Heap & Violay 2021). Finally, the third category of model is a so-called ‘wing crack model’ in which a solid matrix is not considered porous except for starting microcracks which can connect up via growing crack ‘wings’ that propagate from the initial crack tips (Ashby & Sammis 1990; Baud et al. 1996, 2014; Brantut et al. 2012; Wong & Baud 2012; Bernabé & Pec 2022). These wing crack models are only applicable to dense and nominally non-porous rocks and are therefore typically applied to metamorphic or plutonic rocks such as schist or granite (Janach 1977; Wong & Baud 2012) or low-porosity lavas (Heap et al. 2014; Zhu et al. 2016; Heap & Violay 2021). Indeed, Zhu et al. (2016) found that the wing crack model overestimates the strength of basalts with only a porosity of about 0.05 by almost a factor of two, concluding that pore-emanating cracking is likely exerting great influence on their strength.

Here, we focus in on the pore-emanating crack model type because it has particularly wide utility spanning mature sedimentary settings as well as porous volcanic rocks such as lavas and welded ignimbrites (Zhu et al. 2011, 2016; Heap & Violay 2021). The best-known, and most widely used, pore-emanating crack model is the Sammis & Ashby (1986) model which starts from the assumption of circular or spherical pores in a solid matrix. We restrict our analysis of this model to the uniaxial case, simply because there is a wealth of experimental data against which we can compare solutions. However, we propose that once the model is validated for the uniaxial end-member case, the triaxial solutions could be used with more confidence. As will be explored here, this model typically requires a numerical solution to governing equations in order to predict a peak strength for a given stress state (Klein & Reuschlé 2003; Vasseur et al. 2017). However, Zhu et al. (2010) proposed an analytical functional form of the 2-D case of circular pores in a solid matrix that matches the full solution well over a wide range of porosities and input parameters. Here, we provide a simple-to-use, functional form for the 2-D and 3-D case, as well as providing a useful method to estimate the pore length scale and guidance on how to select an appropriate fracture toughness |${K_{{\rm{Ic}}}}$|⁠, vital input parameters in these analytical functional forms, for porous granular rock such as sandstone.

2 THE PORE-EMANATING CRACK MODEL (SAMMIS & ASHBY 1986)

The pore-emanating crack model starts from the geometric assumption that cylindrical (2-D) or spherical (3-D) gas-filled pores are embedded in a solid matrix. When loaded in compression, cracks grow from the pores. The extent of crack growth depends on radius of the pores, the confining pressure and the axial stress applied. The model accommodates the interaction of growing cracks with each other and with the sample edges. The model takes as a primary input the ratio of the minimum compressive stress |${\sigma _3}$| to the maximum compressive stress |${\sigma _1}$|⁠, where that ratio is |$\lambda = {\sigma _3}/{\sigma _1}$|⁠. Here, we focus on the uniaxial and unconfined case (sometimes called ‘simple compression’) where |${\sigma _3} = 0$|⁠, so that |$\lambda = 0$|⁠. The uniaxial strength is a widely used metric for the general strength of rocks and can be used over triaxial deformation metrics as a simple method to test the strength differences between rocks. By assuming |$\lambda = 0$|⁠, the resulting governing equations are simplified substantially. In this section, we lay out those governing equations of the model. Next, we give simple-to-use emulator functions that are analytical, and which match the full model across all practical conditions.

2.1 The full pore-emanating crack model

We consider a circular or spherical pore of radius a and a material with porosity |$\phi $|⁠. In 2-D, |$\phi $| is the area fraction of pores |${\phi _{2{\rm{D}}}} = {A_p}/{A_T}$| and, in 3-D, |$\phi $| is the volume fraction of pores |${\phi _{3{\rm{D}}}} = {V_p}/{V_T}$|⁠, where A and V represent areas and volumes, respectively, and subscript p or T represent the pore contribution and the total, respectively. When a far field stress |${\sigma _1}$| is applied onto a porous material in the absence of confining pressures (⁠|${\sigma _2} = {\sigma _3} = 0$|⁠), any pores in the material will concentrate stress on their surfaces (see Appendix). If |${\sigma _1}$| is axial, then the stresses are maximal and tensile at the pore poles (i.e. the top and bottom), leading to a tendency for cracks to appear and grow from the pore poles. Crack formation and growth occurs to a distance c away from the pore surface and cracks are subject to mode I loading with an associated stress intensity |${K_{\rm{I}}}$| that can be computed by integrating the stress field over the surface of the crack; the details of this integration are given in Sammis & Ashby (1986). The cracks grow until the stress intensity |${K_{\rm{I}}}$| becomes equal to the critical stress intensity |${K_{{\rm{Ic}}}}$| (also called the ‘fracture toughness’). In this case, following Sammis & Ashby (1986), the stable crack length c is related to the applied stress by

where |$\bar c = c/a$| is the dimensionless crack length, |$\bar \sigma = {\sigma _1}\sqrt {\pi a} /{K_{{\rm{Ic}}}}$| is the dimensionless applied stress, |${a_1}$| and |${a_2}$| are constants that depend on whether the 2-D or 3-D case is being considered, and |$\lambda = 0$| is assumed. The first term on the right-hand side of eq. (1) is the crack propagation term, while the second term is the crack interaction term (and hence depends on the pore number density via the porosity |$\phi $|⁠). In 2-D, |${a_1} = 1.1$| and |${a_2} = 3.3$|⁠, while in 3-D, |${a_1} = 0.62$| and |${a_2} = 4.1$|⁠; note these are derived quantities that arise from assumptions made by Sammis & Ashby (1986).

We note that in Sammis & Ashby (1986), the first term on the right-hand side of eq. (1) is strictly dependent on whether the 2-D or 3-D case is considered (via |${a_1}$| and |${a_2}$| that both depend on the dimensions), however the second term (the crack interaction term) appears to only be a 2-D solution and no equivalent 3-D solution is given for that term. Here, we use the crack interaction term as derived by Sammis & Ashby (1986) for all cases. Given that cracks form 2-D planes, we posit that it is possible that 2-D crack interaction is formally valid. Having said that, exploration of a 3-D crack interaction term in future would be a useful advance here.

Using eq. (1), we take the derivative |${\rm d}\bar \sigma /{\rm d}\bar c$| and look for when this derivative is equal to zero, which represents a peak stress achieved during loading. By doing this and rearranging, we find the porosity |$\phi $| that relates to the critical peak stress, and therefore to the peak crack length |${\bar c_c}$|

where |${\bar c_c} = {c_c}/a$| with |${c_c}$| the critical crack length at |${\rm d}\bar \sigma /{\rm d}\bar c = 0$| (see eq. 1). Eq. (2) as given here is a general form that can be adapted to 2-D or 3-D as needed and reduces to the 2-D form given by Zhu et al. (2010) when the appropriate |${a_1}$| and |${a_2}$| are used.

Eq. (2) can be solved numerically to find the maximum crack length |${c_c}$| for a given porosity |$\phi $| and then, by injecting that |${\bar c_c}$| in place of |$\bar c$| in eq. (1), gives the maximum stress |$\bar \sigma = {\bar \sigma _c}$|⁠. This maximum stress |${\sigma _c}$| can be thought of as the UCS. In the Appendix we give a description of some technical issues associated with non-uniqueness of |$\phi ( {{{\bar c}_c}} )$| (these technical issues do not affect the model output).

2.2 Analytical emulator functions in 2-D and 3-D

Using eqs (1) and (2), there is a curve that relates the dimensionless stress |${\bar \sigma _c} = {\sigma _c}\sqrt {\pi a} /{K_{{\rm{Ic}}}}$| to |$\phi $|⁠. This curve is universal in both the 2-D and 3-D form. In Fig. 1 we show that this curve has a power-law form such that |${\bar \sigma _c} = \alpha {\phi ^{ - \beta }}$| where |$\alpha $| and |$\beta $| are dimensionless constants such that |$\alpha $| controls the magnitude of |${\bar \sigma _c}$| at a reference |$\phi $| and |$\beta $| controls the slope |${\rm d}{\bar \sigma _c}/{\rm d}\phi $|⁠. We also note that the UCS is higher in the 3-D case than in the 2-D case, and that this difference is higher at lower porosity (Fig. 1). Based on their power-law form (Fig. 1), we propose the following emulator function

or dimensionally

We fit this emulator function (eqs 3 and 4) to the full solution across all |$\phi $| and find the best-fitting values of |$\alpha $| and |$\beta $| in 2-D and 3-D. This results in |$\alpha = 1.3345 \pm 0.0007$| and |$\beta = 0.4166 \pm 0.0002$| in 2-D; note that this is slightly different than the values given by Zhu et al. (2010). In 3-D, this results in |$\alpha = 1.5697 \pm 0.0012$| and |$\beta = 0.4297 \pm 0.0003$| (Wadsworth et al. 2022). Injecting these best-fitting values into eq. (4) gives us an analytical model for uniaxial compressive strength in 2-D as

and in 3-D as

It is important to note that the 2-D solution is based on exactly the same procedure as employed by Zhu et al. (2010) but the coefficients we obtain are slightly different, suggesting that there are differences in the accuracy of our respective numerical solutions to eq. (2). To provide an example as to the difference between the UCS values predicted in the 2-D and 3-D case, we assume here a |${K_{{\rm{Ic}}}}$| of 0.3 |${\rm{MPa}}.{{\rm{m}}^{1/2}}$|⁠, which is the |${K_{{\rm{Ic}}}}$| typically used for sandstones (Zhang et al. 1990; Wong & Baud 2012) and a pore radius a of 100 |$\mu\rm m$|⁠. For a porous rock with a porosity of 0.1, the 2-D case (eq. 5) provides a UCS of 58.9 MPa and the 3-D case (eq. 6) provides a UCS of 71. 5 MPa. Given that dimensionally this difference is substantial (i.e. +21 per cent in UCS going from 2-D to a proper 3-D assumption), we advocate for the 3-D solution (eq. 6) for general use with rock strength problems.

3 PREVIOUS WORK AND AIMS HEREIN

Constitutive models such as the pore-emanating crack model (Sammis & Ashby 1986) are designed to provide a predictive tool for the strength of porous solids. Deformation models for the crust therefore embed constitutive models in upscaled dynamics models. Despite their importance and wide utility, these models have scarcely been tested against experimental data sets for specific families of porous solids. Where data-driven model performance tests have been conducted, they typically stop short of making predictions of strength and instead show a range of possible model outputs by contouring for |${K_{{\rm{IC}}}}{( {\pi a} )^{ - 1/2}}$| (Heap & Violay 2021). This demonstrates that despite the development of predictive models that can indeed output |${\sigma _c}$| predictions for a given porosity |$\phi $|⁠, these predictions are rarely made concrete. This shortcoming arises from two key problems: (1) the value of |${K_{{\rm{IC}}}}$| that is appropriate for a given scenario is not well known; and (2) while grainsizes in sedimentary rocks are readily measurable, the meaning of the pore size a is also not well known. The first issue is related to the question of whether or not |${K_{{\rm{IC}}}}$| should represent the fracture toughness of the matrix of the grains, or of the cement, or whether a porous |${K_{{\rm{IC}}}}$| should be used in what would be akin to an effective medium approach (Zhu et al. 2011; Noël et al. 2021). The second issue is related to the question of how we define pore sizes in granular sedimentary rocks, and how pore size therefore relates to both grain size and porosity. Given these shortcomings, the principal aims of this contribution are to provide a solution or resolution to these two existing issues. This contribution therefore leads on from Vasseur et al. (2017) who tested a method to find pore sizes in granular media; however, herein we test this approach against a far larger experimental data set, yielding new insights. Similarly, we aim to conclude with an easy-to-use guide and quantitative toolkit to make predictions of strength for any sedimentary rock, and for sandstone in particular.

4 PREDICTING THE CHARACTERISTIC PORE SIZE IN A SEDIMENTARY ROCK

Eqs (5) and (6) are useful models for practical applications to uniaxial compressive strength data in rocks where a pore radius is a key input parameter (e.g. see eq. 6). However, pores are rarely spherical cavities in rocks, and especially granular rocks, and it is often unclear how to constrain the pore radius a beyond ballpark estimates based on 2-D images of rock thin sections. Other techniques such as mercury porosimetry provide the pore throat or pore access radius (Guéguen & Palciauskas 1994), which can differ greatly from the pore radius in some rocks, and more recently developed techniques such as X-ray computed tomography (Ketcham & Carlson 2001) requires access to expensive laboratory equipment. For rocks that form from solid grains, such as in sandstones, many limestones, tuffs and other granular lithologies, the grains are pseudo-spherical solids, the radius of which is typically easy to constrain, whereas the pores are a convolute interparticle fluid phase. And so, to facilitate the ease of estimating UCS using eqs (5) and (6) (Section 2), here we develop tools to derive the relevant pore radius a using a measured particle size or particle size distribution.

Think of a model rock as a pack of spherical grains with radius R. Here, we envisage that compaction and diagenetic processes can be modelled geometrically and approximately as interpenetration of those spherical grains, reducing the bulk porosity below the initial packing porosity (Wadsworth et al. 2016; Vasseur et al. 2017). The so-called nearest neighbour function can be used to give the average interpore length. We can use the probability |$F( r ){\rm d}r$| that a sphere centre lies at a distance between r and |$r + {\rm d}r$| to find the nth moment of |$F( r )$| given by (Torquato et al. 1990)

where a bar above a symbol denotes a parameter normalized by the sphere radius R [i.e. |$\bar r = r/R$| and |$\bar F( {\bar r} ) = F( r )R$|]. The first moment (i.e. |$n = 1$|⁠) gives the mean nearest neighbour distance between spheres |$\bar l \equiv \langle {\bar r} \rangle $|⁠. In our case the spheres of radius R represent the particles and we have an interpore distance termed |${l_1} = {\bar l_1}R$| and an interparticle distance termed |${l_2} = {\bar l_2}R$|⁠. For monodisperse and fully interpenetrable spheres, the nearest neighbour function is (Torquato et al. 1990)

where |$\eta $| represents the sphere reduced density (i.e. the product of sphere number density and an individual sphere volume). Combining eq. (8) with eq. (7), and taking |$n = 1$|⁠, results in an analytical expression

where |${\rm{\Gamma }}$| is the so-called gamma function [available generally as a built-in function in all recent versions of Excel™ as ‘=GAMMA(input)’, where ‘input’ is then 4/3 in our case] which here takes the input 4/3 (returning a value of |$\sim 0.89$|⁠), and |$\eta = - \ln ( {1 - \phi } )$| when |$i = 1$| (the case when |${\bar l_1} = {l_1}/R$|⁠) and |$\eta = - \ln \phi $| when |$i = 2$| (the case when |${\bar l_2} = {l_2}/R$|⁠). The distance |${l_1}$| can be thought of as an interpore distance and the distance |${l_2}$| can be thought of as an interparticle distance. Vasseur et al. (2017) suggested that the latter is a reasonable proxy for an average pore diameter in heterogeneous granular media such as sandstones. In Fig. 2 we plot these two distances normalized by the particle radius and how they vary with porosity for the case of overlapping monodisperse spheres.

For comparison, in Fig. 2 we additionally plot the Lu & Torquato (1992) solution for a characteristic pore radius between overlapping particles |$P( a )$| (here |$\eta = - \ln \phi $|⁠)

where |$\bar a = a/R$| and |$\bar P( {\bar a} ) = P( a )R$|⁠. The nth moment of |$P( a )$| is given by

and the first moment (i.e. |$n = 1$|⁠) gives the mean pore radius |$\langle a \rangle = \langle {\bar a} \rangle R$|⁠.

We have three possible metrics for a ‘pore size’ developed in this section: (1) the distance between any two pores |${l_1}$|⁠, (2) the distance between any two particles |${l_2}$| and (3) the mean size of a spherical pore between a random heterogeneous array of overlapping particles |$\langle a \rangle $|⁠. The distance between any two pores |${l_1}$| is a negative function of porosity, such that it is large at low porosity and diminishes as |$\phi \to 1$| (Fig. 2). This length |${l_1}$| is likely to be an appropriate measure of the effective grain size, which becomes larger as grains are compacted together to low porosity and can exceed the original grain size. Therefore, while our formalism allows us to give that length scale here, it is unlikely to be appropriate for representing pores. Instead, the expectation should be that as porosity increases, the pore size increases. That expected behaviour is the case for both |${l_2}$| and |$\langle a \rangle $| (Fig. 2). Therefore, we anticipate that |${l_2}$| and |$\langle a \rangle $| are viable candidates for a pore radius a in the context of the pore-emanating crack model (eqs 5 and 6; Section 2). In the Supplementary Information we provide an Excel™ implementation of these three metrics so that a user can output them using a measured grain radius R or a measured grain size distribution for their granular rock.

5 THE FRACTURE TOUGHNESS |${K_{{\rm{IC}}}}$|

Inspection of our emulator function form (eq. 4) shows that along with a pore length scale (explored in Section 4), the other parameter that requires constraint is the fracture toughness of the solid matrix |${K_{{\rm{Ic}}}}.$| This parameter is well known for being difficult to measure directly (Barker 1977; Zoback 1978; Bergkvist & Fornerod 1979). Here, we compile existing data sets for |${K_{{\rm{Ic}}}}$|⁠, restricting ourselves to (1) |${K_{{\rm{Ic}}}}$| for minerals that are pertinent to sandstones (quartz, calcite and feldspar), and (2) |${K_{{\rm{Ic}}}}$| for natural bulk sandstones. Therefore, we omit carbonate, metamorphic and igneous rocks. This results in a data set for quartz, calcite and feldspar minerals (Wiederhorn 1969, 1974; Hartley & Wilshaw 1973; Barker 1977, 1979b; Atkinson 1979a; Atkinson & Avdis 1980; Norton & Atkinson 1981; Meredith & Atkinson 1982) and for a wide range of sandstones (Brown et al. 1972; Clifton et al. 1976; Suzuki et al. 1978; Zoback 1978; Bergkvist & Fornerod 1979; Atkinson 1979b, 1980, 1984; Rummel et al. 1980, 1985; Dibb et al. 1983; Rummel & Winter 1983; Winter 1983; Gunsallus & Kulhawy 1984; Meredith et al. 1984; Müller 1984; Senseny & Pfeifle 1984; Atkinson et al. 1985; Chandler et al. 2017; Noël et al. 2021). These data are compiled in Table 1.

The data set compiled here (Table 1) shows that the matrix minerals pertinent to sandstones have |${K_{{\rm{Ic}}}}$| values that vary between 0.19 |${\rm{MPa}}.{{\rm{m}}^{1/2}}$| for calcite (⁠|$10\bar 11$| orientation) and 2.4 |${\rm{MPa}}.{{\rm{m}}^{1/2}}$| for quartz normal to the c-axis. Within the quartz measurements (not including fused quartz) the variability is similarly large varying between 0.31 |${\rm{MPa}}.{{\rm{m}}^{1/2}}$| (⁠|$0001$| orientation) and that same upper limit of 2.4 |${\rm{MPa}}.{{\rm{m}}^{1/2}}$|⁠. This measurement compilation suggests that the widely used value for sandstone matrix solid of 0.3 |${\rm{MPa}}.{{\rm{m}}^{1/2}}$| (Zhang et al. 1990; Baud et al. 2014) may require careful inspection. In the context of the pore-emanating crack model investigated herein, it is clear that it is the mineral (solid matrix) value of |${K_{{\rm{Ic}}}}$| that should be used in eq. (4) and indeed in any of the solutions described in Section 2. However, it is also common to find values for bulk porous rocks such as sandstones (Table 1). When compiled, these data show a strong correlation between |${K_{{\rm{Ic}}}}$| and porosity |$\phi $|⁠, where |${K_{{\rm{Ic}}}}$| diminishes as porosity increases (Fig. 3). Fig. 3 suggests that the widely used value of 0.3 |${\rm{MPa}}.{{\rm{m}}^{1/2}}$| is consistent with a relatively high porosity rock, rather than necessarily representing the matrix value of the solid between the pores from which cracks are emanating during deformation. The increase in |${K_{{\rm{Ic}}}}$| towards |$\phi = 0$| suggests that an extrapolation of the general |${K_{{\rm{Ic}}}}( \phi )$| trend could indicate a matrix solid value that is reliable (which would be consistent with the observation that some of the quartz |${K_{{\rm{Ic}}}}$| values are far higher than 0.3 |${\rm{MPa}}.{{\rm{m}}^{1/2}}$|⁠). To find a useable general value for |${K_{{\rm{Ic}}}}$|⁠, we take the trend of |${K_{{\rm{Ic}}}}( \phi )$| and fit a linear equation |${K_{{\rm{Ic}}}} = {k_1} - {k_2}\phi $| and an exponential equation |${K_{{\rm{Ic}}}} = {k_3}\exp ( { - {k_4}\phi } ).$| For these two functional forms, |${k_1}$| and |${k_3}\,\,$| represent the intercepts for |${K_{{\rm{Ic}}}}$| when |$\phi = 0$|⁠. Our minimization suggests a best-fitting linear form has |${k_1} = 1.58 \pm 0.15\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| and |${k_2} = 4.92 \pm 0.95$|⁠, and a best-fitting exponential form has |${k_3} = 2.15 \pm 0.20\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| and |${k_4} = 6.78 \pm 0.79$| (Fig. 3). The values |$1.58\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| and |$2.15\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| represent the extrapolated matrix solid values when the porosity is zero, and they are consistent with the higher values of the quartz measurements compiled here (Table 1). If these extrapolated values are to be used, then it suggests that the highest |${K_{{\rm{Ic}}}}$| values in a given solid are the limiting values that must be met in fracture experiments or natural fracture scenarios. The most effective |${K_{{\rm{Ic}}}}$| value to use when modelling the strength of rocks is discussed later and a rationale for the result is explored.

6 COMPARISON WITH DATA FOR THE UNIAXIAL COMPRESSIVE STRENGTH OF SANDSTONES

To test the analytical emulator functions for the uniaxial compressive strength (eqs 5 and 6), we compile a data base of laboratory results for the uniaxial compressive strength of sandstones (Bell 1978; Bell & Culshaw 1998; Palchik 1999; Cuss et al. 2003a,b; Demarco et al. 2007; Duda & Renner 2013; Baud et al. 2014; Wasantha et al. 2015; Kim & Changani 2016; Chen et al. 2018; Tang et al. 2018; Heap et al. 2019; Zhang et al. 2020; Noël et al. 2021; Hill et al. 2022; Qi et al. 2022) and synthetic rocks that are designed to represent sandstones (Vasseur et al. 2013, 2017; Fattahpour et al. 2014; Atapour & Mortazavi 2018; Rice-Birchall et al. 2021; Shahsavari & Shakiba 2022; Carbillet et al. 2023). This data set is restricted to room temperature, dry (i.e. not conducted in the presence of water), and all performed at, or close to, the standard axial strain rate for laboratory experiments of 10⁻⁵ s⁻¹. The advantage of the synthetic rocks is that in most cases the grain sizes are known precisely, whereas for the natural sandstones, they are, for the most part, crudely measured from images of 2-D thin sections. We restrict ourselves to data for which the grain size of the natural sandstone is constrained or estimated in the publication containing the data. In Table 2 we give a summary of the references from which these data are taken, and in the Supplementary Information we provide the full data set.

The data compiled here in Table 2 are for the critical peak value of |${\sigma _1}$| at which the samples failed to be load-bearing in uniaxial compression. This is taken to be the uniaxial compressive strength |${\sigma _c}$|⁠. First, we plot the data as the measured |${\sigma _c}$| value (in MPa) as a function of the sample |$\phi $|⁠, which represents the raw data tested herein (Fig. 4). Because each sample has a different value of mean grain radius R (Table 2), it is useful to normalize the measured |${\sigma _c}$| values using |${\bar \sigma _c} = {\sigma _c}\sqrt {\pi \mathcal{L}} /{K_{{\rm{Ic}}}}$| where we pick |$\mathcal{L} = {l_2}$| (Figs 5a and b; using eq. 9) or |$\mathcal{L} = \langle a \rangle $| (Fig. 5c and d; using eq. 11). In Fig. 5 we also plot the emulator functions in 2-D (eq. 5; blue curves) and 3-D (eq. 6; red curves) which are analytical and do not require any minimization to the data. Finally, we also use different values for the |${K_{{\rm{Ic}}}}$| for the natural sandstones: Figs 5(a) and (c) use the standard value of |${K_{{\rm{Ic}}}} = 0.3\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| for natural sandstones (Zhang et al. 1990; Wong & Baud 2012) and Figs 5(b) and (d) use |${K_{{\rm{Ic}}}} = 2.15\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$|⁠, which is the limiting intercept of the exponential fit in Fig. 3. For the synthetic sandstones studied here, for which glass is the main constituent, we use a value of |${K_{{\rm{Ic}}}} = 0.7\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| in Figs 5(a)–(d) (Wiederhorn 1969; Vasseur et al. 2013, 2017). By comparing Fig. 5(a) with Fig. 5(c) or Fig. 5(b) with Fig. 5(d), it is clear that the interparticle length |${l_2}$| (i.e. Figs 5a and b) is the superior descriptor of the pore size in all of the studied samples and that the pore emanating crack model is robust without the need for fitting or minimization to data. Indeed, when using |$\mathcal{L} = \langle a \rangle $| (Figs 5c and d), |${\bar \sigma _c}$| is lowered such that they are far below the model curves. In Fig. 3 we showed that there is evidence that |${K_{{\rm{Ic}}}} = 0.3\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| is a low value for the matrix solid in sandstones. However, using |${K_{{\rm{Ic}}}} = 2.15\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| has the effect of pushing all of the data down below the model curves (Fig. 5) and for this reason, we suggest that this higher solid value of |${K_{{\rm{Ic}}}}$| is not useful when using pore-emanating crack models (discussed further in Section 6).

Finally, we find that, visually, it is clear that the model is a more reasonable description of the uniaxial strength for the lower end of the porosity spectrum (⁠|$\phi \mathbin{\lower.3ex\hbox{$\buildrel<\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} \,\,0.15)$| studied here (Figs 5a and c). This can be seen most clearly when examining the inset to Fig. 5(a). This is discussed in more detail later. To summarize, when |${K_{{\rm{Ic}}}} = 0.3\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| is assumed and if we restrict ourselves to |$\phi \mathbin{\lower.3ex\hbox{$\buildrel<\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} \,\,0.15$| then the 3-D model (eq. 6) fits the data, albeit with substantial spread either side of the model prediction (red curve in Fig. 5(a); discussed later).

7 DISCUSSION

Our results (Fig. 5) show that the pore-emanating crack model emulator functions found here perform satisfactorily when compared with the uniaxial compressive strength of sandstones if the pore length scale is taken to be the calculated distance between two neighbouring particles |${l_2}$| and if the fracture toughness is taken to be |${K_{{\rm{Ic}}}} = 0.3\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| (see Fig. 5a). However, it is also clear that there are remaining issues that warrant discussion. First, at porosities |$\phi \mathbin{\lower.3ex\hbox{$\buildrel>\over {\smash{\scriptstyle\sim}\vphantom{_x}}$}} 0.15$| the data deviate substantially from the model prediction where the prediction overestimates the uniaxial strength. Secondly, even though the fracture toughness required to predict the data has been recommended previously (Zhang et al. 1990), it is substantially lower than what might be expected for the solid matrix in sandstones (Fig. 3). And thirdly, the best performance of the model still has large residuals between model and experiment (Fig. 5a) such that the model is accurate but exhibits poor precision. Here, we discuss these points in turn.

7.1 A ‘crossover porosity’

Carbillet et al. (2023) suggest that there is a critical porosity (or a ‘crossover porosity’) |$\phi ^{\prime}$| below which the pore-emanating crack model holds validity but above which it does not. They present direct microtextural evidence that suggest that for |$\phi > \phi ^{\prime}$|⁠, the micromechanical damage mechanism is one in which individual grains in a sandstone are internally crushed, rather than one in which cracks propagate from pores. Cracks propagating from pores is the crack geometry/mechanism on which the pore-emanating crack model is based and this is only applicable for relatively low porosities |$\phi < \phi ^{\prime}$|⁠. That concept of a microstructural transition from granular model validity at higher relative |$\phi $| to pore-based models at lower relative |$\phi $| has been referred to as a topological inversion of the geometry (Wadsworth et al. 2017). At |$\phi > \phi ^{\prime}$|⁠, the geometry is one in which grains are effectively spheres embedded in a convolute and interconnected intergrain gas phase, whereas at |$\phi < \phi ^{\prime}$|⁠, the geometry is one in which spherical pores are effectively embedded in an interconnected interpore solid phase (Wadsworth et al. 2017).

Our compiled data here seem to suggest that |$\phi ^{\prime} \approx 0.15\,\,$| (Fig. 5). To test this more quantitatively, we take two important steps. First, we investigate the ‘grain crushing model’ (Zhang et al. 1990) that is purported to provide superior accuracy at high porosity compared with the pore-emanating crack model. When normalized using |${\bar \sigma _c}$|⁠, the grain crushing model is

where |$\xi \approx 86$| is a constant related to the Poisson's ratio (Carbillet et al. 2023), E is the Young's modulus of the matrix and |${\bar c_i} = {c_i}/R$| is the initial crack length |${c_i}$| normalized by the particle radius. It is clear that the normalization of the grain crushing model is inconsistent because, unlike with the pore-emanating crack model, even after normalization dimensional material property parameters are retained (eq. 12). Nevertheless, eq. (12) is in a form that can be graphed and compared with our compiled data set.

Carbillet et al. (2023) provide a thorough analysis of eq. (12) with respect to experimental data, which we do not attempt to reproduce here. However, taking |$E = {10^{10}}\,\,{\rm{Pa}},$||${K_{{\rm{Ic}}}} = 0.3\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$|⁠, and an indicative |$R = 100\,\,{\rm{\mu m}}$|⁠, we can plot eq. (12) along with our pore emanating crack model function (Fig. 6). To do so, we choose two bounding values of |$3 \times {10^{ - 10}} \le {\bar c_i} \le 5 \times {10^{ - 10}}$|⁠. It is clear that |${\bar c_i}$| has an influence, but we also acknowledge that |${\bar c_i}$| is wholly unconstrained and therefore is simply used as an adjustable parameter for illustrative purposes. Fig. 6 shows that the crossover porosity |$\phi ^{\prime}\,\,$| at which the two models—grain crushing and pore-emanating cracks—meet is |$\phi ^{\prime} \approx 0.125$|⁠, close to the |$\phi ^{\prime} = 0.15$| identified by eye.

Although this contribution is a meta-analysis and so we do not explore microstructures in post-experimental samples directly (post-experimental microstructure is extremely rare in the studies that contribute to our data base), the fact that the pore-emanating crack model overpredicts the observed uniaxial compressive strength at |$\phi > \phi ^{\prime}$| is consistent with the regime crossover predicted by Carbillet et al. (2023). Those authors advocate for the use of the so-called ‘grain-crushing models’ at higher porosities |$\phi > \phi ^{\prime}$| (Zhang et al. 1990) which tend to predict lower uniaxial compressive strengths when compared with pore-emanating crack models at higher relative porosities.

7.2 Low fracture toughness

The second point worth discussion is the validity of a relatively low fracture toughness |${K_{{\rm{Ic}}}}$| for pore-emanating crack models (Fig. 5a) when the microcracks propagate through the solid matrix in sandstones, very often quartz, the solid matrix should have a relatively high fracture toughness (Fig. 3). The value advocated for here of |${K_{{\rm{Ic}}}} = 0.3\,\,{\rm{MPa}}.{{\rm{m}}^{1/2{\rm{\,\,}}}}$| for sandstones would appear to be consistent with the fracture toughness of relative high porosity sandstone (Fig. 3). At first, this might lead us to conclude that the fractures emanating from pores are propagating through porous matrix in what could be considered a nested effective medium model (Zhang et al. 1990). But this is inconsistent with microstructural analyses on natural and synthetic sandstones which, at |$\phi < \phi ^{\prime}$|⁠, do indeed show that fractures emanating from pores are propagating through the solid matrix and not through multiple pore-solid junctions as if the matrix were an effective porous medium (Carbillet et al. 2023). A key observation is that the uniaxial compressive strength of variably sintered glass beads appears to be well-predicted by our 3-D emulator function at |$\phi < \phi ^{\prime}$| when |${K_{{\rm{Ic}}}} = 0.7\,\,{\rm{MPa}}.{{\rm{m}}^{1/2{\rm{\,\,}}}}\,\,$| is used (Fig. 7a). The value |${K_{{\rm{Ic}}}} = 0.7\,\,{\rm{MPa}}.{{\rm{m}}^{1/2{\rm{\,\,}}}}$| is the valid value for glass (Wiederhorn 1969), which is the solid matrix constituent of the variably sintered glass bead samples. An advantage of those samples is that there is no other material involved in their microstructure other than that glass. By contrast, sandstones have a principal solid matrix phase which is often quartz-dominated, however they also bear an intergrain cement that can be substantially weaker. It is not clear if the cracks exploit weak cement preferentially. Based on these considerations, we have to conclude that the weaker cement material and/or material with substantially lower |${K_{{\rm{Ic}}}}$| than that of intact quartz is involved in the fracturing that emanates from pores at |$\phi < \phi ^{\prime}$|⁠. Therefore, there is an opportunity to extend our understanding of fracture toughness of cement materials.

7.3 Model performance and sandstone heterogeneity

Finally, it is clear from our principal model comparison (Fig. 5) that even in the best case of using |${l_2}$| as the pore length scale and using |${K_{{\rm{Ic}}}} = 0.3\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| (discussed above), there is substantial spread of the data either side of the model prediction. In Fig. 5(a) this spread can be up to |$\sim 10\,\,{\rm{MPa}}$|⁠. Given that this is a large compilation of sandstone data from sandstones worldwide, it is likely that there would be variability in the sandstone mineralogical composition, which could lead to mineralogical controls on the |${\sigma _c}$| that are not accounted for by the pore-crack model. To account for them, the only parameter that would be appropriate to adjust would be |${K_{{\rm{Ic}}}}$|⁠, which is clearly variable from mineral to mineral (Fig. 3; Table 1) so it is possible that a model with a |${K_{{\rm{Ic}}}}$| specific to each sandstone mineralogy could be useful for accounting for the variability in strength. However, given that |${K_{{\rm{Ic}}}}$| for many sandstones is much lower than for the matrix minerals (e.g. quartz; Figs 3 and 7), there is no clear method by which |${K_{{\rm{Ic}}}}$| could be computed for a given mineral mixture. Therefore, a useful test here is to filter the data for quartz content, retaining only the sample data for sandstones with relatively high quartz. Useful filters would be to retain only quartz arenites (⁠|$\ge 95\,\rm per\, cent$| quartz), or sub-arkose (⁠|$> 75\,\rm per\, cent$| quartz) and sub-litharenite (⁠|$> 75\,\rm per\, cent$| quartz). However, filtering for only quartz arenites results in very few data remaining, and filtering for all sub-arkose results in the majority of the data set remaining. Therefore we take a threshold |$> 85\,\rm per\, cent$| quartz, acknowledging that this is similar to taking all sub-arkose sandstones and excluding those with low quartz content. In Fig. 8 we show the result of this approach by presenting a side-by-side comparison of the model applied to all data (Fig. 8a) and the filtered data set (Fig. 8b). With that in mind, it is clear that high quartz content samples are better predicted by the model than low quartz content sandstones.

8 SOME LIMITATIONS AND OUTLOOK

In this contribution, we have tested the uniaxial limit of the Sammis & Ashby (1986) model for compressive strengths of porous rocks. Here, we state some of the key limitations that we have run up against, and an outlook for future work.

1. As identified in Section 2, the crack interaction term derived by Sammis & Ashby (1986) appears to be restricted to the 2-D case, whereas the crack propagation term is specifically given for both the 2-D and 3-D cases. Not having a formally 3-D version of a crack interaction term could be a current shortcoming in solving the 3-D case and the difference between crack interactions in 2-D and 3-D is currently unknown.

2. While the Sammis & Ashby (1986) model is triaxial, here we have chosen to focus on validating the uniaxial case only. This is an expedience because there is a wealth of available data for the uniaxial compressive strength of sandstones. Future work should validate the Sammis & Ashby (1986) model across a range of triaxial stress states.

3. In deriving the pore length scale metric |${l_2}$|⁠, we have used solutions designed for random heterogeneous media, which do not assume a geometry of the ‘pore’ that is being captured by that length scale. However, the Sammis & Ashby (1986) model does explicitly assume a spherical pore—or cavity—in the solid matrix. In natural sandstones, while |${l_2}$| may be a robust measure of the pore space ‘size’, it remains clear that the spherical pore assumption by Sammis & Ashby (1986) is violated by natural samples. Here we argue that low-porosity sandstones will host pores that are closer to the spherical pore assumption. In future, it could be useful to derive stress intensity functions for non-spherical pores closer to the cuspate forms found in natural sandstones across a range of porosity.

9 CONCLUSIONS

We develop analytical approximations to the pore-crack model by Sammis & Ashby (1986) for the 2-D and 3-D case of porous material failure in simple (uniaxial) compression. We show that the models match the full numerical solution across all porosities of interest, and that they retain the physical parameters of fracture toughness and pore length scale that are required to match real data. We provide a downloadable tool to find pore length scales in granular rocks via statistical methods developed for random heterogeneous granular media. We show that the fracture toughness of sandstones and their constituent minerals is highly variable, pointing to substantial uncertainty in this value. By compiling a global data base of sandstone strength data, we are able to select pore length scale choice |${l_2}$| (the average distance between any two grains) and a seemingly universal |${K_{{\rm{Ic}}}} = 0.3\,\,{\rm{MPa}}.{{\rm{m}}^{1/2}}$| for the fracture toughness. These choices cause the model to predict the bulk of the data without any fitting or empirical adjustment. Finally, we breakdown the data base into smaller subsets such as individual data sets that are particularly large or via filtering for high-quartz arenites, in order to understand better what causes the data to spread substantially either side of the model predictions. We find that mineralogy complexity is likely to be a primary cause of variability. In Fig. 9 we give a work flow in the form of a flow chart to help the reader navigate the principal points of our model and choices involved in using the tools we provide.

ACKNOWLEDGEMENTS

We are grateful to Donald B. Dingwell for supporting JV's contribution to this work via ERC Advanced Grant EAVESDROP (no. 834225) and to Lucille Carbillet and Patrick Baud for stimulating discussion relating to Carbillet et al. (2023) and Lucille Carbillet's PhD work in general. We thank M. Ashby for clarifications made about Sammis & Ashby (1986) by email. FBW thanks the staff at Station C3 in the Nußbaumstr. 7 psychiatric facility in Munich for supporting his recovery and creating an environment in which this article could be drafted. MJH acknowledges support from the Institut Universitaire de France (IUF).

AUTHOR CONTRIBUTIONS

FBW and JV conceptualized the study and performed the analysis of the pore-emanating crack model. MJH compiled the published data for the uniaxial compressive strength and provided rock mechanical context throughout. All authors contributed to the writing of the manuscript.

DATA AVAILABILITY

All data used herein and the quantitative tool have been published with figshare (Wadsworth 2024a,b). Additionally, those same resources are given as Supplementary Material available with this article.

REFERENCES

APPENDIX 1: NON-UNIQUENESS IN SOLUTIONS TO THE PORE-EMANATING CRACK MODEL (SAMMIS & ASHBY 1986)

The pore-emanating crack model predicts the porosity |$\phi $| that relates to a given critical crack length |${\bar c_c}$| (eq. 2). However, the full numerical solution to eq. (2) shows that a single porosity relates to up to three different values of critical crack length (see Fig. A1). This problem is not discussed by Sammis & Ashby (1986) nor by Zhu et al. (2010) who find the numerical solution to eq. (2) for the 2-D case. Nevertheless, cross-referencing the result for the strength (see Fig. 1) with the original Sammis & Ashby (1986) strongly suggests that a criterion that we only consider |${\bar c_c} \ge 1$| is valid, thereby restricting ourselves to the right hand branch of Fig. A1. This also avoids the associated issue that some of the output porosity values are unphysical at |$\phi > 1$|⁠. Fig. A1 shows that there is no 3-D solution all the way up to |$\phi = 1$|⁠, which thereby limits the use of our emulator functions (eq. 6) to lower values of porosity; that is consistent with arguments we make for the use of the pore-emanating crack model exclusively for porosities |$\phi < \phi ^{\prime}.$|