Empirical rock physics relationships on carbonate dry-frame elastic properties

Abstract

Laboratory tests were conducted to analyze the ultrasonic velocity response to the pressure change in dry carbonate rocks from the Weyburn oilfield, Canada. Twenty-four samples are from seven wells with helium porosities ranging from 1% to 29%. Thin-section images, SEM and mercury intrusion porosimetry were performed to show their inner structures and pore throat size distributions. P- and S-wave velocities (V_p and V_s) measurements were first done under hydrostatic loading and then while unloading, with confining pressures varying between 3 and 35 MPa. The results indicate that V_p and V_s in these samples follow a linear relation independent of the pressure change. The ratio V_p/V_s is more responsive to pressure change irrespective of the pore volume. One-third of the carbonate samples show abnormal V_p/V_s reduction with the increase in the effective pressure. The pressure dependence of velocities (PDV) of Weyburn carbonate rocks varies widely even for samples from the same formation with similar sedimentary history. Samples with loosely packed crystals and/or relatively large dominant pore diameter have higher PDV. The exponential empirical model V = $A-C\text{e}^{DP_{\text{e}}}$ was tested; therein, V is the elastic wave velocity, P_e is the effective confining pressure, and A, C and D are the best fitting coefficients determined by curve fitting. The model gives good fits for most of the Weyburn carbonate samples. From a statistical point of view, there is no difference between the V_p- and V_s-derived exponential coefficient D.

1 Introduction

The change in stress condition influences the elastic frame moduli of naturally occurring rocks and accordingly results in the change in their elastic wave velocities. The pressure dependence of velocities (PDV) of upper crust rocks is of vital importance for a variety of applications, such as pore pressure investigation of hydrocarbon reservoirs (Xu et al. 2006) and time-lapse monitoring of underground fluid production or injection (Pervukhina et al. 2010; Perozzi et al. 2017a, b; Zeng et al. 2019). Nonlinear response of elastic wave velocities to the effective pressure change has been observed for a long time in a wide variety of rocks, such as mudstones and sandstones (e.g., Han et al. 1986; Freund 1992; Khaksar et al. 1999; Sayers 1999, 2002; Ong et al. 2016; Gao et al. 2019; Li et al. 2019a). From these studies, it is clear that with increasing effective pressure, the elastic wave velocities of rocks first quickly increase and then get slower before reaching a plateau (Fig. 1). Phenomenologically, this is explained by the progressive closure with confining pressure of small aspect ratio and compressible crack-like pores (Walsh 1965; Batzle et al. 1980). For hypothetical elliptical cracks, those with the smallest aspect ratios will close under the lowest pressures. Upon each closure, the material becomes incrementally stiffer and more rigid and its elastic wave speeds correspondingly rise (e.g., the numerical modeling work by Wang et al. 2016). The rate of change depends on the distribution of the crack-like porosity. Ideally, all these cracks are closed once a certain confining pressure is obtained, whereupon further compression results in the almost linear stress–strain relationships for rocks.

Fig. 1

The nonlinear pressure dependence of velocity of porous rocks, and the curve fittings based on Eqs. (1), (2) and (4). The measurements shown here are from a sandstone sample (sample P47, porosity = 0.1437, clay content = 0.13) in Han et al (1986)

For operating convenience, a number of simplified mathematical expressions have been used to provide a descriptive relationship between the wave speeds and the confining pressure. These usually assume a trial equation whose coefficients are calculated in linear or nonlinear least squares regression of the observed velocities versus pressure. The most popular one is a combination of exponential and linear items based on velocity measurements (Eberhart-Phillips et al. 1989):

$$V = A + B\sigma - C{\text{e}}^{{ - DP_{{\text{e}}} }}$$

(1)

where V is the elastic wave velocity, P_e is the effective confining pressure and A, B, C and D are the best fitting coefficients determined by curve fitting. From Eq. (1), the linear item can be ignored at relative low pressure (P_e < 100 MPa) according to the analysis of Shapiro (2003), resulting in:

$$V = A - C{\text{e}}^{{ - DP_{{\text{e}}} }}.$$

(2)

This empirical model gives satisfactory approximations to the measurements of sandstones. Another type of empirical model describing the nonlinearity of rocks comes from the strain–stress measurements (Zimmerman et al. 1986), such as (Melendez-Martinez and Schmitt 2013):

$$\varepsilon = a + b\sigma + c\sigma^{0.5}$$

(3)

where a, b and c are determined by curve fitting. This empirical model provided a superior fit of linear strains with pressure in low-permeability unconventional reservoir shales. For the analysis of PDV, Eq. (3) can be converted to the velocity form with similar variation trend:

$$V = a + bP_{{\text{e}}} + cP_{{\text{e}}}^{0.5}.$$

(4)

Figure 1 shows the fittings of Eqs. (1), (2) and (4) based on the velocities of a sandstone sample, with all giving good performances.

Relative to siliciclastics, carbonate rocks have a much more heterogeneous pore network (Xia et al. 2019), which results in larger variances of the acoustic properties for a given porosity (e.g., Assefa et al. 2003; Kenter et al. 2007; Verwer et al. 2008; Weger et al. 2009; Wang et al. 2015; Li et al. 2019b). First, in carbonate rocks the velocity–porosity relation shows wide scatter (Anselmetti and Eberli 1993; Njiekak and Schmitt 2019; Li et al. 2020). This is caused among others by the composition of the mineral grains, the texture of the rock frame and the pore characteristics. Second, the PDV of carbonate rocks still lack consistent understanding (Wang and Wang 2015). Generally, the velocity–effective pressure curve of carbonate rocks is very flat; the observed hysteresis effect following unloading in carbonates underground indicates low PDV (Huffman 2002). Gomez et al. (2007) analyzed the pore structure and the PDV of 50 carbonate samples. The authors found that the PDV of samples dominated by interparticle/intercrystalline pores is stronger than those of samples displaying mostly stiff pores, such as moldic pores, intraparticle pores and vugs. The observations of Guo et al. (2012) on carbonate outcrops from Sichuan Basin, China, indicated that the PDV of dolomites with crack structures is stronger than that of limestones.

As for the Weyburn reservoir formations investigated in our study, Njiekak and Schmitt (2019) examined the influences of pore characteristics and pore fluids on the stress dependence of velocities on a small set of carbonate samples (two limestones and one dolomite). The study was developed from an earlier extensive laboratory testing of 24 Weyburn samples (Schmitt et al. 2012). On the studied subset of samples, the authors also tested the empirical model of Todd and Simmons (1972) for effective stress coefficient for seismic velocity. Moreover, Melendez-Martinez and Schmitt (2013) investigated the elastic anisotropy change with confining pressure on five samples from the same Weyburn reservoir.

Here, we focus on PDV under dry conditions on the larger set of the Weyburn reservoir samples including 22 carbonate samples and 2 sandstone samples. We will first provide an overview of the geological background and textural features of the 24 studied Weyburn samples. Results of P- and S-wave velocities acquired at different confining pressures will then be given. Using these measurements, the PDV of the dry carbonate rocks will be analyzed in three different ways in order to seek useful empirical trends that may be applied in estimating rock properties under differing conditions: the V_p–V_s linear relation, the change in V_p/V_s ratio and the normalized velocity. Empirical velocity–pressure relations (Eqs. (2) and (4)) will then be tested on the studied carbonate rocks. Finally, the potential controlling factors on PDV for carbonate rocks will be discussed. Overall, this study aims to provide a new set of experimentally measured P- and S-wave speeds in a suite of carbonates, to assess the PDV of these carbonates and the factors controlling this and finally to test some existing empirical velocity–pressure relations that are often employed in the literature.

2 Geological background

The carbonate samples used in our study are from the Weyburn oilfield in western Canada as shown in Fig. 2a. The burial depth of the pay-zone is about 1.5 km within Mississippian carbonates in the Williston Basin (e.g., Burrowes 2001; Whittaker and Rostron 2003; Whittaker et al. 2004). The upper carbonate layer, referred to as the Midale Marly, is largely composed of dolostones. Underneath the Midale Marly is the Midale Vuggy, which consists of limestones. Evaporite formations (i.e., the Midale and the Frobisher evaporites) contribute to the upper and lower seals of the reservoir (Fig. 2b). Above the Weyburn carbonates and evaporite formation is the relatively impermeable Triassic Lower Watrous formation, which is a regional sub-Mesozoic unconformity (e.g., Whittaker et al. 2004).

Fig. 2

modified from Proceedings of the 6th International Conference on Greenhouse Gas Control Technologies, Volume I, Whittaker and Rostron (2003), Geologic storage of CO₂ in a carbonate reservoir within the Williston Basin, Canada: an update, pp. 387, Copyright (2003), with permission from Elsevier

The location of Weyburn oilfield and its stratigraphic section. a The location of Weyburn oilfield, from Underground Gas Storage: Worldwide Experiences and Future Development in the UK and Europe, Volume 313, Riding and Rochelle (2009), Subsurface characterization and geological monitoring of the CO₂ injection operation at Weyburn, Saskatchewan, Canada, pp. 229, with permission granted according to the Geological Society of London fair use policy; b the stratigraphic section of the Weyburn oilfield,

3 Sample characterization

In our study, we are using 24 geologically representative core samples (diameter = 38.1 mm) of the Weyburn oilfield, SE Saskatchewan, collected at depths of about 1313.71 m to 1448.26 m in seven different wells (Fig. 3). Among the 24 studied samples, there are two sandstone samples from the Watrous Red Bed unit, 11 dolomite samples belonging to the Marly, Midale evaporite and Frobisher evaporite units and 11 limestone samples from the Marly and Vuggy units. Four dolomite and one limestone cores were each further cut into three sub-cores (diameter = 25.4 mm, length = 20.1−37.0 mm) normal, oblique, or parallel to the bedding plane. The sub-cores were used for anisotropy studies (see study by Melendez-Martinez and Schmitt 2013). Here, the subscripts V, O and H are attached to their indexes, which stand for vertical, oblique and horizontal, respectively. For example, sample 13-O was cut obliquely to the rock beds; i.e., the main axis of the tested cylindrical plug is oblique to the bedding plane. When no letter is attached to a sample index, it means that the axis of the core sample is parallel to the vertical, i.e., normal to the bedding plane.

Fig. 3

Core pictures of the Weyburn field from seven wells, a from Well 10-11-6-14W2, b from Well 2-16-6-13W2, c from Well 6-1-6-13W2, d from Well 2-10-6-14W2, e from Well 8-36-6-14W2, f from Well 4-23-7-14W2 and g from Well 4-16-7-13W2

The mineral compositions of the samples were determined through XRD and XRF whole rock analysis. Thin-section images and scanning electron microscopy (SEM) of the samples were performed to show their inner structures. Pore throat size distributions were assessed by mercury porosimetry. A He pycnometer was used to determine the grain densities and the rock porosity. The air permeability as measured at room temperature of most of the samples is lower than 2mD. Their porosity from He pycnometer ranges from 1% to 29% with a mean value of about 11%. Further details of these measurements are available in Njiekak et al (2018). In the following, sample characteristics will be described for each formation unit. A summary of the sample characteristics is given in Table 1. The data and characterizations provided here all originate from the earlier reporting by Schmitt et al. (2012).

Table 1 Summary of the sample characteristics for each formation unit. In the third column from the left, the main minerals are in bold character. Other minerals, whenever present, form less than 5% of the rock solid frame

3.1 Watrous red bed unit

The Watrous Red Bed unit is represented by samples 39 and 43. They have a porosity of 7% to 11% and a grain density varying from 2.74 g/cm³ to 2.76 g/cm³, respectively (see Table 4 for the supporting data in Appendix 1). Porosity is mainly intercrystalline and the grain size is up to 100 µm as depicted in Fig. 4a. Quartz and feldspar are the main minerals. Anhydrite, dolomite, calcite and clay-size crystals usually fill the space between the quartz and feldspar grains (Fig. 5a). Low content in anhydritic and carbonate cement is usually indicative of high pore volume and permeability (up to 26 mD) as in sample 43. The dominant pore throat diameter is around 4 µm (Fig. 6a).

Fig. 4

Thin-section images of Weyburn samples from the different units. They were taken under plane-polarized light except images (c) and (d). The latter were taken under crossed Nichols. a Sample 43 from the Watrous Red Bed unit; b sample 29 from the Top Midale Evaporite unit. Discontinuous layers of anhydrite (see arrow) are cutting through its dolomite-rich matrix; c sample 20 from the Marly unit. Note the presence of vugs and channels (see arrows); d sample 31 from the Marly unit. The photomicrograph illustrates anhydrite cement filling pore space; e sample 7 from the Vuggy unit displaying a lime-mud-dominated fabric. Arrows indicate vugs. Vugs of much more larger size are common in some samples from the Vuggy unit as illustrated in Njiekak et al (2018; e.g., their Fig. 3a); f sample 38 from the Frobisher Evaporite unit. Note the dolomite rhombs disseminated in the micritic matrix

Fig. 5

SEM images of representative samples for different units. a Sample 43 from the Watrous Red Bed unit. Image showing anhydrite filling spaces between quartz and feldspar grains; b sample 29 of the Top Midale Evaporite unit. Tightly packed crystals forming the rock solid frame; c sample 20 from the Marly unit. Its frame is characterized by loose arrangement of the dolomite crystals; d sample 31 from the Marly unit. Note the tight arrangement of the dolomite crystals; e sample 7 from the Vuggy unit. The arrow indicates a vug being filled by cement crystals (calcite); f sample 11 and g sample 12 from the Vuggy unit displaying vugs formed following the removal of fossils (likely through dissolution); h sample 38. Its rock matrix is made of relatively tightly packed crystals

Fig. 6

Pore throat distribution obtained through mercury porosimetry for the different units. a Sample 43 from the Watrous Red Bed unit; b sample 45 and c sample 13 from the Top Midale Evaporite unit; d sample 31 and e sample 20 from the Marly unit; f sample 12 and g sample 7 from the Vuggy unit; and h sample 38 from the Frobisher Evaporite unit

3.2 Top midale evaporite unit

Three samples from the Top Midale Evaporite unit were investigated. Overall, these samples have a tightly packed and dolomite-rich matrix, which is crosscut by discontinuous thin anhydritic layers (Figs. 4b and 5b). Their grain density ranges between 2.73 g/cm³ and 2.89 g/cm³ (Table 4), and their porosity, which is intercrystalline, ranges from 1 to 3%. The main pore throat diameter ranges between 0.015 µm and 0.7 µm (Fig. 6b, c). Permeability is lower than 2 mD in the three samples.

3.3 Marly unit

There are ten samples from the Marly units, including eight dolostone and two limestone samples. The grain density varies between 2.64 g/cm³ and 2.83 g/cm³, and the permeability is also smaller than 2 mD (Table 4). Besides some locally distributed vugs or channels (e.g., Fig. 4c), there are mainly intercrystalline pores (Fig. 4d). The main pore throat diameter in those samples varies from 0.02 µm to 1.5 µm (Fig. 6d, e). Based on their porosity, they can be subdivided into two groups:

• 1 The tight Marly samples (Fig. 5d) with porosity up to 8% and composed of samples 9, 16 and 31. Their tight matrixes occasionally have sparse channels and vugs, which are often filled with anhydrite (Fig. 4d).

• 2 The loose Marly samples with porosity ranging from 17% to 27%, such as samples 20, 26 and 46. Their frames commonly display alternating loose domains (Fig. 5c).

3.4 Vuggy unit

The Vuggy unit is represented by eight limestone samples (Table 4). They are mainly wackstones and packstones, except sample 42 which is a grainstone (Njiekak and Schmitt 2019; their Fig. 1e, f and g). Porosity in those samples ranges from 4% to 16% and may be vuggy (Fig. 5f, g), moldic, intergranular and intragranular. Their permeability is generally good, ranging from 16 to 50 mD (Table 4). Recrystallization usually leads to the reduction of porosity by filling the pore space (Fig. 5e). The dominant pore throat diameters in the Vuggy samples vary from 0.01 µm to 8 µm (Fig. 6f, g), and some can reach 35 µm (Njiekak and Schmitt 2019; see sample OL in their Fig. 2, sample OL = sample 42 in this paper).

3.5 Frobisher Evaporite unit

Sample 38 is from a lime-rich interval within the Frobisher Evaporite unit. It has a tight matrix mostly made of calcite and dolomite crystals (Figs. 4f and 5h). Other minerals include quartz and feldspar.

The grain density of sample 38 is 2.73 g/cm³. Its porosity is 9% and is intercrystalline, with the dominant pore throat size around 0.1 µm (Fig. 6h). Its permeability is lower than 2 mD (Table 4).

4 Velocity responses to pressure change

4.1 Velocity measurements

Experimental setup (Fig. 7) and procedures for the P- and S-wave velocity (V_p and V_s) determinations were given in Njiekak et al (2013). Briefly, the tested samples were all cylindrical plugs of 38.1 or 25.4 mm in diameter and of about 20.1 to 58.0 mm in length (Table 4). Once placed inside the pressure vessel, the sample’s pore space was subject to vacuum for about 8 h (pore fluid pressure, P_p = 0 MPa). Then, the hydraulic oil pressure within the vessel (confining pressure, P_c) increased from P_c = 3 MPa to P_c = 35 MPa, stepped by 2 or 3 MPa. During this hydrostatic loading process, the V_p and V_s of the dry sample at each pressure point were measured with the standard pulse transmission at a central frequency equal to 1 MHz. The main source of errors in the velocity measurements is associated with the picking of travel times from the acquired raw waveforms (e.g., Njiekak et al 2013, their Sect. 4.4). Overall, when considering all the sources of uncertainties in our measurements, the percentage errors are up to 0.7% and 0.4% for P- and S-wave velocities, respectively. The velocity measurements were repeated at the same pressure points while unloading. For these measurements done under vacuum dry and hydrostatic conditions, the Terzaghi effective pressure is the same as the confining pressure.

Fig. 7

Schematics of the pressure vessel and the data acquisition system for the P- and S-wave velocities measurements. Two different pump systems controlled the confining and the pore pressures. For the measurements in this paper, the pore pressure is zero

The measurements show that the velocities nonlinearly decrease with the increase in porosity, but it is hard to use one simple relation to fit the wide scatterings (Fig. 8). Such scatter in the velocity–porosity cross-plot for carbonate rocks is consistent with observations made in previous studies (e.g., Anselmetti and Eberli 1993; Assefa et al. 2003; Rogen et al. 2005; Verwer et al. 2008; Weger et al. 2009). Also, a broadly accepted rationalization is that the complex textures and related pore structures of carbonate rocks are responsible for this wide dispersion (Weger et al. 2009; Wang et al. 2015).

Fig. 8

The V_p (circle points) and V_s (triangle points) at P_e = 35 MPa of the dry Weyburn samples whose main axis is normal to the bedding plane. The blue color represents samples whose dominant mineral is dolomite, while red color is for calcite-rich samples. The solid and dashed curves are obtained using differential effective medium concepts (Berryman 1992) with pore aspect ratios of 0.2 and 0.05, respectively. For these calculations, the rock matrix is assumed to be pure calcite (red curves) and dolomite (blue curves), respectively. The arrows indicate the velocities of three samples from the Top Midale Evaporite (samples 13, 29, 45). Sample 45 that contains noticeable amount of quartz and anhydrite according to the XRD data

For comparative purposes, we are using a differential effective medium (DEM) model (Berryman 1992; see detail in Appendix 2) to calculate the velocities of hypothetical rock frames composed entirely of a single carbonate mineral but with different pore aspect ratios (ARs). For our calculations, the rock matrix is assumed to be either pure calcite or dolomite. The values of the input parameters are shown in Table 2. From the DEM modeling, the velocity–porosity relation depends on the pore aspect ratio. The latter defines the stiffness of the rock skeleton. The velocity–porosity relation is nonlinear, especially at low ARs (AR < 0.1). For the same porosity, samples with lower effective AR have lower velocities. Irrespective of the dominant mineral composition, the measured velocities overall fall within the velocity envelopes corresponding to AR = 0.05 and AR = 0.2 (Fig. 8). This means the effective ARs representing the stiffness of the dry skeleton of these samples vary between 0.05 and 0.2, pointing out significant differences in their pressure sensitivity. Sample 45 is dolomite-dominated but, unlike samples 13 and 29 from the same unit and also dolomite-rich, it falls out of the velocity envelopes representing AR = 0.05 and AR = 0.2. A possible explanation for this lies in the presence of clay-like crystals in sample 45 (Fig. 9). Those clay-like crystals, which likely originated from feldspars breakdown, led to lower velocities.

Table 2 Parameters used for DEM model calculation*

Fig. 9

SEM image of sample 45. Note the clay-like crystals wrapping around the other grains (see red arrows). The clay-like crystals are not visible in other samples from the same unit (e.g., sample 29 in Fig. 4b). They may have formed here from the breakdown of feldspars present in the rock

The measured V_p during the hydrostatic loading process on the 24 samples is shown in Fig. 10 (see detailed results in Table 5). Their formation units and the samples' dominant minerals are also reported in the same figure. The measurements during the unloading process are not shown here. In fact, the velocities measured during the unloading process were usually higher than those recorded during the loading process. This is typical for such experiments when there is not enough time left between measurements during the unloading cycle. The observed discrepancy results from the fact that the pace at which pores close during the loading is faster than the pace at which they open during the unloading cycle. Indeed, with enough time left between measurements during the unloading cycle, all the measured velocities will match the “loading” velocities unless the tested sample at some point experiences mechanical damage. We can confirm here that no mechanical alteration of the samples occurred. In fact, the same loading process was applied to a subset of the tested samples at a later stage and similar velocities were obtained.

Fig. 10

The measured P-wave velocity for the 24 dry Weyburn samples during the loading process. The samples are grouped according to their formation units in the graph horizontal axis. For each column, it shows the velocity increases for each sample during the loading process (confining pressure goes up from 3 to 35 MPa). The dominant mineral in each sample is indicated by the color bar above the sample index

As a whole, the measured velocities increase with the effective confining pressure, but the PDV are quite different from each other. Wide variations in PDV also exist among samples of the same formation, i.e., samples originating from a similar depositional environment (Fig. 10). All this makes it difficult to analyze the effects of the mineral composition or the rock texture on the PDV by solely relying on the velocity–effective pressure plot.

4.2 Pressure dependence of velocities

For the pressure dependence of velocities (PDV) analysis, we choose as maximum and minimum effective pressures, 35 and 6 MPa, respectively. This is because, for most of the studied carbonate rocks, high-quality velocity data were obtained from P_c = 6 MPa upward and the P- and S-wave velocities barely changed when the effective pressure was higher than 35 MPa. For consistency, in the following sections only the “vertical” cylindrical plugs are considered in our analyses, unless otherwise stated. For example, sample 13 and sample 13-V are all vertical samples collected from the same core material. However, the former was considered together with vertical samples from other cores in the statistical treatment of our data. Moreover, velocities recorded on samples 13-V, 13-O and 13-H are quite different from those measured on sample 13 (Fig. 10). The discrepancy is likely related to textural differences between samples 13-V, 13-O and 13-H (all cut from a single big rock fragment) and sample 13. In fact, the core from which the four samples were obtained displays strong local textural variations.

4.2.1 V _p –V _s linear relation

Well-known V_p-V_s relations for carbonate rocks were developed in Castagna et al. (1993), wherein the empirical relation for dolomites is linear while that for limestones requires a second polynomial. The measurements of the studied 22 carbonate samples from the Weyburn oilfield are slightly different from those reported in Castagna et al. (1993). Our newly obtained carbonates dataset, regardless of the rock mineralogy, follows a linear V_p-V_s relation (Fig. 11):

$$\left\{ \begin{gathered} V_{{\text{s}}} = 0.4312V_{{\text{p}}} + 0.7134\;\;\;\;\left( {{{{\text{km}}} \mathord{\left/ {\vphantom {{{\text{km}}} {\text{s}}}} \right. \kern-\nulldelimiterspace} {\text{s}}}} \right),\;P_{{\text{c}}} \,= \,6{\text{MPa}} \hfill \\ V_{{\text{s}}} = 0.4101V_{{\text{p}}} + 0.8177\;\;\;\;\left( {{{{\text{km}}} \mathord{\left/ {\vphantom {{{\text{km}}} {\text{s}}}} \right. \kern-\nulldelimiterspace} {\text{s}}}} \right),\;P_{{\text{c}}}\, = \,35{\text{MPa}}. \hfill \\ \end{gathered} \right.$$

(5)

Fig. 11

The linear relations of P- and S-wave velocities recorded on the dry carbonate specimens at P_e = 6 MPa (red color) and P_e = 35 MPa (blue color), respectively. The circle points represent the samples whose dominant mineral is dolomite, while the triangle points represent calcite-rich samples. The two fitting lines nearly overlap, which indicates that the linear relation of V_p and V_s is independent of the pressure condition

Also, linear fittings of the V_p-V_s relation corresponding to different pressure conditions (e.g., at P_e = 6 and 35 MPa; Fig. 11) almost match. According to the assumption of Shapiro (2003), for a given rock, the fitting coefficient D in Eqs. (1) and (2) should be in the first approximation the same for both V_p and V_s, that is:

$$\left\{ \begin{gathered} V_{{\text{p}}} \,=\,A_{{\text{p}}} - C_{{\text{p}}} {\text{e}}^{{ - DP_{{\text{e}}} }} \hfill \\ V_{{\text{s}}} \,= \,A_{{\text{s}}} - C_{{\text{s}}} {\text{e}}^{{ - DP_{{\text{e}}} }} \hfill \\ \end{gathered} \right.$$

(6)

which can be further rewritten as follows:

$$\frac{{V_{{\text{p}}} - A_{{\text{p}}} }}{{V_{{\text{s}}} - A_{{\text{s}}} }} = \frac{{C_{{\text{p}}} }}{{C_{{\text{s}}} }}.$$

(7)

It should be noted that fitting parameters in Eq. (7) are independent of the pressure, which means that this linear relation between P- and S-wave velocities is suitable for different pressure conditions. As a result, the same empirical model for S-wave velocity prediction can be used for certain formation even though the formation depths are quite different.

4.2.2 The change in V _p –V _s ratio

The V_p/V_s is a commonly used attribute for fluid and lithology identification by geophysicists. V_p/V_s variation with effective pressure is also suggested as an indicator of overpressure and shallow water flow (SWF) zones (Dvorkin et al. 1999; Prasad 2002). For sandstones, V_p/V_s decreases with the increase in effective pressure under liquid-saturated conditions, and the opposite is true under dry or gas-saturated conditions (Dvorkin et al. 1999; Dutta 2002). Testing the suitability of this rule in carbonate formations needs to be done with caution because of the complexity of their textures.

In order to quantitatively describe the V_p/V_s variation with effective pressure, the change in V_p–V_s ratio (Δ(V_p/V_s)) in this paper is defined as:

$$\Delta \left( {\frac{{V_{{\text{p}}} }}{{V_{{\text{s}}} }}} \right) = \left. {\frac{{V_{{\text{p}}} }}{{V_{{\text{s}}} }}} \right|_{{P_{{\text{e}}} = {\text{low}}}} - \left. {\frac{{V_{{\text{p}}} }}{{V_{{\text{s}}} }}} \right|_{{P_{{\text{e}}} = {\text{high}}}}.$$

(8)

For dry carbonate samples in our measurements, Δ(V_p/V_s) < 0 means the V_p/V_s normally increases, while Δ(V_p/V_s) > 0 means the V_p/V_s abnormally decreases. As reported at the beginning of Sect. 4.2, the low and high effective pressures considered for the PDV analysis of these carbonate samples are 6 and 35 MPa, respectively. The results are shown in Fig. 12. Unlike the V_p-V_s linear relation discussed in the above section, the V_p/V_s is more responsive to pressure change. According to our measurements, when the main axis of the sample is normal to the bedding plane, the V_p/V_s of 18 carbonate samples normally increases with the pressure increasing from 6 to 35 MPa, while the V_p/V_s of the other three samples slightly decreases. However, when the main axis of the sample is oblique or parallel to the bedding plane, the V_p/V_s of only three carbonate samples normally increases with the pressure, while the V_p/V_s of the other seven sub-cores abnormally decreases. Thus, sample orientation likely played a role in the abnormal behavior of V_p/V_s recorded on one-third of the studied Weyburn carbonate samples. This is quite a high number of samples, and the observed abnormal trends cannot be explained solely by measurement errors.

Fig. 12

Changes in V_p/V_s on dry carbonate rocks following the increase in the effective stress. The histogram subplot represents the number of samples with different levels of V_p/V_s changes. The red points are based on our measurements. Measurements carried out under similar conditions by Fournier et al. (2014) and Hairabian et al. (2014) are also included. Here, pore volume change caused by increasing effective stress is neglected

To put our results into perspective, the measurements on dry carbonate samples published by Fournier et al. (2014) and Hairabian et al. (2014) are also included in Fig. 12. The additional data include 251 carbonate samples that are variably porous and which display a large variety of pore types and textures. It should be emphasized that both Fournier et al. (2014) and Hairabian et al. (2014) described that their samples were “vertically oriented, cylindrical samples.” The uncertainties in their velocity measurements were around 0.5% to 1%.

Overall, the Δ(V_p/V_s) scatters widely and this is independent of the pore volume (Fig. 12). Among the dry vertical-oriented carbonate samples, 229 samples (about 83.9% percent of samples analyzed here) show a Δ(V_p/V_s) smaller than 5% (i.e., |Δ(V_p/V_s)|< 0.05), whereas 101 samples (about 37.0% percent of samples analyzed here) show an abnormal decrease in V_p/V_s with the increase in the effective pressure. These different trends clearly imply that using V_p/V_s as a tool to predict abnormal (pore) pressure in carbonate gas reservoirs requires careful calibration.

4.2.3 Normalized velocity

Since the amount of velocity change is dependent on the pressure range, the quantitative comparison of the PDV of different carbonate samples should be constrained within the same pressure window. For our measurements, the chosen lower and upper effective pressure bounds for the PDV analysis are 6 and 35 MPa, respectively. Indeed, within this pressure range, the P- and S-wave velocities of all the studied rock specimens were reliably collected. To quantitatively analyze the PDV, the normalized velocities (V_pn and V_sn) are defined in the following:

$$V_{{{\text{pn}}}} { = }\frac{{V_{{\text{p}}} \left( {P_{{\text{e}}} = 35\,{\text{MPa}}} \right)}}{{V_{{\text{p}}} \left( {P_{{\text{e}}} = 6\,{\text{MPa}}} \right)}}$$

(9)

and

$$V_{{{\text{sn}}}} = \frac{{V_{{\text{s}}} \left( {P_{{\text{e}}} = 35\,{\text{MPa}}} \right)}}{{V_{{\text{s}}} \left( {P_{{\text{e}}} = 6\,{\text{MPa}}} \right)}}.$$

(10)

As shown in Fig. 1, rock velocity usually increases with pressure; thus, the values of the normalized velocities defined in Eqs. (9) and (10) are higher than 1.0. The normalized velocities quantitatively represent the PDV of these dry carbonate samples because they are calculated at the same pressure window.

As shown in Fig. 13, the normalized velocities do not show clear dependence to the pore volume. V_pn is slightly higher than V_sn. Compared to the linear V_p–V_s relation in Sect. 4.2.1 (Fig. 11), the linear fit between V_pn and V_sn shows stronger scattering, but it is still a significant linear relation. Indeed, the critical linear correlation coefficient (R_{0.01, 20}) for a fitting with 20 degrees of freedom (i.e., 22 data points as in Fig. 13) and with a 0.99 significance level is 0.537. The latter is much smaller than the correlation coefficient R obtained here, which is 0.825. The V_pn and V_sn have similar distributions: (i) V_pn ranges between 1.009 and 1.142 with the mean value of 1.043, while V_sn ranges between 1.004 and 1.095 with the mean value of 1.032; (ii) the V_pn and V_sn of more than 72% samples are smaller than 1.05, which means the increase in velocity is less than 5%. According to the measurements of Freund (1992) on 88 dry clastic rocks, the mean values of V_pn and V_sn calculated by the velocities at the effective pressure of 8 MPa and 40 MPa are 1.196 and 1.127, respectively. Thus, the normalized velocity for carbonate rocks is obviously smaller than that of sandstones at similar pressure range. This means that the PDV of carbonate rocks is much weaker than that of sandstones, and velocities of carbonate rocks are less sensitive to abnormal pore pressure. Pore pressure prediction of carbonate reservoirs based on velocity–effective stress relation may not be as reliable as that in sandstone formations.

Fig. 13

The normalized P- and S-wave velocities with different porosities. a Plotted here are our laboratory measurements on carbonates and the measurements of Freund (1992) on sandstones. b Zooming-in on our measurements (see red rectangle in (a)) to have a detailed view of their statistical distributions. Numbers along the gray lines indicate their slopes. The black line is the linear fitting of the data points. The histograms along the x- and y-axis are the distributions of V_pn and V_sn, respectively

4.3 Empirical velocity–pressure models

To test the suitability of the two types of empirical models reported in Sect. 1 for carbonate rocks, we use the elastic wave velocities measured on the 22 Weyburn carbonates during the loading process. For the exponential model (Eqs. (1) and (2)), the measured velocities at different pressures are fitted by Eq. (2), i.e., without the linear item. This is because our measurements were carried out under relatively low effective pressures. For the two empirical models (Eqs. 2 and 4), we use the same fitting method and the same database to do the fittings. The average value of the squared correlation coefficients (R²) for the fittings of the carbonate samples is higher than 0.97 (Table 3). Predicted velocities on all the carbonate samples based on the two empirical models are compared to our measurements in Fig. 14.

Table 3 The statistics of the squared correlation coefficients (R²) for fittings by different models

Fig. 14

Comparison between the measured velocities and the predictions of the empirical models for all carbonate samples at different effective pressures. a Measurements vs. the exponential model data. b Measurements vs. the power function model data

As shown in Table 3 and Fig. 14, the fitting results indicate that the two models are both effective in describing the PDV for carbonate rocks at relative low pressure (P_e < 60 MPa). Only few samples did not follow the variation trends obtained through these two empirical models. The low quality of the waveforms recorded on those few samples likely explains such discrepancy. Also, velocities measured at different effective pressures on this subset of samples sometimes show very irregular fluctuations. As shown in Fig. 14, the exponential model (Eq. (2)) slightly performs better than the power function model (Eq. (4)). Also, the R² of fittings from the exponential model has slightly higher means and lower standard deviations (Table 3). However, for many applications, the differences between the fitting results of these two models are too small and can be ignored safely.

The coefficient of the exponent (D) in Eq. (2) represents the level of PDV. The fitted values of D derived from V_p and V_s are statistically the same (Fig. 15), regardless of the dominant mineral composition or the orientation of the sample main axis. This is in line with sandstones data listed in Table 2 of Khaksar et al. (1999, see also Fig. 15). Thus, both our observations and those of Khaksar et al. (1999) support the suggestion of Shapiro (2003) that all coefficients D are identical in the first approximation. The D values for our set of carbonate samples are obviously more scattered than those of the sandstones reported in Khaksar et al. (1999), which confirms the complex behavior of the PDV of carbonate rocks.

Fig. 15

The fitting parameter D of the exponential model for P- and S-wave velocities (D-V_p and D-V_s, respectively) for a the vertical samples, here color-coded based on their main minerals (the data of the quartz-rich samples (i.e., sandstones) are from Khaksar et al.. (1999)); and b samples with main axis at different directions to the bedding plane

5 Discussion

5.1 Mineral composition and PDV

While seismic velocity of natural rocks, especially clastic rocks, strongly depends on their pore volumes, it is quite different for carbonate rocks. Neither the V_p-V_s ratio change with pressure (Fig. 12), nor the normalized velocity (Fig. 13) shows apparent correlations with the pore volume. This is because pore volume reduction induced by pressure increase is quite small for most carbonate rocks. The stiffness of the rock skeleton, more specifically some factors controlling it such as the mineral composition and the pore structure, would play more influential roles on PDV than the pore volume.

In order to assess the influence of the mineralogy on the PDV of rocks, the V_pn at P_e = 35 MPa is plotted against the pore volume in Fig. 16; therein, rocks with different dominant minerals are distinguished. For the 24 studied samples, 11 samples are dolomite-rich, 11 samples are calcite-rich and 2 samples are quartz-rich (sandstone). As discussed above, the normalized velocity change seems independent of the pore volume (Fig. 13). Thus, it is not surprising that the V_pn-φ cross-plot is quite scattered (Fig. 16). The linear fitting coefficient is 0.57 and 0.43 for dolomite-rich and calcite-rich samples, respectively, which are smaller than the critical linear correlation coefficient (R_{0.01, 9} = 0.735) for a fitting with nine degrees of freedom (i.e., 11 data points) and 0.99 significance level. There is no noticeable difference in the V_pn-φ between the dolomite- and calcite-rich samples.

Fig. 16

The normalized P-wave velocity for the Weyburn samples considering the difference in mineral compositions. The dashed red line is the linear fit for the dolomite-rich samples. The dashed green line is the linear fit for the calcite-rich samples. The fitting coefficient (R) for all these fittings is shown with corresponding colors. The blue stars represent dry clastic sandstones tested by Freund (1992); the low and high effective pressures considered in their calculation of V_pn (see Eq. 8) are 8 and 40 MPa, respectively

The V_pn-φ cross-plot of the sandstones from our dataset cannot be interpreted with enough certainty due to less data points. Therefore, for more reliability, the measurements of dry clastic sandstones by Freund (1992) are added here (Fig. 16). On the one hand, as a whole the V_pn of the sandstones is much higher than that of carbonate rocks, which shows that the sandstones have stronger PDV than the carbonate rocks. On the other hand, the V_pn-φ cross-plot of sandstones displays a wide scattering. The carbonate rocks data points show some scattering as well but to a smaller degree. In light of experimental studies conducted on sandstones, the presence of clay minerals can be a factor causing the scattering of their PDV (Eberhart-Phillips et al. 1989; Freund 1992; Sams and Andrea 2001). As for carbonate rocks, a similar observation has not been made so far. In our study, overall there is no evidence to suggest that the observed V_pn-φ scattering for carbonate samples is due to non-carbonate minerals (quartz, feldspar, anhydrite) typically present in their solid frames.

5.2 Pore structure and PDV

The V_pn data were grouped by formation unit to assess how pore geometries may influence the PDV under similar sedimentary conditions (Fig. 17). The black dashed lines in the middle subplots of Fig. 17 indicate the average curves for different formations. The left and right subplots are, respectively, the pore throat distributions of samples showing maximum and minimum V_pn in the same formation. It is easy to find that the dominant pore throat diameters for maximum and minimum V_pn have significant differences. For the same formation, the variation of PDV is probably caused by the matrix texture and the pore structure.

Fig. 17

Comparisons of the pressure dependence with respect to pore structures in different formation units. For the middle subplots, different color curves represent different samples. Samples showing maximum and minimum V_pn in each formation unit are indicated by the red and blue colors, respectively, with the sample index marked accordingly in the same color. The black dashed lines indicate the average values of V_pn for samples in the same formation unit. The left and right subplots are, respectively, the pore throat distributions of samples showing maximum (red color) and minimum (blue color) V_pn in the same formation unit

For each unit, the sample with minimum V_pn has a small dominant pore throat diameter. The latter is usually not larger than 0.1 µm (the right subplots of Fig. 17a, b and d), except for the loose Marly unit. A small dominant pore throat diameter is often typical for tightly packed matrixes.

Moreover, the average value of the normalized velocity change in the loose Marly unit is 1.09 at P_e = 35 MPa (Fig. 17c), while it is smaller than 1.06 for other formations (Fig. 17a, b and d). The PDV of the loose Marly unit is significantly larger than those of the other formations due to its loosely packed matrix. Interestingly too, for the loose Marly unit, the samples with the maximum and minimum V_pn, samples 20 and 25, respectively, almost have the same main pore throat diameter (1 to 2 µm; Fig. 17c). Yet, the difference in their pore volume, 27% and 20% for samples 20 and 25, respectively, is possibly the main reason behind their different PDV. Thus, within a carbonate formation made of loosely packed material, pore volume can well be considered as a PDV main controlling factor.

6 Conclusions

Elastic compressional and shear velocities were measured on 24 dry carbonate specimens following loading and unloading processes under hydrostatic condition. The rock samples were selected from seven wells of the Weyburn oilfield, and they display a wide range of pore volume, pore type and depositional texture. Making use of the data collected on the carbonate rocks forming the reservoir, the pressure dependence of velocities (PDV) under dry conditions were investigated. The results indicate that:

• (1) P- and S-wave velocities in the dry Weyburn carbonate specimens follow a linear relation that is independent of the pressure change. The ratio V_p/V_s is more responsive to pressure change irrespective of the pore volume of the carbonate rocks. One-third of the dry Weyburn carbonate samples show abnormal decrease in V_p/V_s at higher effective pressure. Thus, using V_p/V_s to predict abnormal pore pressure may not be reliable for carbonate reservoirs.

• (2) The exponential model V = $A-C\text{e}^{DP_{\text{e}}}$ gives good fits for most of the Weyburn carbonate samples. Statistically, there is no difference between the V_p- and V_s-derived exponential coefficient D.

• (3) The PDV of the carbonate rocks varies widely even for samples from the same formation with similar sedimentary history. The highest PDV are obtained on samples displaying loosely packed crystals and/or a pore network that is mostly controlled by large pores. No evidence points to the PDV being influenced by the presence of non-carbonate minerals.

Appendices

Appendix 1: Supporting data

See Tables 4 and 5.

Table 4 The location, density, porosity and permeability of all samples*

Table 5 The ultrasonic P- and S-wave velocities measured at different confining pressures

Appendix 2: The differential effective medium (DEM) model

The differential effective medium (DEM) model is widely used to predict the effective bulk (K^*) and shear (μ*) moduli of two-phase composites. It is a coupled system of two ordinary differential equations (Berryman 1992):

$$\left\{ \begin{gathered} \left( {1 - y} \right)\frac{{{\text{d}}\left[ {K^{ * } \left( y \right)} \right]}}{{{\text{d}}y}} = \left( {K_{{\text{f}}} - K^{ * } } \right)P^{{ * {\text{f}}}} \left( y \right) \hfill \\ \left( {1 - y} \right)\frac{{{\text{d}}\left[ {\mu^{ * } \left( y \right)} \right]}}{{{\text{d}}y}} = \left( {\mu_{{\text{f}}} - \mu^{ * } } \right)Q^{{ * {\text{f}}}} \left( y \right) \hfill \\ \end{gathered} \right.$$

(11)

with the initial conditions

$$\left\{ \begin{gathered} K^{ * } \left( 0 \right) = K_{{\text{m}}} \hfill \\ \mu^{ * } \left( 0 \right) = \mu_{{\text{m}}} \hfill \\ \end{gathered} \right.$$

(12)

where K and μ are the bulk and shear modulus, respectively, m and f stand for the solid matrix and pore fluid, respectively, and y is the pore volume fraction. The parameters P^*f and Q^*f for elliptical pores of arbitrary aspect ratio (AR) are given by

$$P^{{*{\text{f}}}} = \frac{1}{3}T_{1} ,$$

(13)

and

$$Q^{{\text{*f}}} = \frac{1}{5}\left( {T_{2} - \frac{1}{3}T_{1} } \right).$$

(14)

T₁ and T₂ can be calculated as follows (Berryman 1980):

$$T_{1} = {{3F_{1} } \mathord{\left/ {\vphantom {{3F_{1} } {F_{2} }}} \right. \kern-\nulldelimiterspace} {F_{2} }}$$

(15)

and

$$T_{2} - \frac{1}{3}T_{1} = \frac{2}{{F_{3} }} + \frac{1}{{F_{4} }} + \frac{{F_{4} F_{5} + F_{6} F_{7} - F_{8} F_{9} }}{{F_{2} F_{4} }},$$

(16)

wherein

$$\begin{gathered} F_{1} = 1 + A\left[ {\frac{3}{2}\left( {f + \beta } \right) - R\left( {\frac{3}{2}f + \frac{5}{2}\beta - \frac{4}{3}} \right)} \right], \hfill \\ F_{2} = 1 + A\left[ {1 + \frac{3}{2}\left( {f + \beta } \right) - \frac{R}{2}\left( {3f + 5\beta } \right)} \right] + B\left( {3 - 4R} \right) + \frac{A}{2}\left( {A + 3B} \right)\left( {3 - 4R} \right)\left[ {f + \beta - R\left( {f - \beta + 2\beta^{2} } \right)} \right], \hfill \\ F_{3} = 1 + A\left[ {1 - \left( {f - \frac{3}{2}\beta } \right) + R\left( {f + \beta } \right)} \right], \hfill \\ F_{4} = 1 + \frac{A}{4}\left[ {f + 3\beta - R\left( {f - \beta } \right)} \right], \hfill \\ F_{5} = A\left[ { - f + R\left( {f + \beta - \frac{4}{3}} \right)} \right] + B\beta \left( {3 - 4R} \right), \hfill \\ F_{6} = 1 + A\left[ {1 + f - R\left( {f + \beta } \right)} \right] + B\left( {1 - \beta } \right)\left( {3 - 4R} \right), \hfill \\ F_{7} = 2 + \frac{A}{4}\left[ {3f + 9\beta - R\left( {3f + 5\beta } \right)} \right] + B\beta \left( {3 - 4R} \right), \hfill \\ F_{8} = A\left[ {1 - 2R + \frac{f}{2}\left( {R - 1} \right) + \frac{\beta }{2}\left( {5R - 3} \right)} \right] + B\left( {1 - \beta } \right)\left( {3 - 4R} \right), \hfill \\ F_{9} = A\left[ {\left( {R - 1} \right)f - R\beta } \right] + B\beta \left( {3 - 4R} \right), \hfill \\ \end{gathered}$$

with A, B and R given by:

$$\begin{gathered} A = \frac{{\mu_{{\text{f}}} }}{{\mu^{ * } }} - 1, \hfill \\ B = \frac{1}{3}\left( {\frac{{K_{{\text{f}}} }}{{K^{ * } }} - \frac{{\mu_{{\text{f}}} }}{{\mu^{ * } }}} \right), \hfill \\ R = \frac{{1 - 2\nu^{ * } }}{{2\left( {1 - \nu^{ * } } \right)}}, \hfill \\ \end{gathered}$$

where v^* is the effective Poisson’s ratio of the composite mineral. The functions β and f are given by:

$$\begin{gathered} \beta = \frac{AR}{{\left( {1 - AR^{2} } \right)^{{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-\nulldelimiterspace} 2}}} }}\left[ {\cos^{ - 1} AR - AR\left( {1 - AR^{2} } \right)^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}}} } \right], \hfill \\ f = \frac{{AR^{2} }}{{1 - AR^{2} }}\left( {3\beta - 2} \right). \hfill \\ \end{gathered}$$