Exploring Trends in Microcrack Properties of Sedimentary Rocks: An Audit of Dry and Water Saturated Sandstone Core Velocity–Stress Measurements ()
1. Introduction
Non-linear or stress dependent rock physics models are being applied increasingly to model the influence of stress perturbations due to reservoir production and injection activities on seismic velocities. Laboratory measurement of non-linear rock physical properties of dry core samples can provide valuable information on the stress dependent elastic properties of reservoir rocks [1- 3] and have the potential for up scaling to seismic frequencies [4] as well as relating to static elasticity [5]. More importantly, core measurements can also be used to calibrate rock physics models [6] for the forward prediction of the stress dependence of seismic velocities. Reservoir rocks are seldom under dry conditions and so more realistic characterization of the non-linear rock physical properties should be examined under fluid saturated conditions. However, the number of published studies on saturated core samples are fewer than those on dry samples and this is primarily due to the difficulty of carrying out fluid saturated stress measurements of core samples (e.g., enormous equilibration times are necessary when performing saturated core measurements in comparison with dry sample measurements).
In this paper, we compare the microcrack parameters of the discrete and analytic microstructural stress-dependent model described in [3] for dry and water saturated core. The data used in this study come from the sandstone ultrasonic velocity-stress measurements of [7] and so allow a direct comparison between dry and water saturated microcrack parameters. This work follows from [6] who explore the microcrack properties of over 150 dry-core ultrasonic velocity-stress measurements. Given the influence of fluids within microcracks (and capillary forces) has important implications on the stress sensitivity of reservoir rocks [8], it is necessary to study the effect of fluid saturation on non-linear rock physics model parameters. Thus, the main objective is the paper is to study the influence of fluid saturation fluid on microcrack properties to further calibrate the microstructural non-linear rock physics models discussed in [6].
2. Rock Physics Model
We examine the microcrack properties of two non-linear rock physics models: A discrete microcrack model defined by a secondand a fourth-rank crack density tensor [1-3] and an analytic microcrack model defined by an initial crack density and initial aspect ratio [3,9]. Although the discrete model describes the non-linear dependence of velocity with stress, the input parameters are two tensor quantities that are not necessarily intuitive. The analytic formulation provides a model based on physically intuitive input parameters to forward model the non-linear stress dependence of velocity, yet requires assuming that the microcracks are penny-shaped (i.e., the scalar crack approximation). Although penny-shaped cracks offer intuitive parameterization of the pore space and a reduction in the model complexity, the scalar crack approximation is not totally consistent with ultrasonic core data [6,10]. However, it should be noted that the analytic formulation still captures some of the essential stress dependent behaviour of sedimentary rocks and has utility for forward modelling applications.
It should be noted that the non-linear formulation we examine in this paper is one of many approaches to model the influence of stress on seismic velocity. For instance, [11] present a 1D empirical formulation to describe vertical travel time perturbation due to changes in vertical strain and vertical velocity from 4D seismic data. Authors [12-14] use third-order elasticity theory to characterize 3D stress dependence elasticity and anisotropy. Authors [15-17] introduce nonlinear models consistent with empirically derived phenomenological equations [18]. Our interest in the discrete and analytic models described earlier (and in more detail below) is based on seeking formulations described using few and intuitive effective microstructural model parameters that can be calibrated with available data (e.g., ultrasonic core data).
2.1. Discrete Microcrack Model
Reference [19] adopt the excess compliance approach of [20] to model the influence of stress dependent elasticity due to the deformation of microcracks. The stress dependence and elastic anisotropy is given in terms of an excess compliance ΔS (the inverse of the 3 × 3 × 3 × 3 elasticity tensor C)
(1)
[19,20] where δij is the Kronecker delta and summation convention is being used. The secondand fourth-rank crack density tensors αij and βijkl are defined
(2)
and
(3)
where V is volume and n is the unit normal to the displacement discontinuity set m (i.e., microcrack or grain boundary). 
and 
are the normal and tangential compliances across the microcrack set m having surface area Sm. The effective compliance S of a rock can be expressed
(4)
where S0 is the background (or intact) rock compliance estimated from either mineral composition [21] or high confining stress behaviour [22]. In this paper, we use the high stress approach because we have found using mineral composition does not yield consistent and reliable velocity predictions compare with observation.
The key assumptions for this model are that the microcracks are rotationally invariant and thin. Since ultrasonic measurements are only performed in one direction, we assume the samples are isotropic [22]. Thus, the fourth-rank crack density term can be simplified (i.e., βijkl is isotropic)
(5)
and
(6)
where
and 
(7)
[22]. The scalar N is the number of discontinuities in V, and r is the radius of the crack.
2.2. Analytic Microcrack Model
To enable forward modelling of 4D seismic effects related to perturbations in stresses [23-25], [3] extended the analytic effective medium formulation of [9] to predict ultrasonic anisotropic and stress-dependent velocities. Specifically, the analytic microcrack model introduces initial microcrack aspect ratio and number crack density to predict stress dependence and crack-induced elastic anisotropy. The number crack density is written
(8)
where
(9)
and 
are the effective initial number crack density and effective initial aspect ratio, λi and μi are the Lame constants, and 
is the principal effective stress in the i-th direction. The second-rank microcrack density term is
where 
(10 )
is a normalization factor [26], and 
and 
are the anisotropic intact rock Young’s modulus and Poisson ratio. This derivation yields an expression for the effective elasticity that can model stress-induced elastic anisotropy due to deviatoric stress fields. The key assumptions for this model are that the microcracks are pennyshaped and that the rock does not undergo brittle or plastic deformation.
2.3. Microcrack Properties of Dry Core
Reference [6] compiled a database of over 150 dry-core velocity-stress measurements to explore microcrack properties of the discrete and analytic non-linear rock physic models. Their results indicate that for most lithologies the initial aspect ratio are approximately 0.0005, but can be larger for shales. The initial crack density is sensitive to core damage and consolidation. Most notably, [6] note that the global trend of the compliance ratio 
is not necessarily unity and, for the samples analyzed, is approximately 0.6. This has important implications because, for most sedimentary rocks, the fourth-rank term βijkl is often neglected to enable characterization of the nonlinear stress dependent elasticity based solely from contribution of the second-rank term αij [2,3,27]. Deviations from the scalar crack assumption (i.e., where βijkl is small such that 
= 1) potentially result from several factors, namely presence of fluids with nonzero bulk modulus, clay within cracks, cementation, and complex crack geometries. In this paper, we know the saturating fluid as well as an estimate of the clay content of the rock specimens (see discussion below), but lack quantitative measures of cementation or microcrack geometry.
In [3,6], the analytic model parameter inversion involved implementing a simple grid search over model parameters to minimize the misfit between model predictions and observed ultrasonic data. However, such a simple grid search may not be an efficient method for determining the best-fitting model parameters. For the discrete model parameter inversion in [3,6], a Newton-Ralphson approach was used to minimize the misfit between the model predictions and observations based on derivatives of the elasticity tensor components with respect to model parameters (see Figure 1 in [6], for general workflow). A known limitation of the Newton-Ralphson method is that for nonlinear inversion problems it is often difficult to
Figure 1. Microcrack properties for sandstone sample 10381: Dry (top) and saturated (bottom) sample. In this figure and Figures 2 and 3, the right panel compares the ultrasonic data (Vp is open black circle and Vs is open gray triangle), the crack density inversion results (Vp is solid black circle and Vs is solid gray circle), and the analytical microcrack prediction based on the best fitting initial crack density and initial aspect ratio (Vp is black solid curve and Vs is gray solid curve). The top-left panel shows the BN/BT ratio of the inverted crack densities (open circles) and the mean BN/BT ratio (solid line). The bottom-left panel compares the best fitting crack densities from the crack density inversion (open circles) and the best fitting crack densities from the analytic microcrack prediction (solid curve).
find the true global minimum and so solutions may be biased towards local minima if the initial starting model is not chosen carefully. In this paper, we use the neighbourhood algorithm of [28] to improve the model parameter inversions for both the analytic and discrete formulations. The results of the inversion for the discrete model are exceptional (e.g., see Figures 1-3, where solid symbols represent discrete model predictions and open symbols the ultrasonic data). For the most part, the inversion results for the analytic model are poorer when compared to the discrete model. However, the analytic model predicts the general trend of the stress dependence remarkably well considering it only considers the influence of second-rank crack density tensor (e.g., the solid curves in Figures 1-3).
3. Data
Reference [7] investigated 69 sandstone core samples to examine the influence of fluid saturation on ultrasonic velocities. The measured porosity ranged between 5% and 30%, and the measured clay content ranged between 0% and 50%. The core samples ranged in length between 2 cm and 5 cm, and had diameter of 5.0 cm. It was noted that the dimensions of the core were approximately two orders of magnitude larger than the average grain size.
Figure 2. Microcrack properties for sandstone sample 12677: Dry (top) and saturated (bottom) sample.
Figure 3. Microcrack properties for sandstone sample Fountian B: Dry (top) and saturated (bottom) sample.
The Pand S-wave velocities were measured with the pulse transmission technique based on picking the first arrival peak amplitude. The confining and pore pressure were controlled separately, where the differential pressure limited to 50 MPa and the pore pressure to 1 MPa. The ultrasonic frequency of the Pand S-wave transducers were 1.0 MHz and 0.6 MHz, respectively. Based on the measured velocities, the average dominant wavelengths of the P and S waves were at least five times the mean grain size of the samples. Under dry conditions, the P-wave arrival times were picked to within 0.003 μsec (which equates to ≤1% error in velocity). The S-wave velocity errors were estimated to be less than 2%, except for the poorly consolidated samples at low confining stresses where the velocity errors were up to 3%. The samples were preloaded to 50 MPa and the velocities were measured on the unloading path to reduce the effects of hysteresis. Although hysteresis was observed the magnitude was small (i.e., ≤1%). For water-saturated conditions, the samples were fully saturated with water. However, for the samples having high clay content, the saturating fluid used was brine to minimize chemical alteration effects. [No velocity differences were observed between the water and brine saturation samples.] Velocities were measure during loading and unloading with only minor hysteresis being observed (≤1% for well consolidated samples and ≤2% for poorly consolidated samples).
4. Results
Figures 1-3 show the results of the inversion for three of the sandstone samples. Figure 1 is a sandstone sample with high clay content (46.0%) and low porosity (13.8%). The fit between observation and the discrete and analytic predictions are very good. The dry measurements have a characteristically higher stress dependence compared to the saturated measurements. The 
ratio for the dry measurements show a stress dependence, decreasing from 1.5 to 0.5 with increasing stress, whereas the ratio for saturated measurements is approximately constant at 0.5. Estimates of initial crack density e0 is smaller whereas the initial aspect ratio a0 is greater for the saturated measurement compared to the dry measurement. Figure 2 is a sandstone sample with low clay content 7.0% and high porosity 27.05%. This sample shows similar stress dependent velocity characteristics to the sample shown in Figure 1. However, the 
ratio for the dry measurements is approximately constant around 0.75 and there is minimal change in the predicted initial crack density and initial aspect ratio between the dry and saturated measurements. Figure 3 is a sandstone sample with no clay (0.0%) and moderate porosity (19.8%). This sample displays the same velocity stress dependence (i.e., higher stress dependence for the dry measurement). The 
ratio is also stress dependent, but displaying a concave upward trend with increasing stress compared to the concave downward trend in Figure 1. Although the initial crack density is lower for the saturated sample, there is no change in the initial aspect ratio. [Note, the model parameters are assumed to be isotropic only because the data contain only one Pand one S-wave measurement for each dry and saturated sample. However, the formulation does consider anisotropy in the model parameters if there is sufficient ultrasonic data.]
4.1. Analytic Model Parameters
Figure 4 displays e0 versus a0 estimates for the analytic model for both the dry and saturated measurements. Also shown (inset) are histograms for e0 and a0. For the dry and saturated measurements, the initial aspect ratios show similar clustering centred around 0.0005 and is consistent with that observed by [6] and references within). The initial crack density e0 show more scatter, with values between 0.0 and 0.4 for the dry measurements and 0.0 and 0.25 for the saturated measurements. In Figure 5, initial crack density is plotted with respect to clay content and porosity, and shows no clear relationship