Laboratory constraints on the anisotropic dynamic-to-static ratios for shale's elastic constants: an example from the Duvernay unconventional reservoir

Summary

1 INTRODUCTION

1.1 Context and background

Studies on shale's geomechanical properties are currently in demand (Britt and Schoeffler 2009; Jin et al. 2014; Villamor Lora et al. 2016) due to the expansion of hydraulic fracturing (HF) operations to extract hydrocarbon resources that were previously locked within low permeability shales. HF essentially isolates a borehole segment into which the fracking fluids are injected. The injection causes overpressure that creates fractures growing into the geological materials surrounding the borehole. These fractures connect pockets of hydrocarbon deposits and create extra pore surface areas, resulting in greater hydrocarbon mobility.

Ensuring this process is economically viable (Montgomery & Smith 2010) requires adequate knowledge of the reservoir rocks’ mechanical constitutive models. To achieve ideal HF well completion for production, it is critical to maximize the fracture networks. Though optimization of the engineering operational procedures/techniques can often be beneficial, the geomechanical characteristics of the targeted reservoir still largely dictate the efficacy of fracture creation and subsequent hydrocarbon production following the well completion (Nagel et al. 2013; Singh et al. 2019). Consequently, HF production optimization and forecast require quantitative knowledge of rock masses’ elastic properties (e.g. Settari & Cleary 1984; Dean & Schmidt 2009). Moreover, relevant engineering issues in general subsurface resource development, that includes determination of in situ stress (e.g. Schmitt et al. 2012; Wang et al. 2022), borehole stability (e.g. Ong & Roegiers 1993; Zhang 2013), caprock integrity (e.g. Gor et al. 2013; Vilarrasa et al. 2014), etc., all require the knowledge of geomechanical characteristics of the rock masses.

It is well known that acoustic/seismic wave velocities are functions of material's elastic stiffness and density, and vice versa, wave velocities can be used to invert for these material constants. However, rock's static properties, which describe its response to a long-period, quasi-static external loadings, are generally different from this wave (e.g. seismic and acoustic waves) measured dynamic properties (Cheng & Johnston 1981; Fjær 2009). Subsurface engineering work often requires knowledge of the former static properties. In practice, however, the dynamic properties of the rock are usually more accessible and can be measured with acoustic well logs or via the inversion of active source seismic survey data. Naturally, workers have attempted to correlate the measurements from seismic surveys (e.g. Slatt & Abousleiman 2011; Perez & Marfurt 2013; Sohail et al. 2020) and acoustic loggings (e.g. Venieri et al. 2020, 2021) with rocks’ static elastic constants by introducing a dynamic-to-static ratio.

This demand is particularly urgent for the Canadian Duvernay Unconventional reservoir. The Duvernay Formation is a geological unit covers partially the westernmost part of the Western Canada Sedimentary Basin (Mossop 1988) and hosts one of the promisingly emerging shale oil/gas plays. The spatial extent of the Duvernay Formation covers an area of ∼130 000 km² (Fig. 1) with the net pay thickness reaching as much as 60 m in some parts (Preston et al. 2016). The Duvernay shale gas/oil were previously considered unexploitable owing to Duvernay shale's extremely low permeability (Ghanizadeh et al. 2015)); this deterrence is only recently overcome by the massive adoption of multistage horizontal HF.

Figure 1.

Brief geological context with (a) spatial extent of the Duvernay Formation with locations of the boreholes from where samples are extracted (denoted as green dots); and (b) a vertical cross-section for the volume of the crust where samples are retrieved (denoted with letters A-G and dots). This figure is adapted from Shen et al. (2021).

Open in new tabDownload slide

Despite rising production and drilling activities, at the time of this writing, the Duvernay Formation is still considered to be in its early development stage. Particularly, the East Shale Basin area is still largely at the exploration stage. Geomechanics-related development problems remain one of the capital-intensive issues and researches are urgently in demand.

Though efforts had been expanded, with many different approaches to study the rocks’ mechanical properties, several critical issues have been raised and require attention. The majority of the underground rock mass’ physical properties are anisotropic. A constitutive model built based on an isotropic assumption usually cannot describe rock's mechanical behaviours adequately (e.g. Schoenberg & Sayers 1995; Cho et al. 2012; Vasin et al. 2013; Chandler et al. 2016). On account of their layered structure, the elastic properties of shales are generally assumed to be transversely isotropic with a rotational symmetry axis perpendicular to its isotropic plane (Fig. 2). Past observations from lab experiments (e.g. Sone & Zoback 2013; Meléndez-Martínez & Schmitt 2016; Ong et al. 2016; Gong et al. 2019) have generally supported this assertion well. Further, rock's mechanical elastic ‘constants’ are generally stress-dependent. This is readily observable from the changing velocities of the waves traversing through the sample and the slopes of stress–strain curves of the samples under static loading with a variance of stress. Nevertheless, for small Δσ, the linear Hookean elasticity can sufficiently describe the material's deformation Δε in response to Δσ. Geomechanical measurements on these rocks need to be converted to a ‘compatible’ stress state before being compared.

Figure 2.

An illustration for the TI materials with laminating bedding planes represented by grey lines; and axes x₁–x₂–x₃ denote the coordinate system used in this work.

Open in new tabDownload slide

1.2 TI assumption for shale

This work focuses on measuring Duvernay shale's dynamic and static elastic constants. Here, we briefly recoup the theory of TI elasticity to build the essential context for readers to follow. As illustrated in Fig. 2, the 2-D isotropic plane of the TI material can be denoted by two orthogonal axes x₁ and x₂. The rotational symmetry axis (x₃) is perpendicular to this plane. These axes define a cartesian coordinate system x₁–x₂–x₃ in which TI material's elastic constitutive model can be described by eq. (1) with a total of five material constants (S_ij) known as the compliances:

S can also be transformed into a more often used expression by the petroleum/mining industry practitioners with Young's modulus E, Poisson's ratio υ, and shear modulus G:

and for completeness it needs to be noted that there are three Poisson's ratios (i.e. ν_hh, ν_vh, ν_hv) that may be defined but with the constraint ν_vh/E_v = ν_hv/E_h.

For many sedimentary basins, particularly marine shale plays, the symmetry axis (x₃) is often parallel to gravity. Therefore, under the same configuration, the axis x₃ are often alternatively denoted as the vertical (V) direction and directions parallel to the x₁–x₂ plane are denoted as the horizontal (H) directions.

We also note that, more usually, the material constants measured through wave methods (i.e. dynamic material constants) are reported as the stiffness:

where C simply equals to the inversion of S (i.e. S⁻¹). C_ij are functions of the velocities of different modes of waves travelling through the material at different incident angles (Thomsen 1986).

In this work, we focus on measuring sample's dynamic and static S_ij. We use superscripts D, that stands for dynamic, and S, that represents static, to distinguish them. R_ij refers to the dynamic-to-static ratio S^D_ij/S^S_ij. Some reports refer to such ratio as static-to-dynamic ratio that is defined as the dividend of the static and dynamic Young's modulus E: their static-to-dynamic ratio is equal to our dynamic-to-static ratio because E_V = 1/S₃₃ and E_H = 1/S₁₁ (see eq. 2). At last, in alignment with other geomechanical studies, compressive stresses and strains have a positive sign.

1.3 Objective and scope

In this work, we measure a set of Duvernay shale samples’ elastic (dynamic and static) constants, with full consideration for their anisotropy and pressure dependence and calculate their dynamic-to-static ratios. We also explore the relationship between various geomechanical/geochemical properties and the dynamic-to-static ratio. The practical goal of this work is to provide a set of measurements with interpretations that can facilitate derivation of the static elastic properties of the Duvernay unconventional reservoir at large.

2 EXPERIMENT AND RESULTS

2.1 Samples and preparation

This work conducts experimental measurements on samples extracted from Duvernay Formation cores (samples A, B, C, D and G, see Fig. 3 and details in Table 1).

Figure 3.

Samples cut into prisms ready for ultrasonic experiments. The picture for sample D is taken after the ultrasonic investigation with an ∼1.5-inch diameter H plug drilled for the static loading experiment.

Open in new tabDownload slide

Through sampling processing, we picked rocks that are (1) with no visible fractures and (2) as homogeneous as possible to avoid complications and uncertainties that may arise due to large fractures and sample's own heterogeneity. The detailed characterizations for these samples are included in the Supporting Information.

Samples tested here are cut from the cores. We machine them into the prism shape with each side of the sample having lengths of ∼3 cm to allow sufficient distances for reliable wave velocities measurements. At the direction perpendicular to the symmetry x₃-axis, a cylindrical H plug (1.5 inch diameter, 3–4 inches length) is drilled after the completion of the ultrasonic wave measurements to be subject to static loading tests; a V plug is drilled as close to the location where the prism shaped samples are machined as possible. The residual of the H plug, after its destruction through the static loading tests, are machined into thin sections or powders for geochemical and petrophysical analysis. All experiments are conducted in the samples’ ‘dry’ conditions after being placed in the oven for more than 24 hr at a temperature of less than 70 °C to ensure that no pore fluids are present during the tests. The change of rocks’ properties in this low temperature drying process is minimal and we dry our sample to (1) ensure our results are compatible to other studies that generally are conducted on dry samples; and (2) allow future studies of fluid substitution to account for varying fluid properties. Again, it is important to reiterate that these samples contain no detectable water sensitive clays.

2.2 Measuring dynamic material constants

The dynamic material constants are measured using the ultrasonic pulse transmission method. Matched pairs of piezoelectric ceramic wafers (resonant frequency of 1 MHz) are attached to the prisms in patterns of Figs 4(a) and (b). These piezoelectric wafers, which emit and receive ultrasonic waves, are activated by the sharp fast-rising voltage pulses generated by a pulser/receiver (JSR-PR35). The square wafer emits/receives the shear wave, and the round one emits/receives the compressional wave. We note one additional round wafer [see P45(1) in Fig. 4a] is attached in the 45° surface to account for the beam-skew effects for wave propagation directions away from the rock's symmetry axes (Li et al. 2018). At last, we paint the sample with the Flexane 80 putty that cures after 12 hr (Fig. 4c). The cured putty essentially jackets the samples preventing fluid intrusion when the samples are placed in the hydraulic cell (Fig. 4d) to be measured with hydrostatic pressures as high as 180 MPa.

Figure 4.

(a) Schematics of the hydraulic cell in which the sample is pressurized and measured. (b) Photo of a sample with electric wiring, piezoelectric wafers and strain gauges and (c) a photo taken after painting the sample with impermeable putty coating. (d) Picture of the actual apparatus with hydraulic cell providing up to 180 MPa hydrostatic pressure loading environment.

Open in new tabDownload slide

With the complete waveforms recorded (see example in Figs 5a–f), we determined these waves’ phase velocities (Fig. 5g). We then use equations that are developed on the basis of Thomsen's (1986) equations to obtain samples dynamic compliances S^D_ij (Fig. 5h) which are simply components of C⁻¹. Readers are referred to the Supporting Information for the details of this process.

Figure 5.

(a)–(f) Raw waveforms recorded on a sample for six pairs of piezoelectric wafers (measures P and S waves travel through the sample at different angles) during the entire pressurization/depressurization process; (g) measured wave velocities from each pair of piezoelectric wafers on which basis that (h) material constants S^D_ij are calculated.

Open in new tabDownload slide

2.3 Measuring static elastic constants

The first set of constraints on the static compliances of the sample are obtained on the prismatic sample under hydrostatic loading, at the same time when the dynamic compliances are measured. As shown in Figs 4(a) and (b), foil strain gauges are attached in the surface of our samples in the x₁ and x₃ directions. The strain gauges are wired into a Wheatstone bridge-like system (Melendez-Martinez 2014) that amplifies circuit's voltage change of due to the strain gauges’ varying resistivities in response to the sample's deformation. During the experiment, we increase the hydrostatic pressure at steps of 5 MPa when less than 40 MPa, and then 10 MPa when less than 140 MPa, and then 20 MPa till it reaches the maximum of 180 MPa. Measured strains are manually written down at each pressure point at the same time when ultrasonic waveforms are recorded by the computerized digitizer system (see example in Fig. 6).

Figure 6.

Strains measured in the x₁ (H) and x₃ (V) directions of the prism-shaped sample simultaneous with the ultrasonic wave velocities at hydrostatic pressures.

Open in new tabDownload slide

The slopes of Δε_iversu Δσ_H plots (see Fig. 6), which are defined as KL_i, are functions of the sample's static compliances S^S_ij. By plugging Δσ_hydrostatic = Δσ₁₁ = Δσ₂₂ = Δσ₃₃ and Δσ_ij(i≠j) = 0 into eq. (1), we can have:

and

In some literature KL₁ and KL₃ (e.g. Melendez-Martinez 2014; Ong et al. 2016) are referred to as the ‘linear bulk moduli’ and their inverse 1/KL_i are called the linear bulk compliances. This resembles the bulk modulus/compliances for isotropic material. But, for anisotropic materials, the material's strain under hydrostatic pressure loading is anisotropic. Therefore, we denote this ‘bulk’ moduli/compliance with subscript i that marks its directionality in alignment with the x₁–x₂–x₃ coordinate system shown in Fig. 2. KL^S_i or their inverse (linear static compliances 1/KL^S_i), cannot be used to infer for S^S_ij alone, but they can help validate our static S^S_ij measurements—we will expand on this in the later text.

The second set of measurements for the static compliances are made on two cylindrical plugs (diameter ∼1.5 inch, length ∼3–4 inch): one H plug is cored in the direction parallel to the isotropic x₁–x₂ plane; one V plug is cored parallel to the x₃ symmetry axis. We use Hoek Cell for the tri-axial compression loading tests (Figs 7a and b). We note that H plug is cored from the sample which had undergone dynamic compliances measurements (see example in Fig. 3, sample D) and V plug is cored from the larger bulk of the sample before it is machined down to the prism shape for dynamic compliances measurements. We attached strain gauges on the surfaces of the plugs in a patterns described in Figs 7(c)–(e).

Figure 7.

(a) Picture of the Hoek Cell used for triaxial loading tests and its (b) schematics; the schematics of sample preparation with strain gauges attached to the samples with (c) for a V plug; and (d) and e) for H plug. Pictures of H (f) and V(g) plugs from sample A with strain gauges attached after the static loading tests, samples are destroyed after the loading tests.

Open in new tabDownload slide

For H plug, we first set the radial pressure at a constant 0 MPa and increase the axial pressure while measuring strains and recording axial/raidal pressures; then we hold the axial pressure at a constant value (by shutting off the hydraulic oil valve and do not attempt to actively control it) and increase the radial pressure. At the same time, strains/pressures are measured and recorded.

The measurements made with rising axial pressure (that is equivalent of σ₁₁ in eq. 1) and fixed radial pressure (that is equivalent of σ₃₃ and σ₂₂) allows the calculation of S^S₁₂, S^S₁₁ and S^S₁₃ (see Figs 8c–h, e.g.):

Figure 8.

Stress/strain curves with (a) and (b) the complete loading procedures on the H plug, please note the strains/pressures are reported at a frequency of 1000 measurements per minute; (c)–(h) with axial loading pressure increasing and radial wall pressure set as constant and (i)–(m) with radial wall pressure increasing and not active effort in controlling the axial pressure.

Open in new tabDownload slide

Also, uniquely in this study, adopting the method described in Shen & Playter (2021), we release the constraints on the radial pressure and start increasing it, with no active efforts on controlling axial pressure though to obtain the constraints on S^S₃₃ as function of (1) axial pressure (σ₁₁), (2) radial pressure (σ₂₂) and (3) prior knowledge of S^S₁₃ obtained with eq. (6) in the pressure steps before and after (see Figs 8i–m, e.g.):

We also perform tri-axial loading experiment on the V plug: we record the stress versus strain data series with constant radial pressures and rising axial loading pressures (i.e. Δσ_R = Δσ₁₁ = Δσ₂₂ = 0) to obtain additional sets of measurements on S^S₁₃ and S^S₃₃ as functions of radial (σ₁₁ = σ₂₂) and axial (σ₃₃) pressures:

The strains, reported through linear strain gauges wired into a Wheaston bridge configuration, and pressures, measured by the pressure transducers, are recorded at a frequency of 1000 measurements per minute. We note that as we use hand pumps to increase the pressures (both axial and radial) through hydraulic systems, we aim to set each pressure step at 10 MPa with a loading rate of roughly half-minute per step. Please note the kinks in the stress–strain plots (Fig. 8) are the results of the operator changing hands when using manual pumps to increase the pressure of the hydraulic press. More details on the derivation of eqs (6)–(9) based on which we constrain S^S_ij can be seen in Shen & Playter (2021) and also the attached Supporting Information.

2.4 Calculation of R_ij and the results

Sections 2.2 and 2.3 show the measurements of dynamic material constants S^D_ij and static material constants S^S_ij and 1/KL^S_ij. The challenge present here is that, though, 1/KL^S_i and S^D_ij are measured simultaneously at hydrostatic pressures, linear compliances 1/KL^S_i can not invert for individual S^S_ij directly. Though observation on the variance on RKL_i can provide some insights into the dynamic-to-static ratios, we still need explore method that can calculate R_ij as a function of hydrostatic pressure σ_H by combining measurements of individual S^S_ij that are measured at different pressures using the tri-axial loading experiments (σ₁₁ may not equal to σ₂₂ or σ₃₃), along with 1/KL^S_i determined directly as functions of the hydrostatic pressures.

Here, after cleaning the data, we first perform parametric multivariate nonlinear regression to predict S^S_ij as function of the hydrostatic pressures. The static and dynamic material constants exhibit different levels of stress dependency with varying sensitivities to other stress tensor components. That being said, generally, samples are expected to be stiffer (or less compliant) at greater pressures and with lower sensitivities to the change of stresses. At higher pressures, the stiffness either becomes semi-constant or changes minimally with the increase in pressure. Ideally, a non-parametric regression might better predict the stress values at different pressures and generally does not come with bias (e.g. Dobson & Barnett 2018). But, due to the lack of sufficient data, we performed multivariate nonlinear regression mainly to avoid data overfitting. Nevertheless, we made different attempts on two different nonlinear method (NLM) fit functions (denoted as NLM and NLM1, see below), where S^S_ij is functions of stress states σ_ii(or simply denoted as σ):

and/or

These two functions have two essential characteristics that meet our expectation on the pressure dependence of material constants: (1) material constants change rapidly against pressure at low pressure and (2) the material constants’ sensitivities against pressures becomes negligible at larger pressures.

All of our measurements are obtained from states where shear stress σ_ij(i≠j) are equal to 0, and our goal is to predict S^S_ij at hydrostatic pressures; thus, only σ_{ii =1,2,3} are considered here. Also, honouring the definition of TI that samples respond the same to Δσ₁₁ and Δσ₂₂, σ₁₁ + σ₂₂ is treated as one variable in our modelling process. The variables b₁, b₂, b₃ and b₄ that describe these two predictive functions are calculated using a multivariate regression with Matlab .

It is not uncommon that the nonlinear parametric fitting returns unreasonable data. For such reason, we compare the results by fitting both eqs (10) or (11), keep the more reasonable ones and discard the results that are apparently wrong. For example, in sample A (Fig. 9), we fit the measured σ_ii and S^S_ij with NLM/NLM1 and plot them alone with the produced prediction; we observe that the results from the fitting of NLM1 are more reasonable than NLM. Notably, in Figs 9(c), (f) and (i), NLM1 predicts results closer to our expectation for S^S₁₁, S^S₁₂ and S^S₁₃ at high hydrostatic pressures.

Figure 9.

Examples of nonlinear fitting from sample A with the measured strain (crosses) and strains predicted by NLM (solid dots)/NML1(diamonds) plotted against the pressures. Data shown in (a) and (b) are used for constraining S^S₁₁, (d) and (e) for S^S₁₂, (g) and (h) for S^S₁₃ and (l) and (m) for S^S₃₃. Two curves showing the predicted S^S_ij by NLM or NLM1 with (c) for S^S₁₁, (f) for S^S₁₂, (i) for S^S₁₃ and (n) for S^S₃₃, against hydrostatic pressures.

Open in new tabDownload slide

We would hope the material constants S^S_ij, extrapolated from the non-linear regression using the input of measurements obtained from the cylindrical plugs, are compatible with the S^D_ij from the prism-shaped samples. However, many of the extrapolated S^S_ij does not honour eqs (4) and (5) if we put KL^S_i, measured from the prism shaped samples simultaneously with the dynamic S^D_ij, on the right-hand sides of these two equations. Some contribution factors include: (1) the different stress paths on which S^S_ij (from plugs with anisotropic stress states) and KL^S_i (from prism-shaped samples with hydrostatic stress states) are measured; (2) uncertainties caused by instrumental errors and linear fitting processes during which slopes of the stress/strain plots are obtained and (3) the heterogeneities of the sample itself. Therefore, the results of extrapolated S^S_ij must be validated, and we do so by comparing extrapolated results with the directly measured 1/KL^S₁ and 1/KL^S₃, before being used for calculation of R_ij.

For such purpose, we perform a Monte-Carlo simulation (n = 2e5) with S^S_ij allowed to vary uniformly and randomly within χ₁ per cent of the reported values and are subject to checks and constraints:

1. We calculate their corresponding expression in E and u convention and they are filtered with regular checks with accepted expectations for similar shaly rocks that are 0.5 > u13 > u12 > u31 > 0 and E3 < E1. This expectation mainly comes from the fact that the rock is essentially the composition of many laminated layers and the rock is more compressible in the direction perpendicular to these lamination planes. The layers themselves, which consist of highly compacted minerals, are stiff and isotropic.

2. We incorporate the prism samples' 1/|${KL}_1^s$| and 1/|${KL}_3^s$| measurements. We calculate 1/|${KL_{1^{\prime}}}^s$| and 1/|${KL_{3^{\prime}}}^s$|, using eqs (4) and  (5) with the predicted static compliances S^S_ij (using NLM1 or NLM2) in each round of the Monte-Carlo simulation, only S^S_ij that produces:

   $$\begin{eqnarray} \left| {{1 {\left/ {\vphantom {1 {KL}}} \right. } {KL}}_1^s - {1 {\left/ {\vphantom {1 {KL_{1^{\prime}}^s\left| / \right|{1 {\left/ {\vphantom {1 {KL_{1^{\prime}}^s}}} \right. } {KL_{1^{\prime}}^s}}}}} \right. } {KL_{1^{\prime}}^s\left| / \right|{1 {\left/ {\vphantom {1 {KL_{1^{\prime}}^s}}} \right. } {KL_{1^{\prime}}^s}}}}} \right| < \chi 2\,\mathrm{ per cent} \end{eqnarray}$$

   (12)

and

are kept.

We then use bootstrapping to calculate R_ij. All accepted Monte-Carlo simulation results of S^S_ij, following the check steps eqs (12) and (13), are divided by S^D_ij for the distributions of R_ij at the given hydrostatic pressures. We fit the R_ij distributions with a normal distribution curve and calculated the most probable values (50th percentile) and the range of uncertainties defined as the 25th and 75th percentile of the fitted cumulative probability function (CDF, Figs 10a–d); also for references, RKL₁ and RKL₃ that are calculated by directly measured KL^S_i and KL^D_i (calculated with S^D_ij) are also plotted (Figs 10e and f).

Figure 10.

Summary of the dynamic-to-static ratios constrained in this work. Each curve represents a different sample on which the dynamic-to-static ratios R_ij are measured. The whisker boxes in (a)–(d) represent the 25th and 75th percentiles of the CDF curves.

Open in new tabDownload slide

χ₁ essentially controls how much our individually measured S^S_ij are allowed to vary randomly. χ₂ controls how many of these randomly varying values will be tolerated when compared with KL^S_i that are known to be directly compatible (as they are measured at the same time) with the S^D_ij. Though χ₁ and χ₂ for each sample are determined separately, χ₁ falls with in the range of 10–25 per cent and  χ₂ is 5–15 per cent; for each sample, we try to reach a balance that χ₁ and χ₂ are reasonably small while we still have enough measurements to calculate R_ij for the pressure ranges of 40–140 MPa that we are most interested at. Details on the selection of χ₁ and χ₂ for each sample are attached in the Supporting Information.

To summarize, we use two nonlinear function, with the expected characteristics of the stress dependency of S^S_ij considered, to predict S^S_ij as functions of hydrostatic pressures through which S^D_ij are measured. We use two KL^S_i, that are measured at the same time as the S^D_ij, to ensure that R_ij are calculated from compatible S^S_ij. At last, due to the technical complexities of the experiments and the first few rounds of trials and errors, some samples have parts of the results that are not usable for analysis. Explanations of the issues with samples B and G are included in the Supporting Information.

3 DISCUSSION

3.1 Making compatible static and dynamic measurements

As discussed before, S^S_ij is measured with the boundary conditions that involve setting components of stresses constant. On the other hand, measuring static C_ij, would require the components of strain to be constant. In practice, it is much easier to control the stress via hydraulic valves but nearly impossible to control how much a sample is allowed to deform. Therefore, we report our dynamic material constants in the convention of S so they can be compared with their static counterparts, which are measured as S_ij natively. Nevertheless, making compatible measurements for the samples’ S^S_ij and S^D_ij, which are pressure dependent and anisotropic, is the biggest challenge in developing this work.

Though still rather scarce in quantities, earlier researchers have reported their efforts in constraining R_ij. However, up to the date for this writing, a few critical issues remained challenging. First, at least one ray of waves travelling at 0° < θ < 90°, ideally at the ranges of 30°–60° for robust measurements, is needed for the calculation of S^D₁₁, S^D₁₂, S^D₁₃ and S^D₃₃. Second, making measurements at high pressures (usually >100 MPa) that can provide meaningful analogues for the in-situ application requires sophisticated apparatus. Third, there is a less well-known ‘beam-skew’ effect (see Fig. 4a) that requires attention: the wavefront, when travelling at an incident angle of 0° < θ < 90° through anisotropic medium, deviates from the wave's indecent direction (Li et al. 2020). Having transducers paired in parallel to this direction can lead to underestimation of the wave velocity. How much the beam skews can not be pre-determined without knowledge of the complete S^D_ij.

Researchers have used square samples (e.g. Takahashi 2011; Jin et al. 2018) with multiple pairs of P- and/or S-wave transducers and deformation measurement sensors (e.g. strain gauge) attached/positioned on each side (Fig. 11a). Though, by measuring samples anisotropic deformation, this setting can provide robust S^S_ij measurements (not including S^S₄₄), this setting can not generate waves travelling in the directions oblique to x_i axes. Therefore, measuring dynamic C^D₁₃ that is needed to be subsequently inverted for S^D_ij, is difficult.

Figure 11.

Schematic illustration of approaches to measure the dynamic and static material constants on the same sample with (a) using square sample. (b) Using an H plug and (c) using a V plug. Solid lines in the sample represent the isotropic plane of TI material.

Open in new tabDownload slide

Another popular approach takes the cylindrical plugs with pairs of transducers placed on the curved wall side or axial flat top/bottom sides of the plug (e.g. Gong et al. 2019; Han et al. 2020). In case of an H plug (Fig. 11b), a typical configuration is having transducers placed diametrically on the radial wall side (Wang 2002a, b): this configuration allows waves to travel through the sample at angles of 0° ≤ θ ≤ 90°. Placing all transducers at the circular wall surfaces of the sample also avoids engineering challenges in applying axial loading pressure through hardened steel cap and generating ultrasonic waves using transducers, that is generally fragile, at the same time in the same direction and also at very high pressures. That being said, this configuration can not account for the beam-skew effect: having the transducers paired diametrically will inevitable lead to the underestimation for P(θ) or S(θ) when 0 <θ < 90.

Further, though possible in theory, static deformations measured through similar step loading procedures described in Shen & Playter (2021, also adopted in this study, see section 2.3) allow the determination of S^S₁₁, S^S₁₂, S^S₁₃ and S^S₃₃. Nevertheless, such measurements rely on the assumption that S^S₁₃ varies negligibly during the step increase of the pressure and S^S₃₃ needs be constrained with the prior knowledge of S^S₁₃ measured at the lower and higher pressure steps. The measurement of S^S₃₃ is limited by how much the instrument's side wall (usually made of rubber) can hold hydraulic radial pressures. In our case, we are limited to wall pressure less than 50 MPa—such limitation is also seen in many other researches.

Thus, even for the sake of obtaining enough static measurements for reasonable extrapolation results, a complementary V plug drilled parallel to the material symmetry axis is necessary to obtain S₃₃ for σ₃ > 50 MPa (Fig. 11c). Static loading tests preformed on V plug can report S^S₃₃ and S^S₁₃. However, the dynamic S^D₃₃ can not be measured at the same sample without measuring P(θ) or S(θ) when 0 <θ < 90 which generally cannot be done on the V plug: similar to the square sample, current technology cannot generate waves travelling through the sample in direction oblique to x₃.

In our work, to reach a pressure that is compatible with the stress magnitudes for the Duvernay Formation that can be as high as 140 MPa (Shen et al. 2019, 2021), we put our sample into a hardened steel-made cell that is filled with hydraulic oil and pressurized up to 180 MPa, this pressure is higher than most other reported measurements. We also achieve the measurements of complete S^D_ij with beam-skew effects properly accounted for. Though we still suffer from the deficiency, as encountered by earlier researchers, that individual S^S_ij are not determined the simultaneously with S^D_ij, we rely on statistical modelling techniques to constrain the expected S^S_ij applicable for the hydrostatic pressures at which dynamic S^D_ij are measured. Uniquely in our work, we utilized the measurements of static linear compliances 1/KL^S₁ and 1/KL^S₃ to provide constraints and validation for the statistically extrapolated S^S_ij.

For the calculation of R_ij, We here chose to (1) convert the dynamic C to S; and (2) extrapolate the static S^S_ij measured at varying stress states for the hydrostatic pressures. We take these steps (with details in Section 2.3) mainly considering that (1) despite the complex wiring, in our experiences, measurements of dynamic C_ij are generally more reliable than our static S^S_ij measurements. Thus, converting dynamic C for S will introduce less error. (2) Static S^S_ij is measured at different stress states and extrapolating them consistently for hydrostatic pressures allow easier and more importantly, consistent comparisons.

Also, by aiming for deriving the S^S_ij as a function of the hydrostatic pressure, we reduce the prediction function to be only dependent on one parameter (hydrostatic pressure σ_H), and greatly improve the reliability of our statistical modelling methods. At last, though the in-situ stress is not hydrostatic but at 50–160 MPa, the sensitivity of material constants to stress variation is relatively low, and even lower for R_ij, which is the main goal of this work. Considering the uncertainties in the reported stress values (Shen et al. 2019,2021), constraining R_ij for hydrostatic pressures at ranges between 50–160 MPa should be able to provide meaning analogues for the rocks’ in-situ constitutive behaviours.

At last, it also needs to be stressed that, though the radial pressures are limited by the apparatus (less then 50 MPa), our samples are tested indeed with directional axial pressures as high as 180 + MPa: for H plug the axial pressure is equal to either σ₁₁ or σ₂₂; for V plug, the axial pressure is equal to σ₃₃. Fig. 12 shows an example that our measurements encompasses a wide range of σ₃₃, σ₁₁ and σ₂₂. Though we do have individual σ_ii high enough that is within the ballpark of the magnitudes of hydrostatic pressures in which dynamic S_ij are determined, our S^S_ij lacks directly measured values from stress environments that all three components of σ_ii are all high. Consequently, taking that S^S_ii is a function of σ_ii, our measurements capture the first-order derivative ΔS_ii/Δσ_jj and second-order derivative Δ²S_ii/Δ ²σ_jj but we have no control of Δ²S_ii/(Δσ_jjΔσ_kk) when j≠k.

Figure 12.

The pressure distribution of strains was measured for each experiment used to determine sample A's S^S_ij. Please note that using either H or V plug σ₂ may be equal to σ₁ or σ₃.

Open in new tabDownload slide

We uses two different regression fitting functions that both have two essential characteristics that meet our expectation on the pressure dependence of material constants: S_ij change rapidly (Δ ²S_ii/( Δσ_jj Δσ_kk) is large) against pressure at low pressure but S_ij’s sensitivities against pressures becomes negligible at larger pressures (Δ²S_ii/(Δσ_jjΔσ_kk) approaches 0)—that is, it flattens when σ_ii becomes greater. These characteristics ensemble the expected characteristics of Δ ²S_ii/(Δσ_jjΔσ_kk). After preforming the regression analysis, we arbitrarily pick the results are more likely to be correct. This step is indeed subjective and thus it is followed by our objective validation steps using KL^S_i: such validation ensures that our extrapolated S^S_ij do not contradict with the direct actual measurements. Further, in the discussion below, when interpreting our results, we look onto the values of KL_i or RKL_i to ensure our interpretation of the data is reasonable.

3.2 Notable caveats

It is worth noting that R₄₄ is missing in our discussion. Though S^D₄₄ is readily reported in dynamic experiments, the determination of S^S₄₄ can only be achieved through a plug cut with off-axis angle 0° < ϕ < 90°. In practise, the lack of constraints on S^S₄₄ is a rather isolated issue that can possibly be mitigated in the future work. First, lack of constrained S^S₄₄ does not impact the calculation of E_H,V, υ_hh,hv,vh, which are more critical in addressing issues for horizontal fracture propagation and borehole stability issues of vertical and horizontal wells; but calculating the stabilities of slanted wellbore does require S^S₄₄. It is possible to explore other experimental (e.g. direct shearing experiment) or empirical constraints (e.g. Saint-Venant's equation) that can estimate S^S₄₄ and R₄₄ to obtain independent measurements without the input of the results reported here.

Another unaccounted factor that can result in discrepancies between S^S_ij and S^D_ij is the stress path (Fjær 2009). We observe the KL^S_i derived from the loading (pressure increase) and unloading (pressure decrease) processes are different (see Fig. 6); the wave velocities at these two processes also show minor differences (see Fig. 5). Such hysteresis is generally well known, and it is suspected to be the results of the friction during the closure of the microcracks (Scholz & Hickman 1983; Nishino et al. 1992; David et al. 2012). Nevertheless, here we take the average KL^S_i measured during the loading and unloading processes. We also calculated RKL_i using the measurements from the loading/unloading process separately and compared them the values calculated using the average KL^D/KL^S for the given pressures: see Fig. 13 for sample C and others are in the Supporting Information. Overall, the uncertainties are higher at low pressures where the measurements are stronger hysteresis are observed. RKL₃ also shows greater uncertainties in the low pressures. This is also understandable as most of crack closure and sliding happens in the direction of x₃ which is perpendicular to the bedding plane.

Figure 13.

RKL_i calculated during the loading (pressurization) and unloading (depressurization) process of the experiments and the variances; figures for other samples can be found in the Supporting Information.

Open in new tabDownload slide

We note that the uncertainties of RKL_i are roughly 10 per cent or less for 40–140 MPa. Though we have no active way to measure the S^S_ij during the stress unloading, as the plugs are generally destroyed in the process aiming to obtain their compressive strength, we believe this error, likely lies within the ball park of 10 per cent, are shadowed by the errors from other static S^S_ij measurements and captured by our Bayesian modelling process that introduce variation of ∼10 per cent (χ₁ ∼=0.1, see Section 2.4 and Supporting Information).

3.3 Controlling factors of R_ij

Subsurface resource development requires knowledge of rock's mechanical properties that are representative for large volumes of the rock mass. However, lab experiments only provide a handful of point measurements—the heterogeneity of rock mass generally prevent direct use of these point measurements for actual engineering work. It is of practical use to explore the potential correlation between different parameters that might be useful for predicting this. Further, bearing in mind that correlation can be mere coincidence and does not necessarily reflect the actual correlation, we here also try to reason whether the observed correlations are plausible from the rock physics perspective.

The complete quantitative measurements on samples discussed here include: (1) geochemical mineral composition; (2) petrophysical properties (i.e. porosity and density) and (3) static and dynamic elastic constants that are used in computing the R_ij; we first calculated the correlation coefficients (Weisstein 2006) between each R_ij and the other measured parameters (Fig. 14); for fair comparisons between different parameters, they are all normalized between 0–1.

Figure. 14

Normalized correlation coefficients between the R_ij and other sample characteristics.

Open in new tabDownload slide

3.3.1 Lack of sensitivity against varying pressures

The most immediate observation is the lack of the hydrostatic pressure σ_H versus R_ij correlation. The R_ij plotted against the pressures (see Fig. 10) confirms this observation. Taking R₁₁ as an example, though samples G and A exhibit some rather weak negative correlation against rising pressures, samples B, C and D remain largely constant across the entire pressure spectrum; similar observations are made with other R_ij where most samples report stable R_ij from 60 to 160 MPa.

Intuitively, this observation contrasts with that S^S_ij and S^D_ij are expected to be different functions of σ. However, we note that both S^S_ij and S^D_ij share similar trends in changing faster against rising pressure at low pressure, and both become rather stable when the pressures are greater than 60 MPa. Previous studies by Gong et al. (2019) experimented on cylindrical plugs and reported varying R₁₁ and R₃₃ against the axial pressures. However, Gong et al. (2019) stopped at 50 MPa axial pressure, and the change of R₁₁ and R₃₃ against pressure is not monotonic but peaked at ∼40 MPa and started declining after. At last, the observation of RKL₁ and RKL₃ shows no discernible pattern the can describe a correlation with the stress (see Figs 10e and f), in agreement with earlier works from Meléndez-Martínez & Schmitt (2016) and Ong et al. (2016) that both report minimal fluctuations of KL^S_i/KL^D_i at pressures higher than 40 MPa.

In our work, we observe S^D_ij (Fig. 5) and S^S_ij (Fig. 8) fluctuating with higher magnitudes at lower pressures (i.e. < 60 MPa) – this suggests that R_ij might also be more sensitive to pressure changes at these pressures. Also, S^D_ij’s and S^S_ij’s trends are similar at greater pressures: when they are competing for influences on R_ij it is also plausible that their impacts cancels out. Nevertheless, our most interested pressures are higher than 50 MPa given the in-situ pressures of the Duvernay reservoir are generally higher. We still believe R_ij will be impacted by pressure, but, considering the practical precision needed for engineering practice, this study shows we can effectively treat R_ij as independent of σ for σ_ii > 50 MPa, that is the case of in-situ pressures of the Duvernay reservoir.

3.3.2 Correlation between R_ij and SD_ij

On top of the lack of sensitivity to the pressure variation, we found that R₁₁ is generally not sensitive to other parameters. It is minorly correlated with S^D₁₂ and S^D₁₃ that are associated with material's compactness (see more discussion in the later text). We plot R₁₁ against S^D₁₁ and 1/KL^D₁ (calculated with dynamic S^D_ij using eq. 4) in Figs 15(a) and (b) and do not observe any meaningful correlational trends. R₁₁ generally ranges closer to 1, with sample G's R₁₁ potentially exceeding R₁₁—though it is debatable as Fig. 10 shows the lower limits for sample G's range of R₁₁ is still less than one and Fig. 15 plots only the 50th percentile of the CDF function that is considered the most probable value but not to be taken as the ‘correct’ value. We expect rocks, mainly, in this case, shale, to be TI, primarily owing to their horizontal lamination. The S₁₁ reflects the material's compressibility on the isotropic plane, and most of the material's void spaces are expected to be between the laminated bedding planes. Thus, it is reasonable to expect the dynamic and static elastic properties of the isotropic bedding planes to be similar.

Figure 15.

Scattered R_ijversus S^D_ij plots with linear fit equations.

Open in new tabDownload slide

On the other hand, the R₃₃ for all samples is significantly smaller than 1, with the lowest observation from sample B (R₃₃ as low as 0.2), and other samples’ report values generally range between 0.5 and 0.8—the differences between S^D₃₃ and S^S₃₃ are rather significant. The difference between rock's dynamic and static stiffness results from several physical phenomena, including strain rate, drainage conditions, heterogeneities, anisotropy and strain amplitudes (Fjær 2019) and potentially more relevant in this work, the void spaces in the sample (Cheng & Johnston 1981). In this work, we used a dry sample and measurements were made on the same model, thus, the impacts of drainage conditions and heterogeneities can be excluded. But, the strain amplitudes and rates are significantly different between our dynamic (1 Mhz frequency and negligible amplitude) and static (loading time of several minutes and amplitudes of a few microns) experiments. For rocks that exhibit TI properties, there are generally more horizontally (i.e. parallel to x₁–x₂ plane) oriented void spaces and that directly impact the materials responses to loadings in the perpendicular direction (i.e. x₃ axis) and causing discrepancies between the dynamic and static elastic constants; the material's response to the loadings parallel to these features is less impacted.

This reasoning is supported by the fact that, unlike R₁₁, R₃₃ is correlated with other dynamic material constants. Particularly, strong correlation between : S^D₁₁ versus R₃₃; S^D₃₃ versus R₃₃; S^D₁₂ versus R₁₂; S^D₁₃ versus R₁₃ and S^D₁₂. These curves are also plotted in Figs 15(c)–(h), and we applied linear regression to these data. First, R₃₃ is positively correlated with S^D₁₁ and S^D₃₃. Noting that a larger R₃₃ means that the difference between S^D₃₃ and S^S₃₃ is small, we suspect such correlations maybe caused by the amount of microfractures in the materials, particularly the horizontal fractures known to be present in TI rocks and causing discrepancies between the dynamic and static rock properties (Cheng & Johnston 1981).

Softer materials (larger S₃₃ and S₁₁) generally have more microfractures; fractures will also result in lower velocities and larger S^D₃₃. However, it should be noted that (1) the impact of fractures on velocities is less significant than the static compliances; and (2) S^D_ij are complex functions of velocities of rays of waves from different directions that may have various sensitivities to the preferably horizontally (parallel to x₁–x₂ plane) fractures. Therefor, rocks with less horizontal microfractures might have S_S and S_D more similar to each other. Though S₁₁ reflects the rock's compressibility in the horizontal direction and shale generally has little void spaces on these laminated bedding planes, S^D₁₁ is also impacted by velocities of waves in different angles and modes—S^D₁₁ is a function of C^D₁₁, C^D₃₃, C^D₁₂ and C^D₁₃ (Lubarda & Chen 2008) that are functions of wave velocities (see appendices) of P(90), P(0) and at least one ray of wave with off axis angle 0<θ<90. Thus, it is acceptable to reason that smaller S^D₁₁ that stands for more stiff materials can be used as a proxy to predict for R₁₁ closer to 1.

The reasoning from the paragraphs above is supported by the observed correlation between R₁₂, R₁₃ and S^D₁₂ or S^D₁₃. S₁₂ and S₁₃ describe the deformation, in response to the directional loadings parallel to the x₁, in the direction of x₂ parallel to the isotropic plane and x₃ perpendicular to that plane. Essentially, more compact and stiffer rocks are expected to have less S₁₃ and greater S₁₂.

We also need to note here again that in the convention described in eq. (1), S₁₂ and S₁₃ are expected to be less than 0: in response to directional compression, the material will expand in the perpendicular directions and expansion is denoted as the negative strain in this work. Therefore, for example, when we describe S₁₂ as greater, it means it is closer to 0; a lesser/smaller S₁₂ means it has a larger absolute value considering the negative sign.

We plot R₁₂ that is observed to be correlated positively with S^D₁₂ (Fig. 15g) and negatively with S_D13 (Fig. 15e). We find that R₁₂ can be greater or less than 1 depending on the samples and pressures. Though the reasoning for relatively less variation of R₁₁ across samples studied here should also apply here, R₁₂ varies much greater between 0.5–3. This is because sample G, which has a large R₁₁ (the only R₁₁ reported to be greater than 1), reports R₁₂ as great as 3—this warrants further investigation. Nevertheless, lesser S_D12 and greater S_D13 mean that materials are more compact, which is a plausible reason why R₁₂ is correlated with S_D12 and S_D13. This line of reasoning can also apply to R₁₃, which is negatively correlated with S_D13: smaller SD₁₃ is an indication that samples are stiffer in the x₃ direction (i.e. less expansion due to Poisson's effect) and have a lesser amount of horizontally shaped void spaces. Additionally, we also observe similar trends from the R₃₃ versus 1/KL^D₁ that describe the material's contract in the x₃ direction in response to the hydrostatic pressures and is equal to S^D₁₁ + S^D₁₂ + S^D₁₃. We do not observe similar for 1/KL^D₃ that is only very weakly correlated with S^D₁₂; we suspect that is because 1/KL^D₃ = 2S^D₁₃ + S^D₃₃: with more compact rocks with S^D₃₃ becomes smaller and S^D₁₃ becomes ‘greater’ that means it approaches 0. The impact of rock's compactness on S^D₃₃ and S^D₁₃ might be negated by each other. At last, we note that the aforementioned observation is in alignment with the RKL₁ versus S^D₁₁ curve (Fig. 15i) that shows only weak, but notably positive correlation; RKL₃ versus S^D₃₃ (Fig. 15j) produce strong, positive linear correlation echoing the observations of S^D₃₃ can be a good indicator for the R_ij, particularly for the material constants that are related to the static response to compressions in the x₃ direction.

At last, we also explored potential correlation between the materials dynamic anisotropy defined by the Thomsen parameters γ, δ and ε (Thomsen 1986) that are essentially functions of the dynamic elastic constants. In Fig. 14, we observe correlation between R₁₂ versus ε and R₁₃ versus γ; they are plotted below in Fig. 16. Parameter γ is the fractional difference between the SH (0) and SH(90); while ε for P(0) and P(90). This observation suggesting that shales with greater anisotropy is likely to have greater R₁₂ and R₁₃, but one might not dive too much into this as our ε and γ are functions of S^D_ij and earlier finding have suggested that R₁₂ and R₁₃ are correlated with S^D₁₂ and/or S^D₁₃.

Figure 16.

(a) Scattered plot for R₁₂ versus Thomsen's parameter ε and (b) scattered plot for R₁₃ versus Thomsen's parameter γ. Please note, for clarity, only measurements made at 40, 80 and 120 MPa are plotted.

Open in new tabDownload slide

In this work, we focus on exploring the correlation between R_ij with S^D_ij mainly because this work is trying to explore ways of converting S^D_ij to S^S_ij through estimated R_ij; the constraints of S^S_ij, in our expected use of the results from this work, will be the end goal of the practice rather than the means used for inferring other parameters. Overall, we found that the dynamic constants of the samples can provide a good indicator for materials’ R_ij. Further, another takeaway message from R_ijversus KL^D_i plots is that future study might continue to explore functions of different combinations of S^D_ij as proxies to predict R_ij. This finding is important as, in practice, obtaining dynamic constants for larger volumes of the crust is often much more accessible and our findings can help practitioners to incorporate measurements from various sources to constrain geological materials’ static elastic constants applicable for reservoir scale operations.

3.3.3 Impact of petrophysical properties

Following our observation that indicates the compactness of the rock likely controls the R_ij, we explore the correlation between R_ij and the rock's petrophysical parameters. In Fig. 14, we observe that R₃₃ is highly correlated with samples’ densities and porosities; Figs 17(a) and (b) show the correlation between the bulk density or porosity against the R₃₃. We first observe that the density (Fig. 17a) is inversely proportional to the R₃₃; this seems counterintuitive as we generally expect materials with larger density to be stiffer and more compact. Further, Fig. 17(b) shows that the R₃₃ is also inversely proportional to the porosity. Though one might not read too deeply into this as the measurements from sample B (a possibly outlier) largely dictate the linear correlation. Nevertheless, this observation agrees with the expected results that more void spaces can result in greater discrepancies between the dynamic and static properties of the samples (Cheng & Johnston 1981). That being said, considering that our sample's measured porosities have low variation and are small (all less than 1 per cent with the exception of sample B at 1.2 per cent, still a small value), it is also understandable that the other four samples with porosities way below 1 per cent does not show correlation between R_ij. Similarly, we do not observe linear correlation between RKL₃ and density; only very weak negative correlation is observed for Porosity and RKL₃ (Fig. 17e). Please note, that due to the technical failure (see Supporting Information), sample B does not report static moduli KL^S₃ and thus we do not have R₃₃ reported. Also, for the same reason, we do not observe correlation between samples’ porosities and densities either—it is likely that for the rocks tested here, the variation of density and porosity is so small and they have little influence on the rock's R_ij.

Figure 17.

R₃₃ versus (a) samples bulk density, (b) porosity and (c) total organic carbon. The black line denotes the fitted linear function; and (d) and (e) show the same plots but for RKL₃. Please note, for clarity, only measurements made at 40, 80 and 120 MPa are plotted.

Open in new tabDownload slide

3.3.4 Impact by the mineral composition

Earlier workers have attempted or proposed potential correlations between minerals and rock's mechanical properties that might be useful for inferring for various rock mechanical parameters, however, Fig. 14 shows that the R_ij’s correlation with the mineral contents is low. At the first glance, this observation is counterintuitive as different minerals will have varying physical properties and they composition should change the properties of the rock. Following the philosophy of the effective medium theorem (Choy 2015), the bulk rock's mechanical properties might be an weighted average of that of composite minerals (Wang et al. 2001). Here, for the minerals that are present in our samples, the bulk stiffness of these minerals follows Dolomite > Calcite > Muscovite > Quartz; Pyrite's bulk stiffness is much greater than other minerals but only consists of a very small percentage of the total mass and does not contribute much. If a simple weighted average of each mineral's bulk stiffness are added, we would expect the compliances of our samples to have the order of C < B ∼ A < D ∼G (Katahara 1996; Prasad et al. 2002; Wang et al. 1998) that is different from our measured KL^S₁ (G < A < B < C < D), KL^S₃ (C∼A < D < G), KL^D₁ (B < A ∼C < D∼G) and KL^D₃ (B < C < A < D < G). The authors have argued that the for rocks that are essentially packed granular particles, many other parameters that describe particles cohesion, contact fracture and laminated fraction, etc., all impacts the composite rock's mechanical properties (Allo 2019). After all, we cannot establish that mineral compositions can be used as a useful proxy for R_ij through the observations made in this work alone.

Though we do not manage to establish a meaningful correlation between R_ij and mineral compositions, weak correlation between R₃₃ and TOC is indeed observed (see Fig. 17c); and similar observation can be made on the RKL₃ (Fig. 17f). Again, sample B does not have KL^S₃ and is excluded in this plot. Nevertheless, we do not have a reasonable explanation for how TOC controls R₃₃, and correlation does not necessarily mean causation. Regardless, such reporting is still meaningful as an exploitable HF reservoir generally requires TOC to be greater than 5 per cent, and in our case, R₃₃ is between 0.5–0.8 with TOC between 5 and 9 per cent.

At last, we acknowledge that it is more plausible that each R_ij sensitive multiple parameters and it will be more scientific to preform a more proper analysis with R_ij as a function of more than one parameters; however, with the number of samples available for analysis here (five samples), we will risk overfitting the data and produce R_ij function that does not exist in actual case.

3.4 Implication for field applications

The major motivation for this work is to provide a series of lab measurements that can allow proper conversion of more easily attainable dynamic material constants to the static material constants for larger volume rock masses of the Duvernay unconventional reservoir. Following the conventional approach that is adopted by many earlier studies and to ensure compatibilities with their results, we dry samples and ensure the removal of pore fluids. Reasoning by (Dong et al. 2020), through a review of several literatures, suggest that at temperature lower than 140 °C, the impact of the dehydration on their samples’ property is limited and the organic content of the samples investigated here are lower.

Thus, when doing dynamic-to-static conversions for the material constants that are directly calculated from the in-situ measurements (e.g. sonic logs and seismic surveys), fluid substitution must be carried out. At last, considering, (1) our sample is not very porous (with void spaces generally less than 1 per cent) and (2) the Duvernay shale are known for its low permeability that generally does not allow fluid movement. It is arguable that the rock's property measured here do not differ significantly from the in-situ environment. This reasoning is echoed by earlier works from (Ong et al. 2016) where it is found the rock samples’ measured wave velocities are similar to the measurements calculated from the in-situ borehole sonic logs.

4 CONCLUSION

This paper presented our efforts in constraining the dynamic-to-static ratios, with pressure dependency and anisotropy fully considered, of several samples extracted from the Duvernay unconventional reservoir, a prominent shale play in western Canada that is poised to be a major oil/gas production hotbed.

First, we described a unique approach that determines all the required measurements within a piece of material to mitigate the impact of heterogeneity with complete account for anisotropy. Our approach also provides the measurements for pressures as high as 160 MPa that meets the requirement for exploitation of the Duvernay unconventional reservoir.

Using these methods, we measure five samples that can serve as reference points, with their complete dynamic, static elastic constants and geomechanical/petrophysical characteristics constrained, for the reservoir engineer's use in designing the exploitation strategy for the Duvernay formation. Further analysis of these measurements reveals that: (1) for pressures greater than 50 MPa, R_ij is generally not sensitive to pressure variation; (2) R₁₁ appears to be not impacted by the varying other parameters and generally falls within the range between 0.8–1 and (3) R_i3 and R₁₂ are correlated with the materials’ dynamic elastic constants and dynamic anisotropy (defined by ε, γ), we reasoned this is related to the compactness rock and the horizontal void spaces of the TI sedimentary rocks. (4) For samples studied here, we do not find a correlation between the R_ij and their mineral contents. On top of these findings, the complete data acquired through our work, including the analysis results, are available for access to readers who intend to perform further investigation or more generalized meta-analysis.

At last, we stress that though this work generally focuses on the shale from the Duvernay unconventional reservoir, the workflow described here can be further adapted and potentially improved, with a few improvements suggested in this writing already, for other works that requires accurate determination of a rock's elastic mechanical properties at large.

ACKNOWLEDGEMENTS

LS's contribution is supported by NSFC (grant no. 42227805) and Fundamental Research Funds for the Central Universities, China University of Geosciences (Wuhan) (grant no. CUG240619). The experiment components of this work is supported by Canada NSERC Discovery Grant program to DRS.

AUTHOR CONTRIBUTIONS

Luyi Shen: conducting experiments, formal analysis; conceptualization, methodology, and writing; Douglas Schmitt: conceptualization; funding acquisition, formal analysis, writing-review; Yu-yong Jiao: funding acquisition

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

DATA AVAILABILITY

The experiment data is available at: https://doi.org/10.5683/SP3/CV6WOX

REFERENCES