Stress-Induced Anisotropic Poroelasticity in Westerly Granite

We measured poroelastic properties of cracked granite under triaxial conditions, at elevated confining pressure and a range of differential stresses. Skempton's coefficients and undrained Young's modulus and Poisson's ratio were determined directly by recording in situ fluid pressure during rapid cycles of axial and radial stress. Drained properties were measured both statically and dynamically at ultrasonic frequencies. At a given confining pressure, increasing differential stress leads to the development of elastically transverse isotropy, with symmetry axis aligned with the compression axis. Skempton's coefficients


10.1029/2023JB026909
2 of 18 rocks occurring in fault damage zones can be anisotropic due to pre-existing or stress-induced preferred crack orientations (Rempe et al., 2013).Under these conditions, the Skempton coefficient becomes a tensor.However, full characterization of undrained poroelastic quantities is challenging (Cheng, 2021), and only scarce measurements of the Skempton tensor exist in the literature (Lockner & Beeler, 2003;Makhnenko & Tarokh, 2018).
The Skempton coefficients capture the pore pressure response to stress under undrained conditions expressed by  = 1 3  , where p is pore pressure, and B ij and σ ij are the Skempton coefficient and stress in the i, j direction, respectively.One of the main challenges when measuring poroelastic properties in the laboratory, and specifically undrained parameters such as the Skempton tensor, is that monitoring pore pressure under truly undrained conditions is difficult to achieve.Pore pressure measurements require access to the pore space, which typically disturbs the poroelastic response of the material, and special care is needed to minimize this disturbance.There has been extensive experimental work that produced measurements of isotropic poroelastic parameters (e.g., Berge et al., 1993;Green & Wang, 1986;Hart & Wang, 1995, 2001), where only four independent quantities are sufficient to fully describe the material.For a general anisotropic material the number of independent poroelastic parameters increases to 28.In the special case of transverse isotropy, relevant to cracked rocks under conventional triaxial stress conditions, the number of independent poroelastic parameters reduces to eight.Aoki et al. (1993) produced results for seven of those eight parameters.Using the theoretical framework of Cheng (1997), which provides a complete set of relationships between all anisotropic poroelastic parameters, Lockner and Stanchits (2002) and Lockner and Beeler (2003) measured both undrained and drained transversely isotropic parameters in Berea sandstone, achieving measurements for the two independent principal components of the Skempton tensor.Despite this early success, more systematic measurements remain to be obtained in other rock types, especially in cracked crystalline rocks like granite.
To complement direct experimental measurements, theoretical approaches can be used to determine the exact role of microstructure in the overall effective poroelastic properties of rocks.Wong (2017) proposed a theory where the components of the Skempton tensor are expressed simply as a combination of the elastic moduli of the solid matrix, and two principal components of the crack density tensor defined by Sayers and Kachanov (1995).By combining this approach with independent measurements of crack densities (e.g., by means of wavespeed measurements, see for instance Sayers and Kachanov (1995)), poroelastic properties could be estimated without need of direct laboratory measurements.However, this theory has not yet been thoroughly tested.
Here we experimentally measure drained, anisotropic elastic moduli in cracked granite, and make predictions of undrained quantities, including the Skempton tensor.Those predictions are tested against direct independent measurements of the Skempton tensor.This allows us to directly test the theory developed by Wong (2017).As our test material, we choose Westerly granite, which is to a good approximation elastically isotropic, in the undeformed state.We introduce an initial isotropic network of thermal microcracks, and control the crack density and anisotropy by changing the macroscopic stress state, which induces closing and opening of cracks along the principal stress orientations.The full transversely isotropic stiffness tensor is estimated under dry (equivalent to drained) conditions via ultrasonic wave velocity measurements, from which we also invert for the crack density tensors.Skempton coefficients are measured using in situ, low volume pore pressure transducers (Brantut, 2020;Brantut & Aben, 2021) following step changes in load and confining pressure.A limitation in the method is that we are comparing dynamic and static elastic moduli which are known not to be the same (Paterson & Wong, 2005); despite this approximation, we show that the main trends observed in direct measurements are well captured by predictions based on crack density tensors.

Theoretical Background
In this Section, we briefly summarize key poroelastic relations both in terms of macroscopic (bulk) and micromechanical parameters.The micromechanical relations are also given in terms of the second-order crack density tensor using the formulation of Wong (2017).

Transversely Isotropic Poroelasticity
Based on Biot's theory, the strain-stress relations for a transversely isotropic poroelastic material can be expressed (following Cheng (1997)), with compressive stresses positive, as 10.1029/2023JB026909 3 of 18 where ɛ is the strain tensor, S is the fourth rank compliance matrix, σ is the stress tensor, β σ is the constant stress storage coefficient (units Pa −1 ; denoted by C in Cheng (1997) and by S in Wong ( 2017)), B T is the transpose of the Skempton tensor, p is the pore pressure, and ζ is variation of fluid content (volume of fluid entering solid frame per unit volume of solid frame, units m 3 /m 3 ).In the most general case there are 28 independent poroelastic parameters describing the material: 21 components of S, six components of B, and β σ .For transversely isotropic materials, a set of only eight independent bulk parameters are needed, for instance, where the symmetry axis is the vertical direction (index 3 or z axis): {S 11 , S 12 , S 13 , S 33 , S 44 , B x , B z , β σ } (in Voigt notation).Here, we assumed vertical transverse isotropy, where the symmetry axis is taken as the compression axis.
We follow the micromechanical analysis of Cheng (1997) and make the assumptions of micro-homogeneity and micro-isotropy of the solid skeleton, which allows us to relate bulk properties to properties of the solid and fluid constituents.Specifically, the radial (B x ) and axial (B z ) Skempton coefficients can be expressed respectively as and where S ij are the elements of the compliance matrix (expressed in Voigt notation), and K 0 is the bulk modulus of the solid.The expression for storage is where K is the bulk modulus of the porous material, K fl is the bulk modulus of the saturating fluid, and ϕ is the porosity.The bulk modulus K is expressed in terms of the compliances By expressing the bulk modulus in terms of compliances the transversely isotropic Skempton coefficients are given by and

Relationship to Microstructure
The presence of microcracks increases the compliance of rocks, and a preferential distribution of these cracks leads to an anisotropic elastic response to stress.The effect of microcracks on macroscopic, average elastic properties of materials can be computed using effective medium models.In an homogeneous isotropic elastic material containing a population of noninteracting, penny-shaped microcracks with a transversely isotropic distribution of crack orientations, the effective compliance matrix S ij can be expressed in terms of isotropic solid material Poisson's ratio (ν 0 ) and Young's modulus (E 0 ) and the five independent components of the crack density tensors α and β: α 11 , α 33 , β 1111 , β 1133 , and β 3333 (Sayers & Kachanov, 1995).

10.1029/2023JB026909
4 of 18 The extra compliance due to cracks in rock is given by Sayers and Kachanov (1995) (their Equation 9) where B T and B N are the normal and shear components of the crack compliance tensor, n i is the ith component of the crack normal, V is the volume element containing cracks with surfaces S r .
For a dry circular crack, with radius a, and The change in compliance due to cracks is then where the second rank crack density tensor α ij is defined by where a (m) is the mth crack radius, and the fourth rank tensor β ijkl is The α's depend on the shear crack compliance and the β's depend on the shear and normal crack compliance in such a way that when ν 0 /2 is much smaller than unity, the contribution of β ijkl is small compared to that of α ij and the change in compliance due to cracks can be described by α 11 and α 33 only (Sayers & Kachanov, 1995).
By using Sayers and Kachanov (1995)'s expressions for the effective compliance S ij in Equations 6 and 7; Wong (2017) arrives at the following expressions for the Skempton coefficients for cracked materials: and where the approximation arises from neglecting terms of the order of ν 0 /2 compared to unity.
Here, our goal is to estimate the crack density tensor components α 11 and α 33 by using elastic wave velocity measurements on the dry rock, and comparing the prediction of Equations 14 and 15 to direct measurements under undrained conditions.

Sample Material and Preparation
In this experimental study, we used a cylindrical sample of Westerly granite.The sample, 40 mm in diameter and 100 mm in length, was cored from a large block, with the ends ground flat and parallel to ensure parallelism within a precision of 0.02 mm.To generate open microcracks in the material, the sample was heat treated at room pressure to 600°C (past the α/β phase transition of quartz).The heating rate was 3°C per minute, to minimize thermal gradients (Wang et al., 2013).The sample was left at 600°C for 3 hr, after which the temperature was decreased by 8°C per minute and then left overnight to cool to room temperature.After thermal treatment, the sample expanded from an initial length of 100.00-100.90mm and from a diameter of 40.08-40.40mm, representing a volume change of 2.5%.The porosity of the sample following thermal treatment was 3.2%, measured using the triple weighing method.Thin sections highlight the microcracks and grain boundaries in the thermally cracked sample (Figure 1).

Triaxial Apparatus
The experiment was performed in the large volume triaxial deformation apparatus located in the Rock and Ice Physics Laboratory at University College London (Eccles et al., 2005).Confining pressure is generated by compressing the silicone oil surrounding the jacketed sample using an electric pump, to an accuracy of 0.4 MPa, and measured at the inlet to the pressure vessel with a pressure transducer at 0.01 MPa precision.Axial load is applied by a piston, driven by a servo-hydraulic ram.Load is measured by an external load cell, and corrected for internal friction along the piston seals.Axial displacement is measured by an external Linear Variable Differential Transformer (LVDT), which is corrected for elastic distortion of the loading column (using a calibrated machine stiffness of 480 kNmm −1 ) in order to obtain sample shortening.
An independent pore fluid pressure is imposed at either the top and/or the bottom of the sample via a servo-controlled intensifier with a volume of 44 cm 3 .Steel end-caps accommodate pore fluid ports and disks distribute the fluid across the surface of each end.The pressure is monitored at both ends by pressure transducers.
An LVDT measures the position of the piston in the intensifier allowing for the change in volume of the pore system to be measured.This is used to calculate the change in pore volume of the sample.

Sample Instrumentation
The bottom plug of the pressure vessel, on which the test sample is located, is equipped with high-pressure lead-throughs accommodating up to 56 individual connections to be made between individual sensors and external recording equipment.Here, we used ultrasonic transducers (both P-and Sh-wave) to measure wave velocities (and estimate the sample dynamic compliance S ij ), strain gauges to calculate the static elastic moduli, and differential fluid pressure transducers (Brantut, 2020;Brantut & Aben, 2021) to calculate the pore pressure change from a change in stress (the Skempton coefficients, B x and B z ).
The sample was jacketed in a perforated nitrile sleeve (Figure 2a) and an array of 14 ultrasonic transducers (eight of them sensitive to P-waves and six to Sh-waves) was placed and sealed around the sample in a formation allowing for five P-wave and three Sh-wave ray-paths (Figures 2e and 2f) between signal and receiver transducer pairs.At repeated time intervals throughout the test, surveys were conducted where a 250 V high frequency pulse was sent to each transducer, and transmitted signals on the remaining sensors were amplified to 40 dB and recorded at 50 MHz with digital oscilloscopes.Two pairs of axial and radial 350 Ω strain gauges were positioned around the center of the sample 90° apart.To ensure a smooth surface for bonding, we applied a base layer of epoxy and then polished the area where the strain gauges were glued.The strain gauge signals were amplified, leading to an accuracy of the order of 10 −6 .
Two pore pressure transducers were positioned around the center of the sample 90° apart (Figure 2d), opposite the strain gauges, in direct contact with the surface of the sample.The transducers comprise a stainless steel stem with a 0.2 mm hole connecting the sample surface to a penny-shaped cavity of 3.5 mm radius and 0.2 mm height.
The hole and penny-shaped cavity have a combined dead volume of 2.89 mm 3 that accesses the pore space of the sample but is sealed from the confining medium using an o-ring between the cap and stem, and epoxy around the stem and jacket (Figure 2c).A full description of the transducers can be found in Brantut and Aben (2021).
The pore pressure transducers were individually calibrated for each test by systematically varying pore pressure and confining pressure within a range used in the test, with time allowed from equilibration between each step (Brantut, 2020;Brantut & Aben, 2021).The output (V) in volts was fitted to a linear combination of confining pressure and pore pressure, namely where P c is the confining pressure, p is the pore pressure, and a, b, and c are best fitting coefficients calculated using the least-square method.The coefficients were then used to calculate the pore pressure from the voltage output and measured confining pressure during the test.The relative precision was estimated to be of the order of 0.01 MPa.
The volume of the pore pressure transducers (2.89 mm 3 ) has been designed to be minimized, but it still has the potential to dampen the transient pore pressure depending on the diffusivity of the sample (Brantut & Aben, 2021).
When the pore pressure changes in the sample, the pore pressure within the sensor has to equilibrate with the pressure in the sample by fluid flowing in to or out of the sensor.When the diffusivity of the sample is low, it can take a substantial amount of time for the fluid to flow between the sensor and the sample and to reach equilibrium.
Unfortunately, a loose connection on one of the pore pressure transducers led to it only working intermittently.
For ease of data processing only the transducer that was working throughout was used.
All the mechanical (load, shortening, confining pressure), strain gauge and fluid pressure data were acquired at a sampling rate of 1-5 Hz.

Protocol
We conducted the experiment at two effective pressures (P c − p), 30 and 100 MPa, using the following sequence.
The confining pressure of the dry sample was first increased to 30 MPa, and differential (axial) stress was increased to 60 MPa at a constant strain rate of 10 −5 s −1 and then decreased back to 0. Confining pressure was then increased to 100 MPa, followed by a differential stress cycle up to 200 MPa and down to 0 again (Figure 3a).Differential stress was taken to a maximum of only twice the confining pressure to ensure that the sample remained in the elastic regime, below C′ (Wang et al., 2013).This entire initial sequence was used to determine the dry rock properties.
Subsequently, the sample was saturated by imposing 10 MPa pore pressure at one end, and venting the other end.
After steady flow was established, we considered the sample to be fully saturated.The pore fluid system was closed and pore pressure was equilibrated to 10 MPa.Then, confining pressure was cycled in steps from 20 to 110 MPa (Figure 3b), and permeability was measured (using the constant flow rate method, venting one end of the sample, giving a pore pressure difference of 10 MPa between sample ends) at selected pressures.
At 40 MPa confining pressure, differential stress was increased to 60 MPa (Figure 3b) at a constant strain rate of 10 −5 s −1 and then decreased back to 0, replicating the stress path of the sample when dry at 30 MPa confining pressure.Differential stress was again increased to a maximum of 60 MPa and at selected points axial and radial stress were independently cycled.Rapid stress changes induced an undrained response in the sample, where pore pressure initially changed linearly with stress due to the Skempton effect until it equilibrated in the sample and dead volume (consisting of the pore fluid pipes and distribution disks), the sample was reconnected to the pore pressure intensifier so that pressure diffused to the drained conditions.Therefore, each change in stress was initially undrained and then became fully drained in a sequence similar to the one used in Hart and Wang (2001).
We then repeated this process at 110 MPa confining pressure and selected differential stress values up to a maximum of 200 MPa.

Measurements of Drained Parameters
We measured axial and radial strain during increases and decreases in axial stress at a constant strain rate whilst the sample was under four conditions: (a) dry at 30 MPa confining pressure, (b) dry at 100 MPa confining pressure, (c) saturated at 10 MPa pore pressure and 40 MPa confining pressure, and (d) saturated at 10 MPa pore pressure and 110 MPa confining pressure.Drained Young's modulus and Poisson's ratio were calculated for decreasing stress under the four stress conditions.A quadratic fit (least-squares) of both axial strain with axial stress and axial strain with radial strain was obtained at each stress condition.Separately, the derivatives of the quadratic fit were taken to give the Young's modulus and Poisson's ratio as a function of axial strain.in the saturated sample, confining pressure was cycled, then cycles of axial (black circles) and radial (magenta circles) stress were cycled at increasing differential stress, this stage took approximately 16 days including saturation.The typical magnitude of the rapid stress cycles was 5-8 MPa. Figure 7 shows an example of one rapid axial stress step with the resulting pore pressure response.

Measurements of Undrained Parameters
We measured the change in pore pressure and the changes in axial and radial strain under undrained conditions during imposed constant rate (ramp) changes in axial and radial stress.By assuming fully undrained conditions, the Skempton coefficients were calculated as the linear fit of the pore pressure change with the stress change, and undrained static elastic moduli were calculated from the linear fit of strain with stress change.The typical magnitude of the stress steps was 5-8 MPa.
At each stress level for each of the constant rate changes in axial stress, the final six data points of the ramp were combined together and a linear fit of around 120 points was performed to calculate each value of B z ,   u  , and   u  (see Figure 6).The last points of the ramp were used to ensure a linear increase in stress after piston seal friction was overcome.
In a similar way, for radial stress changes six data points were taken after the stress change was linear (only for increasing steps).However, here the first six points were used as piston seal friction was not an issue for initial changes in radial stress.A linear fit of around 60 points of radial stress and pore pressure was taken for calculated value of B x .

Measurements of Elastic Tensor From Ultrasonic Data
The eight P-wave and six Sh-wave ultrasonic transducers were arranged to give measurements of P-and Sh-wave travel times along five and three different propagation angles with respect to the compression axis, respectively.We first manually picked absolute arrival times for each pulser-receiver pair on a reference set of transmitted waveforms obtained under hydrostatic conditions at the beginning of the experiment.Differences in arrival times between this reference and subsequent surveys were obtained by cross-correlation of P-and Sh-wave trains (see Brantut (2015)).The arrival times were used to compute P-and Sh-wave velocities along the five and three propagation angles, respectively.Because the rock is anisotropic, these velocities are group velocities (i.e., off axis waves are quasi-P and quasi-S), and their relationship with the elastic moduli of the rock is less straightforward than that of phase velocities.
Here, we assumed vertical transverse isotropy, where the symmetry axis is taken as the compression axis.
Our observations consist of eight measured group velocities along eight group angles, and we invert those for the elements of the compliance matrix S ij and eight phase angles.In practice, we follow the same approach as described in Brantut and Petit (2022), and solve the inverse problem d obs = g(m) using a quasi-Newton method (Tarantola, 2005, p. 69), where the model parameters m are the three Thomsen parameters ϵ, δ, and γ (Thomsen, 1986), the vertical P-and Sh-wave velocities, and the eight phase angles.The forward model (function g) first calculates phase velocities from the input Thomsen parameters (Thomsen, 1986, his Equation 10) at given phase angles, and predicts group velocities and group angles (Thomsen, 1986, his Equations 13 and 14).
For the first survey at each confining pressure, a priori model parameters were assumed to be isotropic (Thomsen parameters initialized to zero) and vertical P-and Sh-wave velocities were taken equal to those measured along the horizontal.For subsequent surveys, we used the Thomsen parameters inverted from the previous survey as priors for the inversion, and used the group angles as priors for the phase angles.
The best-fit Thomsen parameters were used to calculate the stiffness matrix (C ijkl ), which was then inverted for the compliance matrix (S ijkl ).The compliances were then used directly to calculate the five crack density parameters (α 11 , α 33 , β 1111 , β 1133 , β 3333 ) by solving the linear system described in Sayers and Kachanov (1995, their Equa tions 16-21).We assumed E 0 = 89 GPa and ν 0 = 0.22 to match intact P-wave velocity of 6.2 km s −1 and intact S-wave velocity of 3.7 km s −1 (Brantut & Petit, 2022).

Wave Velocities in the Dry Material
At 100 MPa confining pressure and low differential stress, the sample has minimal anisotropy, with a maximum difference of 200 m/s between horizontal wave paths and the most vertical wave paths (28°; Figure 4).There was an increase in both P-and Sh-wave velocity at increased differential stress at all angles, but wave velocity increased more at angles closer to vertical.Between 10 and 200 MPa differential stress, sub-vertical Sh-wave velocity (28° angle) increased by 2.3% and increased by less than 1% in the horizontal direction.Therefore, at increased differential stress there was increased anisotropy in the wave velocities.

Characterization of Transport Properties Under Saturated Conditions
Permeability was measured using the constant flow rate method under hydrostatic stress conditions as confining pressure was decreased from the maximum of 110 MPa.The permeability at 110 MPa is around 1 × 10 −18 m 2 increasing to 2 × 10 −17 m 2 at 20 MPa confining pressure (Figure 5).

Measurement of Undrained Quantities
We used the linear undrained part of the time-series of the axial and radial strains and pore pressure changes during the rapid ramps in confining pressure and differential stress to extract directly undrained quantities.We calculated undrained poroelastic parameters by measuring pore pressure, axial strain, and radial strain during cycles of increasing and decreasing ramp changes in stress (Figure 6).As stress is increased, pore pressure increased due to the Skempton effect under undrained conditions, as the duration of the stress increase is much smaller than the time for the pore pressure to diffuse.Stress was then held constant and pore pressure equilibrated within both the pore space of the sample and the dead volume of the pore fluid system.We then re-connected the sample to the pore pressure intensifier to allow pore pressure to equilibrate to the controlled level, producing drained conditions (Figure 7).The whole collection of results is shown in Figure 8.

Axial and Radial Skempton Coefficients
The axial and radial Skempton coefficients were calculated in three ways: (a) Directly using the pore pressure measurements using the process outlined in Section 3.6 (Figures 8a and 8b, solid symbols); (b) calculated from the compliances using Equations 6 and 7 (Figures 8a and 8b, solid lines), and (c) calculated from the crack densities using Equations 14 and 15 (Figures 8a and 8b, dashed lines).Fluid compressibility was taken as 0.45 GPa −1 and porosity was estimated at 1.2% and 0.9% at 30 and 100 MPa effective pressure, respectively.The strain gauges were used to calculate volumetric strain.The change in sample volume was considered to be due to both the change in the porosity (from an initial 3.2%) and the compression of the solid matrix material.
Overall, B z decreases and B x increases with increasing differential stress at both 30 and 100 MPa effective pressure (Figures 8a and 8b).The same trend is for all methods of measurement even though the absolute values differ.
Directly measured B z values (Figures 8a and 8b, triangles) are higher than those calculated from S ij (solid lines) from ultrasonic wave speeds by about 0.1-0.2units.Directly measured B x (circles) are within 0.2 of those calculated using the S ij , with directly measured values lower at 30 MPa effective pressure and similar at 100 MPa effective pressure.The Skempton coefficients calculated using crack densities (dashed lines) are similar to those directly measured for B x at 30 MPa effective pressure and B z at 100 MPa effective pressure, but 0.3-0.4 units higher for B z at 30 MPa effective pressure and B x at 100 MPa effective pressure.

Young's Modulus
At both 30 and 100 MPa effective pressure, undrained and drained axial Young's moduli increase with increasing differential stress (Figures 8c  and 8d).The dynamic Young's modulus (solid line) is higher than the static Young's modulus under both dry (black dotted lines) and saturated (blue  dotted line) conditions, but the difference decreases with increased differential stress, from around 20 to 5 GPa at 30 MPa effective pressure, and from around 10 to 0 GPa at 100 MPa effective pressure.
There is only a small (less than 3 GPa) difference between the static Young's modulus calculated in the dry sample (Figures 8c and 8d; dotted black line) and the static Young's modulus calculated in the saturated sample at the same effective pressure (dotted blue line).

Poisson's Ratio
The drained and undrained axial Poisson's ratios increase with increasing differential stress at 30 MPa and remain around the same with increasing differential stress at 100 MPa effective pressure (Figures 8e and 8f).The static drained Poisson's ratio calculated from strain measurements in the dry sample (black dotted line) is similar to that calculated in the saturated sample (blue dotted line).The dynamic Poisson's ratio (solid black line) is similar to the static at 30 MPa effective pressure, and marginally lower at 100 MPa effective pressure.The undrained Poisson's ratio (solid triangles) is higher than the drained static Poisson's ratio by 0.10 units at 30 MPa effective pressure, but only slightly higher at 100 MPa effective pressure.

Assessment of Model Predictions
The model of Wong (2017) provides a relationship between microstructural quantities (crack density parameters α ij ) and macroscopic poroelastic Skempton coefficients B x and B z (Equations 14 and 15), which has the potential to simplify considerably the estimation of poroelastic parameters in damaged rocks.To validate their model, Wong (2017) calculated crack densities α 11 and α 33 from measured values of B x and B z , and then used the calculated α 11 and α 33 to determine Biot coefficients, which were compared with measured values.Here, we were able to obtain direct measurements of Skempton   (black) is the controlled stress change, the right hand axis (blue) is the pore pressure change which first increases with stress after piston friction has been overcome (inset box), then equilibrates within the pore space and pore system by diffusing from the sample ends into a closed system, and then following the re-connection of the pore fluid system to the intensifier, pore pressure equilibrated at the imposed fluid pressure.The axial Skempton coefficient was calculated from the linear fit between pore pressure change with axial stress change.
coefficients and independent estimates of crack density parameters (Figure 8), and we observe a reasonable general agreement between our data and the model predictions.
An alternative approach to validate the model using our data is to compare our measurements of crack densities (α's) with crack densities calculated from Skempton coefficients (Wong, 2017, their Equation 15).Estimates of α 11 using the Skempton coefficients are around four to five times higher than the estimates using wave velocities (Figure 10, blue vs. black dots).Similarly, estimates of α 33 using the Skempton coefficients are around four times higher than estimates from wave velocities (Figure 10 blue vs. black triangles).The crack densities estimated from the Skempton coefficients appear much larger than what might be expected and show that the discrepancy between the model predictions and the experimental data is larger when analyzed in terms of inferred crack density than directly in terms of Skempton coefficient.
One of the limitations of the model is that the underlying non-interactive effective medium scheme is only appropriate for microstructures dominated by thin microcracks at low concentrations.In Berea sandstone, the crack densities calculated from the Skempton coefficients by Wong (2017) are about one order of magnitude higher than those inverted from wave velocity measurements in the same rock by Sayers and Kachanov (1995), which is around twice the discrepancy seen in our data set on cracked Westerly granite.The high crack densities inferred by Wong (2017) from poroelastic measurements are most likely due to the high porosity of Berea sandstone (ϕ = 21%), which contributes to lowering effective elastic properties but is neglected in the model.The low absolute values of crack densities inferred by Sayers and Kachanov (1995) are unfortunately not directly comparable to those of Wong (2017) in the same rock: in their inversion, Sayers and Kachanov (1995) used matrix elastic properties (E 0 and ν 0 ) that already included the contribution of nonclosable (i.e., equant) pores, so that the resulting crack densities reflect additional variations of elastic moduli around a homogenized porous bulk, rather than absolute values.It is likely that both approaches could be reconciled by using an effective medium model that includes both the contributions of (possibly oriented) microcracks and equant porosity, such as that of Shafiro and Kachanov (1997).
In the derivation of Equations 14 and 15; Wong (2017) used the simplifying assumption that terms of the order of ν 0 /2 could be neglected compared to terms of the order of unity.This assumption simplifies greatly the approximation of the effect of cracks on elastic compliances, which can be captured by a single second order tensor with components α ij .For Westerly granite, ν 0 /2 ≈ 0.11, so errors of the order of 10%-20% would be expected when fitting our wave velocity measurements to crack density parameters, and, eventually, to predicted Skempton coefficients.The differences between our calculations using Equations 6 and 7 and predictions using Equations 14 and 15 are larger than expected (Figures 8a and 8b, solid and dashed lines).At 100 MPa effective pressure, the magnitude of β 1111 was relatively large at around four times the size of α 11 (Figure 9), giving it around 40% of the effect of α 11 in the change in compliance of S 11 .Therefore S 11 in particular was overestimated (by at least three times as much as S 12 , S 13 , and S 33 ) which was the cause of the overestimation of B x .Similarly, at 30 MPa effective pressure, the magnitude of β 3333 was also relatively large at around four times the magnitude of α 33 .However, the magnitude of β 1111 was similar to α 11 so S 11 was not effected as much as the other compliances.
Our estimates of the compliances depend on absolute values of wave velocities that are prone to errors due to the manual picking of the master survey.In addition, estimates of the crack density parameters also depend on the moduli of the solid skeleton (E 0 and ν 0 ) which are not accurately known.Thus, the inverted β ijkl are sensitive to small absolute errors in their estimation and the discrepancy when neglecting β ijkl 's can either be due to errors in calculating β ijkl 's or in the model itself.

Dynamic Versus Static Parameters
Here we calculated static elastic moduli and poroelastic parameters using direct measurements of stress, strain and pore pressure.We compared these to dynamic elastic moduli and poroelastic parameters calculated from the relaxed ultrasonic wave velocities.At both 30 and 100 MPa effective pressure the dynamic Young's modulus is larger than the static Young's modulus by around 10-20 GPa (Figures 8c and 8d, solid line and solid triangles, respectively).At 30 MPa effective pressure the dynamic drained Poisson's ratio is similar to the static drained Poisson's ratio.At 100 MPa effective pressure the dynamic ratio is lower by less than 0.05 units.The dynamic Young's modulus is expected to be higher than the static one in Westerly granite, particularly when thermally cracked, but the expected change in the Poisson's ratio is not known (Blake et al., 2019).
Measurements of static elastic moduli require changes in stress during which cracks progressively close and reopen; whilst the opening and closing of cracks is reversible, the distribution of cracks is not the same at the start and end of the stress step.This is in contrast to the dynamic measurements, using high-frequency but low amplitude waves, where the small changes in strain mean that the crack fabric can be interpreted as frozen.Therefore, during dynamic and static measurements the rock is not in exactly the same microstructural state.
It is difficult to quantify the effect that the differences between static and dynamic moduli have on the calculated values of B x and B z .Equations 6 and 7 depend non-linearly on the compliances.Relative changes in S 11 and S 33 dominate as S 12 and S 13 are smaller in magnitude by a factor of ν x and ν z , respectively.
At higher effective pressures the compliances approach that of the intact material and the numerator and denominator of Equations 6 and 7 become smaller.Therefore, B x and B z become more sensitive to small absolute variations in measurements of compliances.The differences between drained dynamic and static elastic moduli suggests that when we calculate B x and B z using Equations 6 and 7 using dynamic measurements, we would expect a difference from directly measured values.
A decrease in S 33 (so increase in E z ) results in a decrease in B z and an increase in B x .A decrease in S 11 results in an increase in B z and a decrease in B x .Changes in S 11 dominate B x (for B x < 1) and changes in S 33 dominate B z (when B z < 1).Dynamic measurements of E z are greater than static measurements, therefore leading to a smaller value of S 33 .We do not know the difference between dynamic and static measurements of S 11 because we cannot statically measure S 11 in a triaxial apparatus.If we assume that differences between dynamic and static values of S 11 are the same as those of S 33 , then this would mean that dynamic measurements underestimate B x and B z .This is what we see for B z at both 30 and 100 MPa effective pressure.However, B x at 30 MPa effective pressure is larger using dynamic moduli and for B x at 100 MPa the dynamic and static measurements are comparable, when we might expect the dynamic measurements to be smaller.
As B x approaches unity, changes in S 11 become less dominant relative to changes in S 33 .The directly measured B x at 30 MPa effective stress is around 0.9.Therefore, if the difference between dynamic and static S 11 (not known) is not as large as S 33 then we would see an increase in B x compared to the directly measured value.
Another approach to compare dynamic and static measurements is to estimate the compliances using the principal crack densities estimated from the Skempton coefficients (see Figure 10, blue circles and triangles).These estimated compliances can then be compared to those determined from wave velocities and, in the case of the vertical direction (S 33 ), strain gauges (see Figure 11).If the theory and measurements directly aligned, we would expect that the compliance terms calculated from crack densities determined from wave speeds (Figure 11, black dashed lines) and calculated from the Skempton coefficients (Figure 11, blue triangles) to be the same, with a small difference from compliances directly determined from wave velocities due to neglecting β ijkl 's (Figure 11, black lines).However, compliances calculated using the Skempton coefficients are much larger.In particular at 30 MPa effective pressure, S 11 calculated from crack densities estimated from the Skempton coefficients (Figure 11a, blue triangles) are around four times as large as S 11 calculated from crack densities determined from wave velocities (Figure 11a, black dashed line).This discrepancy in S 11 is around twice as large at 100 MPa effective pressure (Figure 11b).S 33 estimated using the Skempton coefficients is around twice as large compared with estimates using the wave velocities (Figures 11c and 11d).

Storage Coefficients
To calculate the Skempton coefficients from the compliances, we required an estimate for storage using The estimated storage capacity decreased with increased differential stress at both 30 and 100 MPa confining pressure estimated using compliances and cracks (Figure 12).However, the change in storage was relatively small, at around a 10% decrease.Storage capacity estimated at 30 MPa confining pressure was around twice that at 100 MPa.The estimated storage capacities was around 1 × 10 −11 Pa −1 higher when estimated using only the crack densities (α ij 's) in place of the compliances (Figures 12a and 12b, dashed lines).The crack densities inverted from the directly measured Skempton coefficients were also used to estimate the storage capacity (Figure 10, blue dots).Estimates of storage were around two to three times higher using the crack densities from Skempton coefficients (Figure 12, blue dots) compared with estimates inverted from wave velocities (Figure 12, solid and dashed lines).
Figure 11.Components of the compliance matrix S 11 and S 33 calculated directly from wave velocities (solid black lines), from wave velocities but ignoring β ijkl 's (dashed black lines), from principal crack densities calculated from the Skempton coefficients (blue triangles), and, for S 33 calculated from strain gauges, the inverse of the Young's modulus shown in Figure 8 (dotted blue line).
Undrained and drained (axial) Young's modulus are related via the axial Skempton coefficient and the storage coefficient Therefore, it is possible to estimate the storage capacity using the measured E u , E d , and B z (Figure 8).The limitations of this approach are clear by observing the difference between the drained and undrained Young's modulus at is not detectable at 30 MPa effective pressure (Figure 8c); implying a storage capacity close to zero.However, the qualitative relationship of a decreasing difference between undrained and drained Young's modulus with increasing differential stress is consistent with the observed decreasing trend in the axial Skempton coefficient (B z ) (Figures 8a and 8b) and estimated storage capacity (Figure 12) with increased differential stress.

Effective Pressure
Ultrasonic wave velocity measurements were made on a dry sample at 30 and 100 MPa confining pressure to compare with drained measurements at 40 and 110 MPa confining pressure and 10 MPa pore pressure.We assume that properties measured under dry and saturated conditions can be compared because they are at the same effective pressure: the effective pressure coefficient is equal to unity.
There are different effective pressure coefficients for different processes (Nur & Byerlee, 1971).Here, we are comparing elastic wave velocities, porosity, and the crack orientations and distributions at different confining pressure and pore pressure combinations.Both porosity and crack closure are governed by an effective pressure coefficient which is likely close to unity (Zimmerman et al., 1986).

Micromechanical Assumptions
To estimate the Skempton coefficients using Equations 6 and 7 we are following the micromechanical analysis of Cheng (1997), which explicitly models the solid, fluid and pores of the material, and makes the assumptions of micro-isotropy and micro-homogeneity.The assumptions of micro-isotropy and micro-homogeneity are not always valid under isotropic conditions (Cheng, 2021;Hart & Wang, 2010;Pimienta et al., 2017).Without these assumptions, for transverse isotropy, we would need to measure the five independent unjacketed elastic moduli, and the directional variation in pore deformation.Neither of these are easily measurable in an anisotropic material.

Anisotropic Poroelastic Response
There have been few published results of directly measured anisotropic Skempton pore pressure coefficients.Aoki et al. (1993) measured B x = 0.51 (parallel to bedding plane) and B z = 0.53 (perpendicular to bedding plane)  17using compliances (solid lines) and crack densities (dashed lines) determined from wave velocities.These estimates were used to estimate the axial and radial Skempton coefficients.Additionally, the crack densities inverted from the Skempton coefficients were used to estimate storage capacity (blue dots).
from shale cored in two directions.Makhnenko and Tarokh (2018) report anisotropy of Skempton coefficients in Berea sandstone but at only one stress state.
The studies which present results for the anisotropic Skempton coefficients as a function of differential stress are those of Lockner and Stanchits (2002) and Lockner and Beeler (2003), and Wong (2017) includes some previously unpublished data from Lockner and Beeler (2003) when validating their model.
We compare the qualitative behavior of B x and B z with increasing differential stress at different effective pressures in both Berea sandstone and thermally cracked Westerly granite (Figure 13).For both rock types, at all confining pressures, B x increases with increased differential stress, and B z decreases.At lower confining pressure, the anisotropy in the Skempton coefficients is higher with both B x higher and B z lower, and the anisotropy increases more rapidly with increased differential stress than at higher confining pressure.The B x are higher in Westerly granite compared with Berea sandstone, even at much higher effective pressures.
At increased differential stress, horizontally aligned cracks close and the vertical compliance decreases whilst the horizontal compliance does not change much.This is the case regardless of the rock type or porosity.

Conclusions
Transversely isotropic drained and undrained elastic moduli and poroelastic parameters have been measured in Westerly granite under increasing differential stress.Direct measurements of pore pressure on the sample surface gave us independent measurements of the two Skempton coefficients.The five independent (dynamic) elastic moduli were inverted from elastic wave velocity data, and used to estimate crack density parameters.We then checked the consistency between undrained properties measured statically and those predicted from crack density parameters in the model of Wong (2017).
The axial Skempton coefficient (B z ) decreases with increased axial stress and the radial coefficient (B x ) increases, with the vertical transverse isotropy induced predominantly through closure of horizontally aligned cracks.The anisotropy in the Skempton coefficients increases with increasing differential stress.Our observations in Westerly granite at two effective pressures, follow the same qualitative trend as those of previous studies in Berea sandstone at multiple effective pressures (Lockner & Beeler, 2003;Wong, 2017).By making assumptions of micro-homogeneity and micro-isotropy in a micromechanical model (Equations 6 and 7) we can provide a reasonable consistency with directly measured Skempton coefficients, with differences between dynamic and static elastic moduli explaining the discrepancies.Furthermore, this qualitative behavior can be predicted by using only two principal components of the second-order crack density tensor (α 11 and α 33 ).

Figure 1 .
Figure 1.Combined reflected and transmitted light images of (left) intact and (right) thermally cracked Westerly granite prior to deformation.This combination was used to highlight microcracks and grain boundaries.

Figure 2 .
Figure 2. (a) A jacketed sample instrumented with ultrasonic and pore pressure transducers, (b) reference frame for the transducer positions, (c) schematic of a pore pressure transducer, (d) locations of pore pressure transducers, (e) locations of P-wave ultrasonic transducers, and (f) locations of Sh-wave ultrasonic transducers.

Figure 3 .
Figure3.A schematic of the stages of the experimental protocol: (a) in the dry sample confining pressure (magenta) was cycled, then differential stress (black) was cycled at two confining pressures, this stage took approximately three days; (b) in the saturated sample, confining pressure was cycled, then cycles of axial (black circles) and radial (magenta circles) stress were cycled at increasing differential stress, this stage took approximately 16 days including saturation.The typical magnitude of the rapid stress cycles was 5-8 MPa.Figure7shows an example of one rapid axial stress step with the resulting pore pressure response.

Figure 4 .
Figure 4. Modeled versus measured group velocities while decreasing differential stress at 100 MPa confining pressure.Using this we verify that the velocities required for the inversion (solid magenta points) are representative of the measured velocities (open black points).Note that the scales for the P-wave velocity and Sh-wave velocity are not the same.

Figure 5 .
Figure5.Permeability (log scale) measured under hydrostatic conditions at selected levels of confining pressure.

Figure 6 .
Figure 6.Example strain and pore pressure response during a step in differential stress at 100 MPa effective confining pressure and 25 MPa differential stress.(a) Pore pressure as a function of differential stress; (b) pore pressure as a function of radial stress; (c) axial strain as a function of axial stress (undrained Young's modulus); (d) radial strain as a function of axial strain (undrained Poisson's ratio).

Figure 7 .
Figure 7.Example of pore pressure response to a stress step.Left hand axis (black) is the controlled stress change, the right hand axis (blue) is the pore pressure change which first increases with stress after piston friction has been overcome (inset box), then equilibrates within the pore space and pore system by diffusing from the sample ends into a closed system, and then following the re-connection of the pore fluid system to the intensifier, pore pressure equilibrated at the imposed fluid pressure.The axial Skempton coefficient was calculated from the linear fit between pore pressure change with axial stress change.

Figure 8 .
Figure 8.At 30 and 100 MPa effective pressure: (a) and (b) show the Skempton coefficients B x and B z calculated in three ways: the triangles (B z ) and circles (B x ) are calculated from direct measurements of pore pressure.The solid lines are calculated from the moduli inverted from the wave velocities and dashed lines are calculated from crack densities.(c) and (d) shows the Young's modulus (E z ): undrained (solid triangles), drained static (dotted line from stress-strain curve: blue is saturated; black is dry), drained dynamic (solid line from wave velocities, dry conditions).(e) and (f) are similar to (c) and (d) but for Poisson's ratio.

Figure 9 .
Figure 9. Crack densities estimated from the compliances determined from wave velocities at 30 MPa confining pressure (left) and 100 MPa (right) assuming dry penny-shaped cracks.

Figure 10 .
Figure 10.Principal components of the second rank crack density estimated from the compliances determined from wave velocities (black) compared with those estimated from directly measured Skempton coefficients (blue) Wong (2017) at 30 MPa confining pressure (left) and 100 MPa (right).

Figure 12 .
Figure 12.Estimates of storage capacity at 30 MPa confining pressure (a) and 100 MPa confining pressure (b) using Equation17using compliances (solid lines) and crack densities (dashed lines) determined from wave velocities.These estimates were used to estimate the axial and radial Skempton coefficients.Additionally, the crack densities inverted from the Skempton coefficients were used to estimate storage capacity (blue dots).

Figure 13 .
Figure13.Comparison of the axial Skempton coefficient (B z , triangles) and the radial Skempton coefficient (B x , circles) for (a) Berea sandstone (solid symbols) from Wong (2017) (using data fromLockner and Beeler (2003)) and (b) Westerly granite from this experiment.Matching colors on each subfigure are at the same confining pressure which is given by the annotated value.