Stress path dependence of time-lapse seismic effects in shales: experimental comparison at seismic and ultrasonic frequencies

SUMMARY

Changes in pore pressure within geological reservoirs, due to, for example hydrocarbon production, CO₂ and energy storage or wastewater disposal, may cause substantial stress changes in the overburden, altering propagation velocities of elastic waves. The corresponding time-shifts are detected and quantified using time-lapse (4-D) seismic analysis. To invert seismic time-shifts for changes in stress and strains, stress sensitivity of rocks is studied in laboratory experiments on core plugs. Such measurements are typically conducted at ultrasonic frequencies. However, previous studies indicate that the stress sensitivity of velocities at seismic frequencies could be higher than that at ultrasonic frequencies. Therefore, calibration based on laboratory ultrasonic data may lead to inaccurate prediction of stresses and strains when applied to 4-D seismic data. To study the influence of frequency on stress sensitivity of acoustic wave velocities, a series of laboratory experiments was performed on two overburden shales with different petrophysical properties. In a low-frequency apparatus—a triaxial pressure cell that combines measurements at low (seismic) and high (ultrasonic) frequencies—the shale samples underwent stress changes with different ratios of horizontal to vertical stress amplitudes to mimic stress variations across the overburden. High-frequency velocity changes were directly recorded, while low-frequency velocity changes were obtained indirectly from the elastic parameters measured at seismic frequency by applying a rock physics inversion using third-order elasticity model. The experiments were conducted at undrained conditions, a representative state for reservoir overburden composed of shales. The results suggest that the stress sensitivities and strain sensitivities (R-factor) of P-wave velocities could be two to four times greater at seismic frequencies than at ultrasonic frequencies. Furthermore, it was found that the previously reported linear relation between the stress sensitivity and stress-path parameter (horizontal/vertical stress change) at ultrasonic frequencies also holds for seismic frequencies. We discuss the theoretical background for frequency-dependent stress sensitivity of wave velocities that supports the experimental findings. The effect of frequency when using laboratory-based calibrations should be taken into account when inverting time-lapse seismic data for changes in stresses and strains.

1 INTRODUCTION

Reservoir depletion by hydrocarbon production or inflation by injection of CO₂, water and other fluids results in stress changes in and around the reservoir. Stress changes and related strains may cause elastic properties of rocks to change, leading to changes in wave velocities, which are reflected by changes in the two-way traveltimes of seismic waves in time-lapse (4-D) seismic surveys. Due to the large thickness of the overburden (typically several kilometres), stress and strain-induced velocity changes can accumulate to significant and detectable time-shifts in the seismic signal (Barkved & Kristiansen 2005; Hatchell & Bourne 2005; Røste et al. 2006, 2007; MacBeth et al. 2018). Therefore, understanding the stress and strain sensitivities of the elastic properties of overburden rocks, in particular shales, is essential for the interpretation of 4-D seismic data used for detection of undepleted pockets of hydrocarbons, monitoring of reservoir fluid movement and caprock integrity (e.g. during waterflooding or CO₂ storage). Shale rocks are emphasized due to their abundance as overburden material and therefore their significant role in controlling the mechanical behaviour and seismic response of the subsurface. Beyond applications in hydrocarbon recovery and CO₂ storage monitoring, our findings are also relevant for assessing earthquake source regions, monitoring and prediction of landslide, assessing actively deforming geological provinces and geothermal reservoir management (Malehmir et al. 2016; Wagner & Uhlemann 2021).

Inverting changes in seismic signals for geomechanical changes in shales requires an understanding of how their seismic velocities depend on stress and strain. Multiple laboratory studies focusing on wave velocity stress sensitivities in shales have been published over the last few decades (Jones & Wang 1981; Johnston 1987; Johnston & Christensen 1995; Vernik & Liu 1997; Hornby 1998; Wang 2002; Dewhurst & Siggins 2006; Sarout et al. 2007; Pervukhina et al. 2008; Dewhurst et al. 2011; Szewczyk et al. 2017; Holt et al. 2018). However, only a few of the published tests were done under undrained conditions. For time-lapse seismic analysis, it is important to obtain undrained properties of low-permeability shale formations (typical permeability in the order of nanoDarcy, 10⁻²¹ m² (Horsrud et al. 1998). Utilizing classical consolidation theory (Biot 1941; Fjaer et al. 2021) and empirical data, Duda et al. (2023) provide estimates for diffusion times and argue that the timescales for achieving pore pressure equilibrium within the significant (for seismic monitoring) volumes in overburden shales extend well beyond the operational duration of typical reservoirs. Drained laboratory measurements are also of interest for understanding physical behaviour of shales, however, due to long consolidation time of field cores, drained stress-change measurements are challenging to complete within a reasonable time frame.

To quantify the effect of stress changes in different principal directions of stress, a stress-path parameter |$\kappa$| is defined as a ratio of minimum horizontal stress (⁠|${\sigma }_{\mathrm{h}}$|⁠) change to vertical stress (⁠|${\sigma }_{\mathrm{v}}$|⁠) change:

The majority of the published laboratory studies report only isotropic stress paths (⁠|$ \Delta \sigma _{\mathrm{h}} = \Delta \sigma _{\mathrm{v}}$|⁠, hence |$\kappa =1$|⁠), while a few consider triaxial stress paths (⁠|$ \Delta \sigma _{\mathrm{h}} = 0$| while |$ \Delta \sigma _{\mathrm{v}} \ne 0$|⁠, so |$\kappa =0$|⁠) without providing comparison to isotropic ones (Dewhurst et al. 2011; Sarout et al. 2007; Pervukhina et al. 2008). However, the stress path in the overburden deviates from pure isotropic or triaxial. For example, using the classical model by Geertsma (1973a, b) which assumes no elastic contrast between a disk shaped reservoir and its surroundings, one finds that the overburden is mostly following stress paths obeying a condition of constant mean stress (⁠|$ \Delta \sigma _{\mathrm{mean}} = 2\Delta \sigma _{\mathrm{h}} + \Delta \sigma _{\mathrm{v}} = 0$| as |$ \Delta \sigma _{\mathrm{v}} = -2\Delta \sigma _{\mathrm{h}}$|⁠, and therefore |$\kappa =-0.5$|⁠). The actual overburden stress path may however be quite different, as it depends on geometry and orientation of the reservoir, on contrast of elastic properties between overburden and reservoir, on elastic anisotropy and on non-elastic and non-linear mechanisms (Yan et al. 2023; Lozovyi et al. 2020). These factors can only be handled properly by numerical modelling, for example, using finite element method. The complex stress-strain fields in the overburden resulting from pore pressure changes in the reservoir are observed in field case studies (De Gennaro et al. 2008). Thus, there is a need to understand how wave velocities depend on stress changes for different values of stress-path parameter.

To our knowledge the only systematic experimental study of the stress-path dependence of ultrasonic velocities of different preserved field shales was reported by Holt et al. (2018). Yet, it was not clear if the observed dependencies would also hold for seismic frequencies. The laboratory velocity measurements are conventionally done at ultrasonic frequencies. To our knowledge, only Szewczyk et al. (2017) showed velocity stress dependence measured on shales at seismic frequencies. These authors point out that the stress sensitivity of Mancos shale at seismic frequencies might be significantly higher than at ultrasonic frequencies (by up to a factor of 2.5). Yet, only isotropic stress changes were applied, and all tests were performed at atmospheric pore pressure.

Several types of stress sensitivity models are available in the literature. Grain contact models, effective medium models, or phenomenological models can be used to predict velocity changes with stress. For shales, conventional grain-contact models like the Hertz–Mindlin (Hertz 1882; Mindlin 1949; Mindlin & Deresiewicz 1953) model are not applicable since the shape of clay particle is far from spherical. Effective medium models consider softening of the rock due to crack- and pore-like inclusions with different orientations and grain boundaries (Sayers 1999; Shapiro & Kaselow 2005; Fjaer 2006; Sayers 2006; Pervukhina et al. 2008; Ciz & Shapiro 2009). Phenomenological models, account for non-linear elastic effects by introducing higher-order elastic constants fulfilling given symmetry conditions. A simple phenomenological model involving only three non-linear coefficients was presented by Prioul et al. (2004). To overcome Prioul’s model limitation of isotropic strain sensitivity of the velocities, Duda et al. (2020) introduced a non-linearly elastic model based on earlier work by Fuck & Tsvankin (2009), which uses eight third-order elastic coefficients to fully characterized strain sensitivity of a TI medium. Bakk et al. (2024) further extended this model to account for the impact of shear strains on seismic velocities. Holt et al. (2018) suggested an empirically based model where velocity changes are linked to principal stress and pore pressure changes and the stress sensitivity is linearized with respect to stress-path parameter |$\kappa$|⁠.

We report the results of laboratory experiments with two fully saturated overburden shales with different clay content, porosity and depth. During the tests, samples were loaded to equivalent effective in situ stress and subjected to undrained loading-unloading cycles using stress-path parameters, |$\kappa$|⁠, in a range from −0.5 to 1. Elastic properties at both seismic and ultrasonic frequencies were quantified at each stage of the tests. Then, high-frequency velocity changes were directly computed from the ultrasonic data. Low-frequency velocity changes, due to limitations of the technique, were obtained indirectly from the elastic parameters measured at seismic frequency by applying a rock physics inversion using third-order elasticity model. As a result, the stress and strain sensitivities of axial P-wave velocity as a function of stress-path parameter (⁠|$\kappa$|⁠) were compared at seismic and ultrasonic frequencies.

The paper is organized as follows: first, we describe the sample materials and laboratory experimental methods; secondly, we detail the data analysis and data processing procedures, including uncertainty estimates for the experiments and data processing; thirdly, we present the experimental results and the results of seismic data processing; fourthly, we discuss the seismic data processing aspects, the relevance of the results for seismic inversion, and provide theoretical insights into the frequency dependence of wave velocities and their stress sensitivities; finally, we conclude the paper.

2 MATERIALS AND EXPERIMENTAL METHODS

2.1 Samples

Two Middle Miocene shales that represent overburden rocks covering hydrocarbon reservoirs at different geographical locations were used for the laboratory studies. After extraction from the well bore, rock cores were preserved by sealing them in an airtight and inert environment to protect the sample material from atmospheric conditions, moisture, contamination, and other factors that could alter its original state. The preservation process ensures the samples maintain their original saturation, preventing any changes in the rock’s physical and chemical properties. Best practices for shale core preservation and handling are discussed in Ewy (2015). In our laboratory, cylindrical core plugs with diameter of 25.4 mm and length of 48–52 mm were drilled and trimmed using Marcol oil as cooling agent and air barrier. After preparation, the shale plugs were stored in Marcol oil until the tests were started. The tested shales received code names T and M shale for consistency with previously published studies (Lozovyi & Bauer 2019a; Bakk et al. 2024). Petrophysical properties of the tested shales (Table 1) indicate difference between M and T shales in density, porosity and clay content which is attributed to different depth and geological stress conditions. Stresses applied in the low-frequency tests (Table 1) were selected to be equivalent to in situ through using the same net stresses at reduced pore pressures. Net stress, |$\sigma ^{\prime }$|⁠, is defined as the difference between total external stress, |$\sigma$| and pore pressure, |$p_\mathrm{f}$|⁠:

also known as Terzaghi’s effective stress. Under the tests, pore pressure was applied to the samples using a hydraulic system that was filled with equivalent to the in situ brines in terms of chemical activity to minimize osmotic ion exchange between specimens and the pore pressure fluid.

Both T and M shale core plug samples were tested in the low-frequency apparatus to obtain seismic and ultrasonic elastic properties. Additional ultrasonic data set from a test on M shale in a standard triaxial cell and under in situ stress condition was used to complement our data (Appendix  C). A comparison of the ultrasonic data of the low-frequency test (measured at in situ net stress) and complementary test (measured at in situ stress) did not reveal considerable differences, neither in absolute velocities, nor in stress sensitivities. Thus, the use of effective net stress is considered an adequate condition for studying stress path dependence of elastic properties.

2.2 Experimental setup

In this study, we used a low-frequency apparatus (Fig. 1), as described in detail by Szewczyk et al. (2016). Additional information as well as error analysis is found in Lozovyi & Bauer (2019b). This pressure cell was developed for measuring dynamic elastic moduli of samples at seismic frequencies and ultrasonic P- and S-wave velocities. In the cell, rock samples with dimensions of 25.4 mm in diameter and |$\sim$|50 mm in length can be subjected to various stress states by independently controlling axial and radial stresses, as well as pore pressure. Strains are measured in both axial and radial direction with four biaxial resistive strain gauges, and additionally by a set of three LVDT measuring strain in the axial direction (see Fig. 1, right-hand panel). The sample is sealed from the confining fluid by a rubber sleeve. The dead volume of the pore-fluid system under undrained conditions is equal to approximately 2 ml, which is several times smaller than typical sample’s pore volume, and considered sufficiently small for reliable measurement of pore pressure in undrained experiments.

2.3 Seismic and ultrasonic measurements

The principle of low-frequency (seismic) measurement is to apply small sinusoidal axial stress oscillations of a certain frequency, and record axial and radial strain modulations. A piezoelectric actuator (Fig. 1, left-hand panel) is controlled by a sinusoidal signal and creates force oscillations with a given frequency and amplitude. Frequency is limited to the range of 0.1–150 Hz. The resulting force modulation is measured by a piezoelectric force sensor, and the strain modulations are measured by strain gauges attached to the surface of the sample. Strain-modulation amplitudes are kept below |$10^{-6}$| to stay in the elastic regime—for higher strain amplitudes, the dynamic stiffness is impacted by inelastic effects (Winkler et al. 1979; Tutuncu et al. 1998; Lozovyi et al. 2017). Dynamic Young’s modulus, |$E$|⁠, and Poisson’s ratio, |$\nu$|⁠, are obtained from the stress- and strain-modulation amplitudes by the following relations:

where |${\varepsilon }_{\mathrm{ax}}$| is axial strain and |${\varepsilon }_{\mathrm{r}}$| is radial strain.

A conventional pulse-transmission technique (Krautkrämer & Krautkrämer 1977) is used for measuring ultrasonic P- and S-wave velocities in the axial direction. Ultrasonic P- and S-wave transducers and receivers with a central frequency of 500 kHz are mounted inside titanium endcaps of the axial piston.

2.4 Tests

For a cylindrical core sample, stress-path parameter, |$\kappa$|⁠, is defined as the ratio of radial stress change (⁠|${\sigma }_{\mathrm{r}}$|⁠) to axial stress change (⁠|${\sigma }_{\mathrm{ax}}$|⁠): |$\kappa =\Delta {{\sigma }_{\mathrm{r}}}/\text{ }\Delta {{\sigma }_{\mathrm{ax}}}$| (equivalent to eq. 1). In our experiments we operate with four types of stress path—undrained loading-unloading cycles defined by different stress-path parameter |$\kappa$|⁠. The stress paths are referred as follow:

• isotropic (⁠|$\kappa =1$|⁠)

• uniaxial strain (⁠|$\kappa =K_{0}$|⁠)

• triaxial (⁠|$\kappa =0$|⁠)

• constant mean stress (CMS) (⁠|$\kappa =-0.5$|⁠).

The |$K_{0}$| value is individually determined for each sample at uniaxial strain conditions (radial stresses are adjusted in real-time to provide no lateral deformation of the sample during axial loading, also known as oedometric conditions).

The low-frequency apparatus does not allow for a complete stiffness characterization of a transversely isotropic (TI) rock using a single core plug, neither at seismic, nor at ultrasonic frequencies. Thus, for each of the shales, three cylindrical samples drilled at angles of 0°, 45° and 90° relative to the bedding-plane normal (Fig. 2) were tested. A set of elastic properties measured from three core plug orientations reveal the five independent dynamic stiffness parameters (see Section 3 for a detailed explanation).

The 0|$^\circ$| samples were used for stress (and strain) sensitivity studies. Samples were first loaded to the respective in situ net stress state (see Table 1). After consolidation, undrained cycles with axial stress amplitude of |$\sim$|5 MPa were conducted while confining pressure was adjusted according to the selected stress path (an example of experimental timeline is shown in Fig. 3). Additional tests with 45° and 90° oriented samples were conducted for both shale types at initial stress states only.

Ultrasonic velocities were recorded continuously during the whole test, whereas the low-frequency measurements, due to technical limitations of the data acquisition system, were done at the following stages for each stress cycle: (i) before the loading cycle; (ii) at loaded state, after consolidatio and (iii) after unloading and consolidation. The low-frequency measurements are shown with the numbers 1–12 in Fig. 3.

3 DATA ANALYSIS

In the low-frequency apparatus, two different sets of dynamic properties are obtained from seismic and ultrasonic measurements. While the low-frequency technique provides Young’s moduli, |$E$|⁠, and Poisson’s ratios, |$\nu$|⁠, the ultrasonic technique gives P- and S-wave velocities. Therefore, a comparison of the obtained results requires conversion between the different parameters. To this end, this chapter provides description of the following issues: how the conversion between measured elastic parameters, |$E$| and |$\nu$|⁠, and wave velocities is done, related measurement errors, definitions of stress and strain sensitivities, explanation of how the stress sensitivity of Young’s modulus and Poisson’s ratio measured at seismic frequency is converted (under certain assumptions described below) into the stress sensitivity of P-wave velocities.

3.1 Linking seismic and ultrasonic measurements

Shales are commonly considered as transversely isotropic (TI) media. Stiffness tensor of a TI medium contains five independent stiffness parameters and can be written in the form of a 6 |$\times$| 6 matrix [using Voigt notation; Nye (1985)]:

where |$C_{11}$|⁠, |$C_{33}$|⁠, |$C_{44}$|⁠, |$C_{66}$| and |$C_{13}$| are the five independent stiffness parameters. On a single sample, however, the low-frequency apparatus allows for measuring only up to three independent elastic parameters (depending on sample orientation). To completely characterize TI media, measurements with three differently oriented samples (Fig. 2) are required. The parameters measured from the following sample orientations are

1. 0° sample – |$E_\mathrm{V},\nu _{\mathrm{VH}}$|⁠.

2. 45° sample – |$E_{45}$|⁠.

3. 90° sample – |$E_\mathrm{H}$|⁠, |$\nu _{\mathrm{HV}}$|⁠.

Here, the first index denotes the direction of uniaxial stress: V (vertical, perpendicular to bedding, this direction is also referred as ‘ax’—axial in the text), H (horizontal, parallel to bedding, can be also referred as ‘r’—radial) and 45 (under 45° to bedding); the second index in Poison’s ratio denotes the direction of transverse strain.

The measured elastic parameters can be converted into stiffness matrix parameters using the following equations (Bower 2010):

The stiffness parameters relate to wave velocities as follows:

Here, |$V_{\mathrm{PV}}$| and |$V_{\mathrm{PH}}$| denote velocities of P waves perpendicular to bedding and parallel to bedding, respectively, V|$_{SV}$| is the S-wave velocity perpendicular to bedding and |$V_{\mathrm{SH}}$| is the S-wave velocity parallel to bedding with the polarization within the bedding plane.

3.2 Errors

This section provides a summary regarding errors of the low-frequency measurement and data processing. For a comprehensive error analysis of the low-frequency apparatus, we refer to Lozovyi & Bauer (2019b).

3.2.1 Directly measured elastic parameters at low frequencies

The systematic error in Young’s modulus at low frequencies is in the range from −4 to −1.2 per cent, and the systematic error in Poisson’s ratio at low frequencies in the range from −0.016 to −0.008. Meaning that both quantities are systematically underestimated. Random error within a single test in Young’s modulus is |$\pm$|0.7 per cent, and in Poisson’s ratio is |$\pm$|0.008. Therefore, for different stress paths applied to a sample, relative changes of elastic parameters at low frequencies can be measured with high accuracy.

When the full set of core plug orientations is tested (Fig. 2), four of the low-frequency measured parameters are used to confirm TI symmetry of the medium (Holt 2016):

Deviations from unity of the left-hand term in eq. (7) will indicate deviations from TI symmetry or differences in properties of 0° and 90° oriented samples caused by heterogeneities (within a single samples and from one sample to another) as well as experimental error (from one measurement to another).

3.2.2 Conversion to velocities at low frequencies

Although the number of independent stiffness coefficients for TI media is five, the low-frequency setup allows for measuring six parameters (see Section 3.1). A least-squares fit routine additionally constrained by the requirement of positive elastic energy (Nye 1985) is used to estimate the five independent stiffness parameters from the six measured parameters (using relations from eq. 5). The estimated error in P-wave velocity at seismic frequencies (derived from the directly measured elastic parameters with differently oriented rock samples) can be within 10 per cent, but sometimes even higher for some specific rocks.

3.2.3 Ultrasonic velocities

The measured ultrasonic velocities have a systematic error of less than 1 per cent, and changes in velocity within a single test can be measured within less than 0.1 per cent (random error).

3.2.4 Elastic parameters obtained from ultrasonic velocities

The systematic error in inverted Young’s moduli from P- and S-wave velocities measured in axial direction for 0°, 45° and 90°-oriented samples (using eqs 5 and 6) may be up to 10 per cent, depending on the degree of rock heterogeneity. The corresponding errors for Poisson’s ratios are within |$\pm$|0.07. The errors are particularly sensitive to uncertainties in velocity measured with the 45° sample (as pointed out by Lozovyi & Bauer (2019b)).

3.3 Verification of stress sensitivity measurements

To ensure that the test system does not influence measured stress sensitivities of elastic parameters at seismic frequencies, a reference experiment was conducted using a linear-elastic reference sample. The experiment revealed no significant dependence of Young’s modulus or Poisson’s ratio on stress, with observed changes falling within an experimental error. Details and results of the verification test are provided in the Appendix  A.

3.4 Stress and strain sensitivities

For the stress and strain sensitivity analysis, the following parameters are acquired in the tests: changes in stresses (⁠|$\Delta {\sigma }$|⁠), strains (⁠|${\Delta {\varepsilon }}$|⁠), ultrasonic velocities (⁠|${\Delta {V}}$|⁠) and low-frequency moduli (⁠|${\Delta {E}}$| and |${\Delta {\nu }}$|⁠). In addition, reference values of ultrasonic velocity and low-frequency stiffness at initial stress state are needed. In the data analysis, strains, ultrasonic velocities, and low-frequency stiffness parameters are calculated for the unloading parts of the undrained cycles. The reason is that some irreversible rock deformations may occur during the loading stage that sometimes results in higher stress sensitivity than for unloading stage (Johnston 1987). For both T and M shales, a repeated isotropic loading cycle was performed at the end of the test in order to confirm that the properties of the sample were not affected in a significant way by the different stress cycles. For the data analysis, results obtained from the two isotropic stress cycles were averaged.

For the analysis of 4-D seismic data, the vertical-strain sensitivity of wave velocities (R-factor) is often used. The strain sensitivity (R-factor) is expressed as follows (Hatchell & Bourne 2005; Røste et al. 2006):

Holt et al. (2018) define vertical-stress sensitivity of wave velocities (S-factor):

For shales under undrained conditions, the stress sensitivity (S-factor) exhibits a linear dependence on stress-path parameter, |$\kappa$|⁠. In contrast, strain sensitivity (R-factor) demonstrates a nonlinear dependence on |$\kappa$| with strain sensitivity effect typically being more pronounced (Holt et al. 2018). Stress and strain sensitivities are linked by:

where the ratio |$\Delta {\sigma }_{\mathrm{ax}}/{\Delta {\varepsilon }_{\mathrm{ax}}}$| is a stiffness parameter taking different forms depending on the stress-path parameter |$\kappa$| - for triaxial loading (⁠|$\Delta {\sigma }_{\mathrm{r}}=0$|⁠), |$\Delta {\sigma }_{\mathrm{ax}}/{\Delta {\varepsilon }_{\mathrm{ax}}}= E_\mathrm{V}$| (vertical Young’s modulus), and for the uniaxial-strain stress-path parameter (⁠|$\Delta {\varepsilon }_{\mathrm{r}}=0$|⁠), |$\Delta {\sigma }_{\mathrm{ax}}/{\Delta {\varepsilon }_{\mathrm{ax}}}= C_{33}$| (vertical P-wave modulus).

To analyse the results of the low-frequency measurements we define vertical stress sensitivities of Young’s modulus (⁠|$S_{E}$|⁠) and Poisson’s ratio (⁠|$S_{\nu }$|⁠), directly measured at seismic frequencies:

3.5 Estimating high- and low-frequency strain sensitivity parameters

In our analysis, we focus on the stress sensitivity of the vertical P wave |$V_{PV}$|⁠, as this is the most relevant parameter for most geophysical applications. The vertical P-wave modulus, |$C_{33}$| (eq. 6), is estimated differently for the ultrasonic and low-frequency experimental data. In the ultrasonic technique, we measure P-wave velocities directly, and hence |$C_{33}$| can be estimated immediately without additional conversions. In the low-frequency technique, |$C_{33}$|⁠, can be also measured directly at seismic frequencies by synchronized modulation of both the axial and radial stress suppressing radial strain (Lozovyi & Bauer 2019b). However, this technique is not suited for measuring the stress sensitivity of |$C_{33}$|⁠, as the measurement error is larger than typical stress-induced changes in |$C_{33}$|⁠. Instead, we estimate |$C_{33}$| from dynamic Young’s moduli and Poisson’s ratios, which can be measured with high accuracy, making it possible to detect even small stress-induced parameter changes between the measurements taken at the reference and modified stress states.

The general algorithm for estimating the anisotropic strain-sensitivity parameters at low (seismic) and high (ultrasonic) frequencies is visualized in Fig. 4. High-frequency third-order elastic coefficients and directly computed from ultrasonic wave velocity measurement. Low-frequency third-order elastic coefficients are obtained using a rock physics inversion method. Inverting eq. (5), we relate |${C}_{\mathrm{ij}}$| elastic parameters to Young’s modulus and Poisson’s ratio measured for the 0° sample that is used for stress sensitivity study:

The system of eqs (12) and (13) is underdetermined, making it necessary to utilize a strain-dependent TI model (Duda et al. 2020), which allows for estimation of the strain sensitivity of the tested anisotropic material. The model is based on the third-order elasticity theory (Thurston 1974; Prioul et al. 2004; Fuck & Tsvankin 2009) and its general form is:

which is expanded into (using eight independent third-order parameters |$c_{\mathrm{ijk}}$|⁠)

where the compound coefficients are (Duda et al. 2020)

In the high-frequency limit, we estimate only the two relevant |$c_{\mathrm{ijk}}$| coefficients, that is |$c_{133}^{\mathrm{high\ f}}$| and |$c_{333}^{\mathrm{high\ f}}$|⁠, directly from the static strain and the ultrasonic P-wave measurements using eqs (6) and (15). In the low-frequency limit, the changes of the measured parameters are expressed as

where parameters |$E_{\mathrm{V}}^0$| and |$\nu _{\mathrm{VH}}^0$| are the averaged parameters repetitively measured at the reference stress state and |$\Delta E_{\mathrm{V}}$|⁠, and |$\Delta \nu _{\mathrm{VH}}$| are differences in low-frequency dynamic parameters measured before and after the loading/unloading step.

To carry out the inversion of the third-order elastic parameters, eqs (17) and (18) are combined with eqs (12), (13) and (14). Due to significant discrepancies in the values of Poisson’s ratio and Young’s modulus, the inversion routine, based on the Nelder–Mead method (Nelder & Mead 1965), requires normalization and an additional empirically estimated parameter to balance the contributions of these two parameters to the penalty function. The inversion method is described in more detail in Appendix  B.

4 RESULTS

4.1 Dispersion

Young’s moduli, |$E_{\mathrm{V}}$|⁠, and Poisson’s ratios, |$\nu _{\mathrm{VH}}$|⁠, were obtained in the experiments and plotted as a function of frequency (Fig. 5) for both tested shales at their initial stress (given in Table 1). Young’s modulus dispersion within the seismic frequency range (from 1 to  100 Hz) of about 5–6 per cent is similar for both shales. Poisson’s ratio in the same frequency range is almost unchanged for both shales. Young’s modulus dispersion from lower seismic (1 Hz) to ultrasonic (500 kHz) frequency is about 30 per cent for T shale and about 60 per cent for M shale. Poisson’s ratio dispersion from 1 Hz to 500 kHz is +0.02 for T shale and −0.05 for M shale. Vertical P-wave velocity results plotted as a function of frequency (Fig. 6) indicate that M shale does not exhibit any significant dispersion for the whole accessible frequency range, while T shale exhibits 10 per cent increase in velocity from 1 Hz to 500 kHz. Results of seismic and ultrasonic measurements at the initial stress state are summarized in the form of stiffness matrix parameters, |$C_{\mathrm{ij}}$| (Table 2).

4.2 Stress and strain sensitivity

The directly estimated high-frequency third-order elastic coefficients |$c_{133}^{\mathrm{high\ f}}$| and |$c_{333}^{\mathrm{high\ f}}$| are 449 and 485 GPa for T shale, and 139 and 148 GPa for M shale. The normalized stress sensitivities (S-factors, eq. 9) for axial ultrasonic P-wave velocities of T and M shales are plotted as a function of stress-path parameter (Fig. 7). Measured stress sensitivities are compared to fitted stress sensitivities given by the non-linearly elastic model for transversely isotropic medium. The chosen non-linear model (Duda et al. 2020) provides a good fit to the experimental data for both shales, giving an increase in stress sensitivity values with stress path. The small deviations between the measured and the modelled values might have been explained by velocity and strain measurements errors and by potential accumulation of non-elastic strains.

To obtain low-frequency third-order parameters |$c_{\mathrm{ijk}}^{\mathrm{low\ f}}$| for T and M shales inversion procedure described in Section 3.5 was used. In the case of T shale, three stress paths were used for the inversion (⁠|$\kappa =\lbrace 1, 0, -0.5\rbrace$|⁠). In the case of M shale, the results of the low-frequency experiment consist of only two stress paths (⁠|$\kappa =\lbrace 0, -0.5\rbrace$|⁠). However, Duda et al. (2020) indicates that data points obtained for stress-path parameters |$\kappa =-0.5$| and |$\kappa =0$| are suitable to estimate the third-order elastic parameters and to model behaviour of the rock in a wider stress paths range. All eight (four per sample) resulting distributions of |$c_{\mathrm{ijk}}^{\mathrm{low\ f}}$| parameters can be accurately approximated with normal distribution curves (Fig. 8).

The parameters of the distributions of coefficients |$c_{133}^{\mathrm{low\ f}}$| and |$c_{333}^{\mathrm{low\ f}}$| relevant for modelling of |$C_{33}$| are summarized in Table 3. All the estimated statistical measures of the low-frequency third-order elastic tensor elements are larger than the high-frequency third-order parameters estimated directly from the ultrasonic measurements. These results indicate larger strain sensitivity of the low-frequency seismic velocities. Small differences between the average, the median and the mode values confirm the natural character of the distribution and suggest a strong tendency of data to cluster in the direct proximity of these central values. For later analysis, we will use a range of the inverted |$c_{133}^{\mathrm{low\ f}}$| and |$c_{333}^{\mathrm{low\ f}}$| parameters given by their average values and corresponding standard deviations, that is |$c_{\mathrm{ijk}}^{\mathrm{ave}} \pm s_{\mathrm{ijk}}$|⁠.

The results indicate that T shale exhibits strain sensitivity of |$C_{33}$| several times larger than M shale at both high and low frequencies (Table 3). At high frequency values of |$c_{133}$| and |$c_{333}$| are comparable, but at low frequency |$c_{133}$| is 1.5 and 2 times higher than |$c_{333}$| for T and M shales, respectively. Furthermore, at low frequency, similar qualitative and quantitative relationships between the third-order elastic parameters are held for both rocks. In general, parameters |${c}_{133}$| are smaller than |${c}_{333}$| and their distributions are of comparable width. For both rocks, parameters |${c}_{113}-{c}_{336}$| are smaller than |${c}_{111}^{*}-{c}_{166}^{*}$|⁠, and in both cases parameter |${c}_{111}^{*}-{c}_{166}^{*}$| seems to be the least constrained (i.e. has the largest standard deviation) of all inverted coefficients.

To check how well the model using the average values of the inverted third-order elastic coefficients |$c_{\mathrm{ijk}}$| has captured the experimental low-frequency data, the measured seismic stress sensitivities of Young’s moduli and Poisson’s ratios are plotted together with the modelled stress sensitivities (see Fig. 9). We obtained a good fit of modelled and experimentally estimated sensitivity values for both shales and observed that for Young’s moduli the sensitivity increases, whereas for Poisson’s ratio the sensitivity decreases with stress path.

To visualize the difference in response to applied stress paths between the frequency limits, we plotted measured high-frequency and estimated low-frequency |$C_{33}$| changes for 5 MPa change in axial stress (Fig. 10). The low-frequency stiffness changes were estimated with the entire range of the third-order elastic parameters given by their average values and their standard deviation (Table 3). For both shales, the estimated changes of the low-frequency stiffness parameter |$C_{33}$| are significantly larger than their measured high-frequency counterparts. Even though the standard deviation of the inverted low frequency |$c_{133}$| and |$c_{333}$| parameters are comparable for both shales, in |$\Delta C_{33}$| results M shale have significantly higher standard deviation error due to its lower stiffness that results in higher strain yielded to the same amount of stress compared to T shale.

Next, we use |${c}_{133}$| and |${c}_{333}$| to estimate stress sensitivities of vertical seismic P-wave velocity (eq. 9), as shown in Fig. 11. The stress sensitivity of velocities for both T and M shales is larger in the seismic frequency range. On average, the ratio between seismic and ultrasonic velocity stress sensitivity of T shale is 2.6, and for M shale, the average ratio is 3.7.

5 DISCUSSION

5.1 Aspects of low-frequency laboratory data inversion for seismic velocities

The low-frequency experiments, inversion of the third-order elastic coefficients and subsequent modelling of the vertical P-wave velocity changes indicate high consistency of the results. First, for both considered rocks and all analysed stress paths, the strain sensitivity of the low-frequency signal is higher than the sensitivity of ultrasonic signal velocities. In both rocks, this relative discrepancy between the strain sensitivities increases with decreasing stress-path parameter value (Fig. 12). Secondly, the distribution of the inverted third-order elastic parameters |$c_\mathrm{ijk}$| indicate high robustness of the method and good quality of the results—the approximately normal distribution indicates the existence of a global minimum of the objective function, and hence of a probable set of unique physical values of parameters |$c_\mathrm{ijk}$| describing the behaviour of the rock. This observation is supported by small differences between mean, median and modal values of the inversion results. Third, the low-frequency strain sensitivities (R-factors) obtained with mean values of the third-order elastic coefficients (⁠|$R$| = 6–34) are consistent with the strain sensitivity (R-factor) values derived from field studies (⁠|$R$| = 4–35) reported by MacBeth et al. (2018).

The complexity of the inversion workflow and introduction of the additional penalty coefficient in eq. (B5) remain the subject of discussion. We consider improving the quality of the direct experimental measurements of the low-frequency P-wave modulus to be the most promising candidate to simplify and generally improve the data processing routine. In this case, parameters |${c}_{133}$| and |${c}_{333}$| would be estimated directly from measured changes of |${C}_{33}$|⁠, as in eq. (15), and not indirectly through |$\Delta E_{\mathrm{V}}$| and |$\Delta \nu _{\mathrm{VH}}$|⁠.

The use of the Gaussian distribution to describe the inversion results not only allows us to preliminarily verify their quality, but also offers an advantage in data processing. If we consider the average value of |${c}_\mathrm{{ijk}}$| parameters as the most representative for rock properties, the distance from this central value given by the standard deviation logically serves as a candidate for a quality threshold, thereby limiting the number of considered solutions. Such standard deviation analysis enables a deeper exploration of the relationships between |${c}_\mathrm{{ijk}}$| parameters obtained for both studied overburden shales. The relative order and magnitude of average values, as well as the relative width of the Gaussian distributions, are similar for both rocks, thus supporting our claim that the results for both samples are comparable. To reduce the number of independent model parameters and further simplify the method, additional experimental studies are required. These studies might reveal further correlations between the second and third-order parameters of the rocks.

5.2 Relevance of the results for seismic inversion

The frequency of ultrasonic measurement, the most common technique to assess the dynamic elastic properties of rock samples in the laboratory, is several orders of magnitude higher than that used in field seismic surveys. Our work is aimed to bridge the gap between strain sensitivities (R-factors) estimated from ultrasonic experimental data (e.g. Bakk et al. 2024) and those inverted from field data (e.g. MacBeth et al. 2018). The accounting for frequency dependence of strain sensitivity would enable more detailed and nuanced 4-D seismic-based strain and stress inversion studies (such as those by Hatchell & Bourne 2005; Herwanger & Horne 2005, 2009; Hodgson 2009; Røste et al. 2015; Røste & Ke 2017). These studies could then account for differences between various rocks in the overburden and the impact of stress path variations above the reservoir, as postulated by Toomey et al. (2017). Moreover, by using the results obtained with high-frequency parameters by Bakk et al. (2024), the analysis of the low-frequency strain and stress velocity sensitivity can be extended to examine the impact of the offset and azimuth of the seismic data. This analysis would further improve the quality of the links between seismic and geomechanical data.

Our results indicate that the stress sensitivities of P-wave velocities are systematically higher at seismic frequencies than at ultrasonic frequencies. This finding is consistent with previous laboratory results (Szewczyk et al. 2017), where a difference between ultrasonic and seismic stress sensitivities by a factor of up to 2.5 was reported for partially saturated Mancos shale under an isotropic stress path. Our empirical comparison of the measured ultrasonic and the estimated seismic sensitivities suggests a two to four times higher sensitivity of the low-frequency velocities. The test conditions were as close as possible to in situ: samples were preserved, subjected to equivalent in situ net stresses, and undrained conditions were applied. The effect of temperature was not considered in this work, as laboratory measurements were performed at room temperature. The test protocol consisted of subjecting the samples to different stress paths with |$\kappa$| ranging from −0.5 to 1 in order to simulate the stress conditions of an overburden in a realistic reservoir depletion scenario. For the isotropic stress path, the estimated seismic stress sensitivity is nearly 2 times higher than the ultrasonic for both T and M shales (Fig. 12). As stress-path parameter |$\kappa$| decrease (deviatoric stress increase), the difference between the seismic and ultrasonic stress sensitivity rise to an average factor of 4 for constant mean stress condition.

Stress sensitivity parameter (S-factor), the primary subject of our study, is a practical parameter for predicting the velocity change for any given stress path due to observed linearity of stress sensitivity with respect to the stress path (Fig. 11). Using eq. (10), our results can be readily transformed into strain sensitivity (R-factor), through the application of static stiffness properties. In the realm of 4-D seismic data interpretation, strain sensitivity (R-factor) is often derived and can be a useful tool for addressing changes in overburden stress. The strain sensitivities (R-factors) for both T and M shales increase non-linearly from |$\kappa$| of −0.5 to 1 at seismic and ultrasonic frequencies (Fig. 13). Strain sensitivity (R-factor) is a non-linear function of the stress-path parameter (Holt et al. 2017) due to its dependence on the strain, which itself is stress path dependent. Strain sensitivity (R-factor) for T shale is generally higher than that for M shale, which is opposite than for stress sensitivity (S-factor) where M shale exhibit higher stress sensitivity than T shale. The higher strain sensitivity for T shale is explained by much higher stiffness of T shale—the same amount of stress change result in significantly smaller strains that, according to eq. (15), results in less change in |$C_{33}$| and corresponding P-wave velocity.

The presented results suggest that inversion of 4-D seismic data by applying stress sensitivities obtained from ultrasonic measurements in the laboratory could result in significant overestimation of stress changes (or strains) in the overburden if the frequency dependence of stress sensitivities is ignored. Stress changes pose a threat for the geomechanical integrity of the cap rock situated above reservoirs, whether they are depleting or inflating. Additionally, if ultrasonic stress sensitivities of velocities are used for 4-D-seismic feasibility studies, estimated time-lapse time-shifts would be underestimated, which could potentially result in the conclusion that 4-D seismic is not feasible, while in reality, time-shifts could be high enough to be detected in 4-D seismic surveys.

5.3 Frequency-dependent wave velocities

Our measurements reveal that shales exhibit elastic dispersion within frequencies between 1 and 100 Hz and further up to ultrasonic frequency (500 kHz). Holt et al. (2023) argued that the mechanisms responsible for seismic dispersion in shales are complex and go beyond established dispersion theories like Biot flow, squirt flow or patchiness, and may be attributed to viscoelasticity of adsorbed (or bound) water on mineral surfaces and inside clay mineral platelets (Holt & Fær 2003; Holt & Kolstø 2017). The transition frequency between low- and high-frequency regimes in shales appear to be between |$10^3$| and |$10^5$| Hz (Duranti et al. 2005; Szewczyk et al. 2018), but since no continuous measurements are available over the entire frequency range, it is not possible to be more accurate for any specific shale. Squirt (or local) flow occurs between wide pores and low aspect ratio features like micro-cracks, regions between asperities at grain contacts, or compliant and elongated pores. Pores in shales have aspect ratios 0.02–0.05 (Xu & White 1995). The transition frequency |$f_{\mathrm{c}}^{\mathrm{sqt}}$| for squirt flow is |$f_{\mathrm{c}}^{\mathrm{sqt}}\approx (2/\pi ) (K/\eta ){{\alpha }^{3}}$|⁠, where |$\alpha$| is aspect ratio, |$\eta$| is viscosity, and |$K$| is the bulk modulus of the background medium (O’Connell & Budiansky 1977). With |$K$| in the range 10–30 GPa and water viscosity of 0.001 |$\text{Pa}\cdot {}\text{s}$|⁠, the aspect ratio needs to be in the vicinity of |$10^{-3}$| for dispersion to fall within the observed frequency range. Larger aspect ratio cracks or pores lead to transition above the ultrasonic frequency used in the laboratory. It is however known that water in the 2–3 molecular layers closest to solid surfaces may have strongly enhanced viscosities (Antognozzi et al. 2001; Goertz et al. 2007; Ulcinas et al. 2011): In several experimental studies water viscosities from 10 to |$10^{7}$|  P have been measured within 1 nm from mica surfaces. Thus, for pores in shales with |$0.02 \lt {\alpha } \lt 0.05$|⁠, the transition frequency for squirt flow may easily fall in the range 100 Hz–10 kHz. Clearly, a spectrum of aspect ratios as well as a spectrum of pore sizes may smear out this transition, so that it may extend from the seismic to the ultrasonic range. In addition, the smaller pores may also, if filled by highly viscous adsorbed (‘bound’) water, be sources of another dispersion mechanism known as viscous shear relaxation. The associated transition frequency is |$f_{c}^{\mathrm{visc\ shear}}\approx (1 / 2\pi ) (G/ \eta ) \alpha$| were |$G$| is the shear modulus of the host material surrounding the pore. With |$G$| between 5 and 15 GPa, the transition frequency for viscous shear relaxation may be from seismic to ultrasonic frequencies for cracks as well as normal pores with thickness smaller than 2–3 nm and water viscosities |$10^3 - 10^5$|⁠. Thus, the viscoelastic nature of bound water may explain the dispersion in shales. For stress excitations above the transition frequency (e.g. at ultrasonic frequencies), the water molecules at grain contacts and in low aspect ratio pores become immobile, which means that the bound water effectively acts as ‘cement’ between the grains giving a high rock stiffness. At low frequencies (e.g. at seismic frequencies) local flow of water occurs between low aspect ratio pores and larger and stiffer pores during stress modulations, resulting in reduced rock stiffness. In shales, the local flow is controlled not only by the aspect ratio of the compliant pores, but also by their size: the smallest ones will contain the water with highest viscosity and will contribute most to the velocity increase at high frequencies.

5.4 Frequency-dependent stress sensitivity of wave velocities

Building on the discussion of frequency-dependent wave velocities, we discuss physical mechanisms that can be responsible for frequency-dependent stress sensitivity. Under undrained conditions, stress sensitivity is mainly controlled by soft elongated pores and microcracks, which also, as discussed above, contribute to dispersion by squirt flow and viscous shear relaxation. Conventional crack models (Fjaer et al. 2021), in general form, predict relative changes in stiffness as

where |$Q_{\mathrm{ij}}$| is a coefficient given by the properties of the solid frame material, |$\xi$| is crack density parameter that relate to size and number of cracks (Hudson 1981), |$D$| is a drainage parameter, which depends on fluid and solid properties, aspect ratio of elongated cracks and is also controlled by frequency (Fjaer 2006). We assume all pores are oriented with their normal along the symmetry axis, which is a reasonable assumption for layered rock as shale, so that |$\xi \equiv \xi _{\mathrm{ij}}$|⁠. For low frequencies, fluid flow between cracks and surrounding pores is permitted (Thomsen 1985), and |$D \rightarrow 1$|⁠. For high frequencies, there is less time available for fluid flow, and in the ultimate limit, |$D \rightarrow 0$|⁠. The change in crack density is controlled by stress change (Fjaer 2006):

where |$\xi _{0}$| is the initial crack density, |$\Delta \sigma ^{\prime }$| and |$\Delta \tau$| are changes in effective normal and shear stress, respectively, |$a$| and |$b$| are stress sensitive coefficients that can be calibrated from the laboratory experiments. The effective stress change, |$\Delta \sigma ^{\prime }$|⁠, controlling crack density is assumed equal to the net stress (eq. 2) change :

The undrained pore pressure change is given by Skempton’s law, containing pore pressure parameters |$A$| and |$B$| (Skempton 1954). The combination of the stress changes (eq. 21) and Skempton’s equation leads to stress path dependency of the stress sensitivities (Holt et al. 2018). For isotropic stress change, relative change of P-wave modulus can be derived as

In this case, the pore pressure change is given by Skempton’s |$B$| through |$\Delta p_{\mathrm{f}}=B\Delta \sigma$|⁠, so the resulting stress sensitivity is

Because |$D$| decreases when frequency is increased from seismic to ultrasonic, hence the stress sensitivity given by |${\Delta C_{33}}/{C_{33}}$| also decreases. Under a non-isotropic stress path, changes in pore pressure will differ, and the effects of shear stresses will also become apparent. (eq. 20). Different values of the coefficients |$a$| and |$b$| will lead to different stress sensitivity. Both positive and negative slopes of stress sensitivity vs. stress-path parameter |$\kappa$| (as observed in Fig. 11) are possible.

Accounting for the impact of pore fluid bulk modulus and bulk density change as a result of pore pressure change, it is evident that the effective stress coefficient controlling P- and S-wave velocities will be different from 1. For P waves, the effective stress coefficient is |$\lt 1$|⁠, while for S waves the effective stress coefficient may be smaller but, in some cases, also larger than 1. Mavko & Vanorio (2010) demonstrated theoretically that the squirt flow dispersion mechanism causes the effective stress coefficient to depend on frequency of wave velocity. For a transition from ultrasonic to seismic frequency, the effective stress coefficient decreases and becomes closer to 1. The decrease of the effective stress coefficient also increases stress sensitivity, as long as we relate it to the total stress change.

6 CONCLUSIONS

We quantified stress-induced wave velocity changes for seismic and ultrasonic frequencies in laboratory tests conducted on shale cores. We directly recorded changes in elastic moduli at seismic frequencies and velocity changes at ultrasonic frequency as the samples underwent stress variations with different combinations of radial/axial stress changes. These combinations were designed to mimic stress paths in different overburden locations. At seismic frequencies, velocities could not be measured directly, therefore we applied a third-order elasticity model to derive P-wave velocity changes from changes in dynamic Young’s moduli and Poisson’s ratios measured at seismic frequencies.

The tested shales exhibited significant differences in stiffness, porosity, and stress sensitivities. However, the samples showed comparable dependence of velocity stress- and strain-sensitivities on frequency. Our results indicate a factor of 2–4 higher stress sensitivity of P-wave velocities at seismic frequencies than at ultrasonic frequencies for all applied stress paths. The same results apply to the strain sensitivity (R-factor). We also observed that the ratio of seismic/ultrasonic stress and strain sensitivities systematically depends on the stress path. However, stress sensitivities, unlike strain sensitivities, can be sufficiently well generalized by linear relation between the stress-path parameter and stress sensitivity at both seismic and ultrasonic frequencies. Therefore, in the context of rock physics modelling, the use of stress sensitivity (S-factor) is more practical than the use of strain sensitivity (R-factor).

Our experimental observations, supported by theoretical considerations, suggest that for saturated shales, the stress sensitivity may be higher at seismic frequencies than at ultrasonic frequencies, however, further studies are needed to confirm this. We recommend to consider the effect of frequency in geophysical analysis, in particular, for the inversion of time-lapse seismic data for strains and stress changes in shale formations.

ACKNOWLEDGEMENTS

The authors would like to acknowledge the following contributions: Stian Rørheim (for setting up a reference test on PEEK), Jørn Stenebråten (for support in laboratory testing) and Audun Bakk (for discussions). The authors acknowledge financial support from The Research Council of Norway (Grant no. 234074/E30), AkerBP, DONG Energy, ENGIE, Maersk and Total through the PETROMAKS 2 KPN-project: Shale Rock Physics: Improved seismic monitoring for increased recovery at SINTEF. This work was supported by the Research Council of Norway, AkerBP, Equinor, Shell and Vår Energi (Grant no. 294369; PETROMAKS 2) through the research project Improved Prediction of Stress and Pore-pressure Changes in the Overburden for Infill Drilling. The authors also acknowledge funding from the Research Council of Norway through PETROMAKS 2 program researcher project (Grant no. 301910).

AUTHOR CONTRIBUTIONS

S. Lozovyi (Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing–original draft, Writing–review & editing); M. Duda (Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing–original draft, Writing–review & editing); A. Bauer (Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing–original draft, Writing–review & editing) and R. M. Holt (Conceptualization, Funding acquisition, Investigation, Resources, Supervision, Writing–original draft, Writing–review & editing).

DATA AVAILABILITY

Data supporting this study are openly available at https://zenodo.org/records/12781507.

References

APPENDIX A: VERIFICATION OF LOW-FREQUENCY MEASUREMENTS

To verify that observed low-frequency changes of elastic parameters in rock are attributable to the material’s stress sensitivity rather than the test system, a reference experiment was conducted using a polyether ether ketone (PEEK) sample. We chose PEEK because its response can be approximated by linear elastic behaviour for the stress changes applied in our experiment (Hamdan & Swallowe 1996), implying that its stiffness remains constant regardless of stress magnitude or stress path. Following a similar test protocol as outlined in Fig. 3, the experiment encompassed isotropic, triaxial, and constant mean stress loading-unloading cycles with a stress magnitude of 5 MPa. Seven low-frequency measurement sweeps (1–144 Hz) were performed: one measurement for each of the three loaded states, and four measurements before and after each stress cycle. The reference experiment did not reveal systematic stress or stress-path dependence for Young’s modulus (Fig. A1a) or for Poisson’s ratio (Fig. A1b). Normalized changes in Young’s modulus (⁠|$\Delta E/E/\Delta {\sigma }_{\mathrm{ax}}$|⁠) and Poisson’s ratio (⁠|$\Delta \nu /\nu /\Delta {\sigma }_{\mathrm{ax}}$|⁠) are within an experimental error of 1.5 GPa|$^{-1}$| for both quantities (Figs A1c and d).

APPENDIX B: INVERSION OF LOW-FREQUENCY PARAMETERS

Eqs (17) and (18) together with (12), (13) and (14) yield a set of two non-linear equations describing changes in stiffness for unloading segment of each stress-path parameter |$l$|⁠, which we normalized by the difference between the largest and the smallest measured values of their corresponding parameter:

where |$C_{\mathrm{ij}}^{0}$| are estimated through measurements on three differently oriented samples (0°, 45° and 90°), |$\Delta C_{\mathrm{ij}}$| are given by eqs (14)–(16), and |$a$| = |$b$| = 0. Stiffness coefficients |${C}_{11}$| and |${C}_{66}$| are grouped together, and hence they are treated as a single effective parameter, yielding:

Eqs (15), (B1), (B2) and (B3) indicate that four third-order coefficients are needed to describe strain sensitivity of |${E}_{\mathrm{V}}$| and |${\nu }_{\mathrm{VH}}$|⁠, that is |${c}_{133}^{\mathrm{low f}}$|⁠, |${c}_{333}^{\mathrm{low f}}$|⁠, |$c_{111}^{\mathrm{*,low f}}-c_{166}^{\mathrm{*, low f}}$| and |${c}_{113}^{\mathrm{low f}}-{c}_{366}^{\mathrm{low f}}$|⁠. However, although we need to estimate all of them to obtain a physically realistic set of third-order elastic coefficients, we need only |${c}_{133}^{\mathrm{low f}}$| and |${c}_{333}^{\mathrm{low f}}$| to model the behaviour of |$C_{33}$| at low frequencies with eq. (15).

The objective function |$J$| minimized by the inversion algorithm takes form of a least squares error

where |$N$| is the total number of considered loading or unloading stress-cycles (isotropic, triaxial or constant mean stress), |$a^l$| and |$b^l$| are given by the eqs (B1) and (B2), respectively, and |$Y^l$| is a penalty function as follow

The penalty function is introduced since we look for the smallest possible |$\Delta C_{\mathrm{ij}}$| explaining observed |$\Delta E_{\mathrm{V}}$| and |$\Delta \nu _{\mathrm{VH}}$|⁠. The coefficient |$2.6 10^{-10}$| in eq. (B5) was determined empirically—it has the highest possible value allowing for replication of |$\Delta E_{\mathrm{V}}$| and |$\Delta \nu _{\mathrm{VH}}$| (too high value of this coefficient would not allow |${C}_{ij}$| to change enough to replicated the measured values) which also keeps the distribution of the results free of artifacts (i.e. sharp distribution edges—details in Section 4.2).

We use Monte Carlo method to generate starting values of the optimized |$c_{\mathrm{ijk}}$| coefficients required to initiate the Nelder–Mead inversion. To determine the range of possible starting values, we use the high-frequency third-order elastic coefficients estimated directly from the ultrasonic data. For |$c_{133}$| and |$c_{333}$|⁠, Monte Carlo method generates initial values ranging from 0.5 to 2.0 of their high frequency values, that is 0.5|${c}_{133}^{\mathrm{high f}}$|–2|${c}_{133}^{\mathrm{high f}}$| and 0.5|${c}_{333}^{\mathrm{high f}}$|–2|${c}_{333}^{\mathrm{high f}}$|⁠. For the remaining two parameters, we have no a priori information on their high-frequency values, therefore we assume that they do not differ significantly from |${c}_{133}^{\mathrm{high f}}$| and |${c}_{333}^{\mathrm{high f}}$| (Duda et al. 2020). Hence, Monte Carlo method is limited by the average of the known high-frequency parameters, that is generates starting values between 0.25 |$({c}_{133}^{\mathrm{high f}}+{c}_{333}^{\mathrm{high f}})$| and |$({c}_{133}^{\mathrm{high f}}+{c}_{333}^{\mathrm{high f}})$|⁠. We run the inversion 500 000 times. We limit the number of considered solutions, that is sets of low-frequency third-order elastic coefficients, using a relative error and a coefficients distribution criteria (described in more details in the Sections B1 and B2).

B1 Inversion-results filtering: fit-error criterion

First, to remove the inversion results yielding inferior fit of the modelled changes of |$E$| and |$\nu$| to the experimentally measured values, we disregard 10 per cent of the inversion results (i.e. sets of |${c}_{133}^{\mathrm{low f}}$|⁠, |${c}_{333}^{\mathrm{low f}}$|⁠, |$c_{111}^{\mathrm{*,low f}}-c_{166}^{\mathrm{*, low f}}$| and |${c}_{113}^{\mathrm{low f}}-{c}_{366}^{\mathrm{low f}}$|⁠) with the highest fit error values.

B2 Inversion-results filtering: coefficients-distribution criterion

The distribution of the remaining low-frequency |$c_{\mathrm{ijk}}$| coefficients was approximated by an average (⁠|$c_{\mathrm{ijk}}^{\mathrm{ave}}$|⁠) and a standard deviation (⁠|$s_{\mathrm{ijk}}$|⁠) values, usually used to describe normal distributions. To visualize the behaviour of |$C_{33}$| in low frequencies we used the inversion results for which both |$c_{133}$| and |$c_{333}$| were within |$c_{\mathrm{ijk}}^{\mathrm{ave}} \pm s_{\mathrm{ijk}}$| limits, as we expect them to be the closest to the real physical values of the third-order elastic coefficients in the tested samples.

APPENDIX C: COMPLEMENTARY ULTRASONIC DATA SET

For T shale, the ultrasonic and the low-frequency measurements were taken during the same test. However, for M shale, low-frequency measurements were complemented with the ultrasonic data that come from a different test (the experimental setup and the sample are described in Holt et al. (2018). The complementary data set was used because the ultrasonic measurements during the isotropic stress cycle in the test in the low-frequency apparatus were not reliable, whereas the stress sensitivities for the triaxial and constant mean stress paths were comparable. Notice that the additional test on M shale was performed at the same net stresses but with higher absolute values of stresses and pore pressure, comparable to the in situ stress state. Additionally, the complimentary ultrasonic velocity measurements includes uniaxial strain stress-path parameter(⁠|$\kappa = 0.76$|⁠).

Both the ultrasonic test data and the low-frequency test data for M shale are supplemented with readings of strains in axial and radial directions. In data analysis, we use strains, Poisson’s ratio, and Young’s modulus measurements from the low-frequency experiments (stress-path parameters: |$\kappa$| = −0.5 and 0) to estimate the low-frequency third-order elastic coefficients |$c_{\mathrm{ijk}}$|⁠. Then, we use the estimated low-frequency parameters |$c_{133}^\mathrm{{low f}}$| and |$c_{333}^\mathrm{{low f}}$| together with strains measured during the high-frequency experiment to model potential changes of the low-frequency stiffness |$C_{33}$| during the high-frequency test. This methodology allows us to compare the low- and high-frequency stiffness changes between the two tests and to extrapolate the modelling results obtained for the low-frequency range to a wider range of stress paths—from |$\kappa = \lbrace -0.5, 0 \rbrace$| in the low-frequency experiment to |$\kappa = \lbrace -0.5, 0, 0.76 ,1\rbrace$| in the high-frequency test. This modelling step mimics a situation in which results from a low-frequency test are used to model the response of the rock in a field-development scenario in which the stress paths range can be significantly wider than in the stress conditions examined in laboratory conditions.