Experimental Investigation on Static and Dynamic Bulk Moduli of Dry and Fluid-Saturated Porous Sandstones

Abstract

Knowledge of pressure-dependent static and dynamic moduli of porous reservoir rocks is of key importance for evaluating geological setting of a reservoir in geo-energy applications. We examined experimentally the evolution of static and dynamic bulk moduli for porous Bentheim sandstone with increasing confining pressure up to about 190 MPa under dry and water-saturated conditions. The static bulk moduli (K_s) were estimated from stress–volumetric strain curves while dynamic bulk moduli (K_d) were derived from the changes in ultrasonic P- and S- wave velocities (~ 1 MHz) along different traces, which were monitored simultaneously during the entire deformation. In conjunction with published data of other porous sandstones (Berea, Navajo and Weber sandstones), our results reveal that the ratio between dynamic and static bulk moduli (K_d/K_s) reduces rapidly from about 1.5 − 2.0 at ambient pressure to about 1.1 at high pressure under dry conditions and from about 2.0 − 4.0 to about 1.5 under water-saturated conditions, respectively. We interpret such a pressure-dependent reduction by closure of narrow (compliant) cracks, highlighting that K_d/K_s is positively correlated with the amount of narrow cracks. Above the crack closure pressure, where equant (stiff) pores dominate the void space, K_d/K_s is almost constant. The enhanced difference between dynamic and static bulk moduli under water saturation compared to dry conditions is possibly caused by high pore pressure that is locally maintained if measured using high-frequency ultrasonic wave velocities. In our experiments, the pressure dependence of dynamic bulk modulus of water-saturated Bentheim sandstone at effective pressures above 5 MPa can be roughly predicted by both the effective medium theory (Mori–Tanaka scheme) and the squirt-flow model. Static bulk moduli are found to be more sensitive to narrow cracks than dynamic bulk moduli for porous sandstones under dry and water-saturated conditions.

1 Introduction

The elastic moduli of rocks are important for many geotechnical applications, such as in design and construction of rock engineering projects, and in interpretation of seismic data in geophysical explorations. There are, in general, two main methods used to determine the elastic moduli of rocks. Specifically, static moduli are derived from stress–strain relations in quasi-static rock compression tests performed on intact samples with very low loading rate (i.e., static method). Alternatively, dynamic moduli of rocks are determined using measurements of ultrasonic wave velocities transmitting through the samples (Cheng and Johnston 1981). Although dynamic elastic moduli are typically estimated using high-frequency ultrasonic wave velocities measurements, it is estimates of static elastic moduli that are widely employed for evaluating rock deformation, especially in rock excavation projects. However, static and dynamic moduli of rock samples are commonly found to be different, thus understanding the underlying physical causes for the observed difference is crucial for utilizing both static and dynamic measurements.

For nonporous, homogeneous materials like metal and glass, dynamic and static moduli are experimentally found to be equal (Simmons and Brace 1965; Ledbetter 1993). In contrast, the dynamic modulus is almost always larger than the static one for many rocks (Simmons and Brace 1965; King 1969; Cheng and Johnston 1981; Coyner 1984; Fjær 2009, 2019). Recently, Fjær (2019) summarized some potential reasons for the discrepancy between static and dynamic moduli of rocks, including the different strain rates, drainage conditions, as well as the heterogeneity of rock properties.

For dry porous rocks, measurements reveal that the difference between dynamic and static moduli is largely suppressed by increasing pressure with values converging at high pressure (Simmons and Brace 1965; King 1969, 1983; Cheng and Johnston 1981; David et al. 2013). Although the difference in frequency and induced strain amplitude can well predict that the dynamic modulus is larger than the static modulus (Martin and Haupt 1994; Fjær 2009), it appears to fail in accounting for the observed pressure-dependent discrepancy between static and dynamic bulk moduli of dry porous rocks. Three-dimensional microstructural observations on porous sandstones reveal that the whole void space can be divided into three distinct categories: equant nodal pores located at four-grain vertices, tubular pores at three-grain edges, and narrow inter-granular cracks at two-grain interfaces (Bernabe 1991; Fredrich et al. 1995). Void structures are highly variable, but to first order, the total porosity of porous sandstones may be divided into stiff porosity (i.e., the former two types of void space aforementioned) and compliant porosity (i.e., crack porosity) (Shapiro 2003; Wang et al. 2020a). Under pressure, stiff pores exhibit little deformation. Considering compliant cracks that can be easily closed with increasing pressure (Walsh 1965a), the correlation between dynamic and static moduli of porous sandstones is very likely related to the void structure.

For fluid-saturated porous rocks, the frequency-dependent interaction between the rock skeleton and the pore fluid may affect the dynamic bulk modulus measured by ultrasonic wave. When the frequency is high, wave-induced fluid pressure variation may not equilibrate at the representative elementary volume (REV) scale. The fluid can then be considered to be partly mobile or even immobile and thus, higher pressures are maintained in pore structures leading to a higher compressional stiffness of fluid-saturated rocks (O’Connell and Budiansky 1977; Dvorkin et al. 1994, 1995; Gurevich et al. 2010). Conversely, static bulk moduli of fluid-saturated rock samples measured during the quasi-static loading are similar to moduli measured at oven-dry conditions because the fluid pressure has sufficient time to equilibrate in the REVs (i.e., drained conditions) (Gurevich et al. 2010). Therefore, the difference between dynamic and static bulk moduli is expected to be enhanced by fluid saturation. The measurements of ultrasonic frequency P- and S-wave velocities of fluid-saturated porous sandstones with increasing pressure indicate that dynamic bulk moduli are considerably increased at low pressures, and show only minor changes at high pressures (Coyner 1984; Fortin et al. 2007; David et al. 2013). This suggests that the different mechanical responses of void structures to pressure may lead to the distinct pore pressure distributions, and ultimately affect the macroscopic dynamic and static bulk moduli of fluid-saturated porous sandstones. However, the evolution of void structures with pressure was mostly inverted from bulk moduli based on several empirical assumptions on void geometry and on the relation between bulk compressibility and external pressure (Zimmerman 1990; Shapiro 2003; David and Zimmerman 2012; Zhang et al. 2019), rather than experimental determination. Hence, the quantitative relation between static and dynamic bulk moduli of porous sandstones and real evolution of void structures still remains ambiguous.

The purpose of this study is to improve our understanding of static and dynamic bulk moduli of porous sandstones under dry and fluid-saturated conditions from a perspective of laboratory investigation. A series of hydrostatic compression tests under dry and water-saturated conditions were performed on porous Bentheim sandstone samples, and the static and dynamic bulk moduli were continuously measured as a function of pressure. Based on the measured evolution of pore space with increasing pressure, we quantitatively evaluated the role of different pore geometries on the discrepancy between static and dynamic bulk moduli of dry and water-saturated Bentheim sandstone samples.

2 Materials

The sample material used in this study is porous Bentheim sandstone from a shallow marine formation deposited during the Lower Cretaceous. The samples were retrieved from outcrops near Bentheim, Germany. The material is a homogeneous, isotropic porous sandstone, composed of 96.5% quartz, 2% feldspar and 1.5% kaolinite. As a result of depositional environment, Bentheim sandstone is composed of well-sorted mostly rounded to sub-rounded quartz grains (Fig. 1). The diameter of quartz grains varies from 0.05 to 0.55 mm, with median grain size ranging between 0.20 and 0.33 mm. The average initial connected porosity of samples used in this study is about 0.233 and the permeability measured at ambient pressure is about 1 Darcy (Wang et al. 2020a). In addition to interconnected equant pores surrounded by adjacent grains, many narrow microcracks connected to equant pores are observed mostly at grain–grain contacts (inter-granular cracks) or partly as intra-granular cracks, as illustrated in Fig. 1.

Fig. 1

Thin-section optical micrographs of Bentheim sandstone impregnated with blue epoxy shown in different magnifications (a, b). Scanning electron microscopy (SEM) images of grains in Bentheim sandstone (c, d). Thin inter-granular cracks (red arrows) at gain-grain contact and equant pores (orange arrows) surrounded by several neighbouring grains are clearly observed. The cracks and pores are well interconnected forming a continuous network

3 Methods

3.1 Experimental Procedures

3.1.1 Sample Preparations and Strain Measurements

The size of cylindrical samples cored from a single block was 50 mm in diameter and 100 mm in length. The samples were kept in an oven at 50 °C temperature for at least 24 h before testing. Two pairs of orthogonally oriented electric resistance strain gages (Tokyo Sokki TML-FCB) allowed monitoring axial stain ε₁ and radial strain ε₃. Each strain gage (120 Ω) was mounted in a one-fourth Wheatstone bridge. The volumetric strain ε_b (ε_b = ΔV_b/V_bⁱ where V_bⁱ and ΔV_b are the initial bulk volume of the rock sample and the changed bulk volume of the rock sample with pressure, respectively) is calculated using the expression ε_b = ε₁ + 2ε₃, in which the axial strain ε₁ is the average value of two vertically oriented strain gages, whereas 2ε₃ represents the sum of tangential strain values measured by two strain gages oriented in circumferential direction (see Fig. 11a in Appendix A). Two pairs of orthogonally oriented strain gages were placed in axial and circumferential directions onto a thin layer of fast curing adhesive that was applied to fill up the surface pores. Neoprene jackets were used to isolate the samples from the oil confining medium. It is worth noting that compressive stress and compressive volumetric strain are considered positive in this study.

3.1.2 Ultrasonic Wave Measurements

To record ultrasonic wave velocities along different traces, twelve P-wave and four polarized S-wave piezoelectric transducers (PZTs) contained in brass housings were glued directly to the surface of the rock and sealed in the neoprene jacket using a two-component epoxy (see Fig. 11b in Appendix A). Piezoelectric P- and S-wave sensors have a resonant frequency of ~ 1 MHz. After assembling, jacketed samples were placed in the oven again for at least 12 h at 50 °C for drying.

Six out of twelve P-wave PZTs and two of four polarized S-wave PZTs, were used as ultrasonic transmitters. These transducers were used to emit every 30 s a rectangular electrical pulse with 100 V amplitude and 3 μs duration. Meanwhile, the time span between each ultrasonic transmission was 5 ms. Ultrasonic pulses were recorded by the remaining transducers forming 36 different transmitter–receiver traces for P-wave velocity measurements as well as 4 different S-wave traces. Transducer signals were amplified by 40 dB using Physical Acoustic Corporation (PAC) preamplifiers equipped with 100 kHz high-pass filters. Ultrasonic signals were stored in a 16-channel transient recording system with an amplitude resolution of 16 bit at a sampling rate of 10 MHz, corresponding to ± 0.1 μs uncertainty for P-wave and S-wave arrival times. The arrival time of ultrasonic waveform was picked using a series of picking algorithms including the Akaike information criterion (Wang et al. 2020b). Net travel time through the samples was determined by means of correcting the arrival time for delay in the brass housing or in the loading plates. Also, the ultrasonic velocity measurements were corrected for the deformation of samples affecting the distances between PZTs.

3.1.3 Mechanical Compression Tests

We carried out a series of hydrostatic compression tests summarized in Table 1. Experiments were performed in a servo-hydraulic 4600 kN loading frame (MTS) equipped with a pore pressure system (Quizix 6000) (see Fig. 11c in Appendix A). All tests were performed at room temperature on dry and water-saturated samples.

Table 1 Summary of conducted experiments

In two hydrostatic compression tests performed on oven-dry specimens, samples were installed in the pressure vessel and subjected to a vacuum of about 10^–2 bar at 2 MPa confining pressure for 12 h prior to starting the experiments and zero readings were taken from the strain gages. Subsequently, the hydrostatic pressure was increased up to 190 MPa at a rate of 1 MPa/min. The maximum confining pressure was maintained for 1 min, followed by hydrostatic unloading at a similar rate of 1 MPa/min. Vacuum was maintained during the entire experiment.

For two repeated experiments conducted at drained conditions, dry samples were saturated with distilled water at pore pressure of P_p = 2 MPa and confining pressure of P_c = 3 MPa for several hours until samples were completely water-saturated. Afterwards, pore pressure was kept constant at 2 MPa during the entire deformation of sample at the imposed pressurization rate of 1 MPa/min. Based on the sample initial permeability (κ ≈10^–12 m²), the sample length (L = 100 mm), water viscosity (η ≈ 10^–3 Pa s) and bulk compressibility of water (C_f≈ 0.5 GPa⁻¹), the characteristic diffusion time t_c for fluid to equilibrate after perturbations across the sample may be computed by t_c = L²ηC_f /κ (Mavko et al. 2009; Wang et al. 2020b). The estimated diffusion time t_c < 5 × 10^–3 s, far shorter than the experimental duration, indicates that the fluid pressure within the sample equilibrates rapidly, and thus the complete drained condition is met. With changing confining pressure and constant fluid pressure, fluid volume changes are continuously monitored by a Quizix pump volumometer with an accuracy of about ± 0.001 cm³. The volume changes are assumed to reflect changes in volume of interconnected void space. The concept of Terzaghi effective pressure (i.e., P_e = P_c − P_p) was adopted to report the results of drained tests.

3.2 Calculations of Static and Dynamic Bulk Moduli

3.2.1 Static Bulk Moduli

In this study, we monitored changes in static and dynamic bulk moduli of Bentheim sandstone with changing hydrostatic pressures. In general, during initial hydrostatic loading of the samples, a non-linear increase in volumetric strain with pressure is observed, which is attributed to the progressive closure of pre-existing compliant microcracks or the non-linear contact law between grains (Walsh 1965a; Mavko et al. 2009; Wang et al. 2020a). The static bulk modulus K_s is defined as the local slope of the confining pressure (P_c) vs. volumetric strain (ε_b) curve, as expressed by

$$K_{{\text{s}}} = \frac{{dP_{{\text{c}}} }}{{d\varepsilon_{b} }}.$$

(1)

To calculate K_s, we first smoothened raw data by applying moving average window of 21 points with a sampling rate of 2 Hz (i.e., the reading record at every 1/120 MPa interval for confining pressure). Afterwards, a high-order polynomial was fitted globally to the entire pressure–volumetric strain curve, and calculated the derivative of the confining pressure‒volumetric strain curve to obtain K_s.

3.2.2 Dynamic Bulk Moduli

The dynamic bulk modulus K_d is obtained from the P- and S-wave velocities using the common expression for isotropic rocks:

$$K_{d} = \rho \left( {W_{p}^{2} - \frac{4}{3}W_{s}^{2} } \right),$$

(2)

where ρ is the bulk density of rock sample. P-wave velocity W_p was taken as the average of five horizontal wave velocities denoted by W_p,radial and one axial wave velocity denoted by W_p,axial. S-wave velocity W_s was calculated using the mean value of all measured horizontally and vertically polarized S-wave velocities (denoted by W_sh and W_sv, respectively). The measurements of ultrasonic wave velocity were corrected for the changes in the distance among the PZTs associated with deformation. The bulk density of water-saturated samples is estimated by ρ = ϕρ_w + (1‒ ϕ)ρ_g where ρ_w is water density, ϕ is porosity and ρ_g is the density of grains (≈ 2.61 g/cm³). Furthermore, we also corrected the changes in densities of water and grains with increasing pressure, and the resulting changes were found to be negligible.

4 Results

4.1 Sample Deformation and Ultrasonic Wave Velocities

For dry and drained compression tests, axial strain, radial stain and volumetric strain increase nonlinearly at low pressure below ~ 60 MPa, and then grow linearly at higher pressure (Fig. 2). Similarly, P-wave and S-wave velocities for dry and water-saturated samples increase nonlinearly at pressure up to ~ 60 MPa, followed by a minor rise with increasing pressure (Fig. 3). Good agreement between axial strain and radial strain suggests that the Bentheim sandstone samples are deformed isotropically (Fig. 2), supported by the observed minor difference between P-wave velocities transmitting in vertical and horizontal directions (Fig. 3). Thus, for Bentheim sandstones, we may rule out the potential influence of rock property heterogeneity on the difference between static and dynamic bulk moduli. In addition, the irreversible volumetric strain after unloading accounts for about 5% of overall strain. For water-saturated samples, the P-wave velocities are about 10% larger than dry samples at P_e = 2 MPa and about 2% at P_e = 180 MPa, respectively. In contrast, corresponding S-wave velocities for water-saturated samples are slightly lower than for dry samples.

Fig. 2

The measured axial strain, radial strain (a) and volumetric strain of bulk sample (b) with confining pressure in Dry Test 1. The measured axial strain, radial strain (c) and volumetric strain of bulk sample (d) with effective pressure in Drained Test 1. Solid curves indicate loading phase while dashed ones represent unloading path. For simplicity, here we do not show the results of Dry Test 2 and Drained Test 2, which yielded quite similar results

Fig. 3

Measured P- and S-wave velocities as a function of pressure in a Dry Test 1 and b Drained Test 1, respectively. Solid curves indicate loading phase while dashed ones represent unloading path. The error bar for P-wave velocities indicates standard deviation of W_p,axial and W_p,radialwhereas the error bar for S-wave velocities denotes the standard deviation of W_sh and W_sv

4.2 Total Porosity and Crack Porosity

The induced grain cracking for Bentheim sandstone under hydrostatic loading is not expected to occur because the corresponding critical pressure denoting the onset of grain crushing and pore collapse was experimentally found to be ~ 440 MPa (Tembe et al. 2008), far more than our applied maximum pressure. Hence, the whole void space of Bentheim sandstone is expected to be squeezed in response to hydrostatic pressure, leading to a gradual reduction of total porosity. The measured drained water volume (ΔV_p) is increased nonlinearly to about 1.5 mL at P_e = 60 MPa (Fig. 4a), and the following linear evolution of ΔV_p with pressure reveals a linear response of remaining stiff pores to hydrostatic pressure.

Fig. 4

Evolution of a drained water volume, b total porosity ϕ and c crack porosity ϕ_c with increasing effective pressure during loading in drained compression tests. The crack porosity ϕ_c is determined by subtracting linear evolution of stiff porosity ϕ_s from total porosity ϕ (Wang et al. 2020a). Error bars indicate the maximum and minimum of two drained compression tests

The evolution of total porosity ϕ is given by ϕ = V_p/V_b = ϕ_i(1‒ε_p)/(1‒ε_b) where V_b and V_p are the pressure-dependent bulk volume of the rock sample and the pressure-dependent volume of void space, respectively. ϕ_i is the initial porosity and ε_p = ΔV_p/V_pⁱ (V_pⁱ is the initial volume of void space within the sample) is the volumetric strain of void space. As stated above, we assume that total porosity may be represented by a combination of compliant cracks and stiff pores. Based on the linear extrapolation of stiff porosity (ϕ_s) at high pressures (all cracks are closed), the crack porosity (ϕ_c) at low pressures can be estimated by subtracting the stiff porosity from total porosity (Wang et al. 2020a) (Fig. 4b). Apparently, the initial crack porosity constitutes only a very small fraction of the total porosity and decreases rapidly from about 0.47% at P_e = 1 MPa to zero at P_e ≈ 60 MPa (Fig. 4c). Beyond the closure of all cracks, stiff porosity decreases linearly with increasing pressure. The crack closure pressure P_cl is defined as a characteristic hydrostatic pressure above which almost all cracks are closed (Walsh 1965a). The crack closure pressure for Bentheim sandstone is ~ 60 MPa (Fig. 4b, c), confirmed by the observed evolution of P- and S-wave velocities with pressure. Since the deformation of void space at dry conditions is expected to be similar to that in drained tests at the same pressure (Mavko et al. 2009), hereafter the porosity data measured in drained tests were used to characterize the deformation of void space in dry compression tests.

4.3 Static and Dynamic Bulk Moduli of Dry and Water-Saturated Porous Sandstone Samples

Since the Bentheim sandstone samples were mainly deformed elastically over an entire range of applied pressure, here we only reported the static and dynamic bulk moduli measured during loading. For hydrostatic compression tests performed on dry specimens, static and dynamic bulk moduli increased significantly from K_s ≈ 5 GPa and K_d ≈ 8 GPa, respectively, at P_e= 2 MPa (Fig. 5a), to K_s ≈ 12 GPa and K_d ≈ 13 GPa at P_e≈ 60 MPa, followed by a slight linear increase at higher pressure. Linking the pressure-dependent bulk moduli to the evolution of crack porosity and stiff porosity with increasing pressure allows us to estimate the modification of bulk moduli with changing void space (Fig. 5b). Linear regression of the evolution of dynamic and static bulk moduli of dry samples with respect to crack porosity reveals an average increase of about 1.5 − 2.0 GPa per 0.1% reduction of crack porosity (Fig. 5b). Dynamic bulk modulus is found to be larger than static bulk modulus over an entire range of pressure, and their difference becomes narrow with decreasing crack porosity.

Fig. 5

a The evolution of dynamic bulk modulus K_d and static bulk modulus K_s with increasing pressure during loading for dry and water-saturated samples, respectively. b Dynamic bulk modulus K_d and static bulk modulus K_s as a function of crack porosity and stiff porosity reduction during loading for dry and water-saturated samples, respectively. Error bars indicate the maximum and minimum of two repeated tests. Note that the stiff porosity at pressure of 60 MPa is regarded as a reference value for the calculation of subsequent stiff porosity reduction

Similar to tests performed at dry conditions, static and dynamic bulk moduli of water-saturated samples at drained conditions increase rapidly with increasing effective pressures (Fig. 5a). At pressures above ~ 60 MPa, only a slow linear increase is found. In addition, we observed that the static modulus of water-saturated samples is close to that of dry tests. Static and dynamic bulk moduli measured in drained tests increase significantly as crack porosity decreases, but after the complete closure of compliant cracks, only a modest change is observed (Fig. 5b).

5 Discussion

5.1 Comparison of Measured Crack Porosity with Inverted Crack Porosity from Bulk Moduli

The estimation of evolution of crack porosity with pressure can be theoretically inverted from measured bulk moduli based on some empirical assumptions on void geometry and on the relation between bulk compressibility and external pressure (Zimmerman 1990; Shapiro 2003; David and Zimmerman 2012). Here, we compared the experimentally measured crack porosity in this study with that inverted from the measured static bulk moduli. As suggested by Zimmermann (Zimmerman 1990), static bulk compressibility C_bc (i.e., C_bc = 1/K_s) of porous rocks is frequently fit by exponentially decreasing functions with respect to pressure, expressed by $C_{bc} = C_{bc}^{\infty } + \left( {C_{bc}^{i} - C_{bc}^{\infty } } \right)e^{{{{ - P} \mathord{\left/ {\vphantom {{ - P} {P^{ * } }}} \right. \kern-\nulldelimiterspace} {P^{ * } }}}}$ where $C_{bc}^{i}$ denotes the initial value of bulk compressibility, $C_{bc}^{\infty }$ indicates the value at high pressures, and P^* is a characteristic pressure. In this case, the initial crack porosity $\phi_{c}^{i}$ is analytically derived by $\phi_{c}^{i} = \left( {C_{bc}^{i} - C_{bc}^{\infty } } \right)P^{ * }$ when assuming that smooth-walled thin cracks characterized by a linear closure law in response to pressure are randomly distributed into an elastic and homogeneous solid matrix (Zimmerman 1990). Using the measured static bulk moduli of dry and water-saturated Bentheim sandstone (see Sect. 4.3), the initial crack porosities for dry and water-saturated Bentheim sandstones are estimated to be about 0.11% and 0.12%, respectively. These values are quite similar, but they are smaller than the measured initial crack porosity of about 0.47%. This disagreement may be attributed to the assumption of continuous linear closure of smooth-walled cracks used in the applied model. Experimental observations have shown that asperities are irregularly distributed over grain surfaces in porous sandstones (Bernabe 1991; Fredrich et al. 1995; Wang et al. 2020a), which is expected to modify the closure law of inter-granular cracks. When the asperities come into contact, further closure of inter-granular cracks is increasingly difficult, as opposed to the assumed continuous linear closure law for smooth-walled cracks (Kachanov and Sevostianov 2005; Wang et al. 2020a). Accordingly, the theoretically predicted initial crack porosity likely underestimates the real crack porosity, as supported by our experimental results.

5.2 Static and Dynamic Bulk Moduli of Dry Porous Rocks: Crack Density Inversion

From Fig. 5b, bulk moduli of porous sandstone may depend on the presence of compliant cracks. To evaluate the effects of pore geometry on bulk moduli of dry samples, we used the effective medium theory (Mori–Tanaka scheme) (Mori and Tanaka 1973; Benveniste 1987; David and Zimmerman 2012) to derive crack density from static and dynamic bulk moduli of dry Bentheim sandstone, respectively, by assuming that the elastically isotropic rock contains a population of randomly oriented oblate spheroidal cracks with various aspect ratios (aspect ratio is defined as the ratio of the half length of short axis to the crack radius) in conjunction with stiff spheroidal pores having a fixed aspect ratio (inverted from high-pressure data) generally between 0.01 and 1. The crack density Γ is defined as $\Gamma = {{N\langle a^{3} \rangle } \mathord{\left/ {\vphantom {{N\langle a^{3} \rangle } V}} \right. \kern-\nulldelimiterspace} V}$ where N is the number of oblate spheroidal cracks with radius a in a representative elementary volume V and the angle brackets represent an average. The size of a rosette of strain gages applied for the measurement of quasi-static deformation of Bentheim sandstone samples is about 10 mm in diameter, far larger than the general size of REV [e.g., approximately 7 grain diameters for homogeneous granular materials (Costanza-Robinson et al. 2011)]. The wavelength of ultrasonic waves is estimated to be about ten times larger than the average grain diameter of Bentheim sandstone, suggesting that the measured dynamic bulk modulus represents the mechanical behavior of the whole sample. The specific inversion procedures are given in Appendix B. The aspect ratios of stiff pores deduced from dry static and dynamic bulk moduli at high pressure where only stiff pores remain, are very similar with estimated values of about 0.099 and 0.114, respectively (see Table 3 in Appendix B). In contrast, the crack density estimated from dry static bulk modulus is about 0.6 at P_c = 3 MPa, approximately two times larger than the value inferred from dry dynamic bulk modulus (Fig. 6). Since both static and dynamic bulk moduli are measured on the same rock, the same amount of compliant cracks should be applied for the measurements of static and dynamic bulk moduli. However, crack density estimated from static modulus is higher than estimated from dynamic modulus, and changes more rapidly at low pressures (Fig. 6), as observed in other sandstones (Pervukhina et al. 2010; Zhang et al. 2019). The wave velocity of propagation of an elastic pulse is thought to be less influenced by pre-existing cracks than the deformation of an entire rock sample caused by quasi-static loading (Simmons and Brace 1965; Walsh 1965b). The numerical modelling demonstrates that almost pure elastic deformation occurs at grain contacts during the transmission of ultrasonic waves, whereas quasi-static loading may result in some plastic deformation at grain-grain contacts (Li and Fjær 2012), as observed in our experiments by the occurrence of minor irreversible (plastic) deformation after unloading (Fig. 2).

Fig. 6

Calculated crack densities Γ from experimentally determined dynamic and static bulk moduli of dry Bentheim sandstone samples using effective medium theory (Mori-Tanaka scheme)

Furthermore, together with published experimental data for other porous reservoir sandstones (physical properties and compositions are given in Table 2), the ratios of dynamic bulk moduli to static bulk moduli (K_d/K_s) versus effective pressure and versus crack porosity are displayed in Fig. 7. We clearly found that the ratios between dynamic and static bulk moduli (K_d/K_s) for porous sandstones reduce rapidly from about 1.5−2.0 at ambient pressure to about 1.1 at crack closure pressure under dry conditions and from about 2.0−4.0 to about 1.5 under water-saturated conditions, respectively, followed by almost constant values. Interestingly, the differences in K_d/K_s between different dry porous sandstones are small and tend to vanish with increasing pressure, irrespective of mineral compositions and grain sizes. Likewise, crack closure pressures of about 50−60 MPa for dry samples are similar for the different rock types irrespective of total porosities and grain sizes (Table 2). This suggests that within a range of grain sizes and total porosities, the differences between dynamic and static bulk moduli remain similar for these porous sandstones and are mainly affected by the presence of narrow cracks. It should be emphasized that the crack porosity of other sandstones shown in Fig. 7 was deduced from measured bulk moduli, which may slightly overestimate the influence of cracks (see Sect. 5.1).

Table 2 Physical properties and composition of sandstones in this study.

Fig. 7

a Dynamic-static bulk moduli ratios K_d/K_s as a function of pressure for dry and water-saturated sandstone samples. b Dynamic-static bulk moduli ratios K_d/K_s as functions of crack porosity and stiff porosity reduction for dry and water-saturated sandstone samples. Note that the static and dynamic bulk moduli of Berea, Navajo and Weber sandstone samples were experimentally measured by Coyner (1984) while the corresponding crack porosities were inverted from the measured static bulk moduli by Cheng and Johnston (1981) for Berea sandstone and by Zhang et al. (2019) for Navajo and Weber sandstones

5.3 Effect of Fluid Saturation on Static and Dynamic Bulk Moduli of Porous Rocks

P- wave velocities increase provided that the porous sandstones are water-saturated, which has frequently been reported (King 1966; Gregory 1976; Coyner 1984; Fortin et al. 2007; David and Zimmerman 2012; David et al. 2013). In our experiments, measured W_p of water-saturated Bentheim sandstone is about 10% higher than the dry W_p at initial loading and about 2% larger at high effective pressures when compliant cracks are closed (Fig. 3). In contrast, S-wave velocities for wet Bentheim sandstone are slightly lower than for dry samples (Fig. 3), which has been observed also for other types of sandstones (Nur et al. 1980; Coyner 1984; Winkler 1985; Fortin et al. 2007). The saturated shear modulus (G) is thought to be equal to the dry shear modulus (Gassmann 1951), resulting in lower W_s for saturated rocks compared to dry rocks because of increased density in water-saturated conditions (i.e., $W_{s} = (G/\rho)^{1/2}$). Consequently, increasing W_p and decreasing W_s jointly enlarge the dynamic bulk moduli of water-saturated samples (see Eq. 2) compared to dry samples, as shown in Fig. 5.

When the fluid-saturated porous sandstone is exposed to an oscillating stress field over a wide range of frequency, the deformation of the rock frame might cause fluid pressure variation and further induce fluid flow occurring at different scales: global fluid diffusion at whole rock sample scale or local within a REV (Gurevich et al. 2010; Pimienta et al. 2015). Drained, undrained, and unrelaxed fluid pressure states (O’Connell and Budiansky 1977; Cleary 1978) are expected to occur with increasing loading frequency separated by two characteristic frequencies (Fig. 8). In the drained regime, the induced fluid pressure gradient has sufficient time to equilibrate through global flow between the void structures with different aspect ratios. Thus, the drained bulk modulus of a fluid-saturated rock measured from the stress–strain curve under drained conditions resembles its dry bulk modulus measured at similar loading conditions. In this sense, the static bulk modulus of a fluid-saturated rock is expected to be equivalent to its drained bulk modulus and to its dry static bulk modulus. This is clearly supported by our observation that the static bulk modulus of water-saturated Bentheim sandstone is similar to static bulk modulus of dry samples across the entire range of pressures (Fig. 5a). A slightly lower static bulk modulus compared to dry tests may be due to adsorption of water molecules on silicate mineral surfaces, resulting in a reduction in surface free energy of grains when exposed to water (Tutuncu and Sharma 1992; Pimienta et al. 2014). With increasing loading frequency, fluid exchange only occurs between void structures inside a REV and this state is called the undrained regime, which may be described by Gassmann’s model. In contrast, if the frequency is sufficiently high, like ultrasonic frequency in the laboratory, the fluid pressure may not equilibrate within a REV. In this case, only so-called squirt flow takes place between neighbouring compliant cracks and stiff pores, and local squirt flow may be even inhibited leading to the unrelaxed regime, where each individual void behaves like isolated one. In general, the estimates of the characteristic frequencies f₁ and f₂, which correspond to the drained/undrained transition and the undrained/unrelaxed transition, respectively, are given as (O’Connell and Budiansky 1977; Cleary 1978):

$$f_{1} = \frac{{4\kappa K_{{{\text{drained}}}} }}{{\eta L^{2} }},$$

(3a)

$$f_{2} = \frac{{\alpha^{3} K_{m} }}{\eta },$$

(3b)

where κ is the rock permeability, K_drained is the drained bulk modulus, K_m is the bulk modulus of rock matrix and α is the average aspect ratio of compliant cracks. As suggested by Walsh (1965a), the average aspect ratio of compliant cracks can be estimated from $\alpha = {{4(1 - \nu_{m}^{2} )P_{{{\text{cl}}}} } \mathord{\left/ {\vphantom {{4(1 - \nu_{m}^{2} )P_{{{\text{cl}}}} } {\left[ {3\pi (1 - 2\nu_{m} )K_{m} } \right]}}} \right. \kern-\nulldelimiterspace} {\left[ {3\pi (1 - 2\nu_{m} )K_{m} } \right]}}$ where P_cl is the crack closing pressure and ν_m is the Poisson’s ratio of the rock matrix, respectively. Using P_cl ≈ 60 MPa and measured K_m ≈ 32 GPa and ν_m ≈ 0.085 for Bentheim sandstone (Wang et al. 2020a), α is about 10^–3. Further, taking a permeability κ of about 10^–12 m² and water viscosity η of around 10^–3 Pa·s, the estimates of f₁ and f₂ are about 2 kHz and 36 kHz, respectively. As a consequence, the static bulk modulus (applied frequency ~ 4 × 10^–5 Hz) and the dynamic bulk modulus (ultrasonic wave frequency ~ 1 MHz) of water-saturated Bentheim sandstone measured at low pressures may represent the drained and unrelaxed bulk moduli, respectively. At high pressures, where only stiff pores with an average aspect ratio of about 0.1 remain (see Sect. 5.2), f₂ is estimated to be ~ 36 GHz. This suggests that the dynamic bulk modulus of water-saturated Bentheim sandstone obtained at high pressures may theoretically lie in the undrained regime.

Fig. 8

Schematic illustration of three stress-induced fluid flow regimes associated with loading frequency: drained, undrained and unrelaxed regimes [adapted from O’Connell and Budiansky (1977)]. Each transition between two adjacent regimes corresponds to a critical frequency (f₁ and f₂, respectively). A combination of complaint crack and stiff pore embedded in the representative elementary volume (REV) is considered. The distribution of stress-induced pore pressure within the compliant crack and stiff pore is reflected by the legend in which a darker colour corresponds to a higher pore pressure. The direction of fluid flow at transitional frequencies due to an imposed pore pressure gradient is indicated by the arrows

5.4 The Difference Between Static and Dynamic Bulk Moduli of Fluid-Saturated Porous Rocks: A Micromechanical Model

It is clearly observed that the ratio between dynamic and static bulk moduli of water-saturated sandstones at low pressures is considerably larger than that at high pressures (Fig. 7a), suggesting that the void geometry may affect the difference between the static and dynamic bulk moduli of fluid-saturated porous rocks. To illustrate the influence of void geometry on the difference between dynamic and static bulk moduli of fluid-saturated samples, we consider a micromechanical model in which a hypothetical rock is composed of a solid phase containing randomly oriented fluid-filled spheroidal voids with similar aspect ratio and similar initial pore pressure. To account for the different void geometries from narrow cracks to equant pores, the aspect ratio of spheroidal pores is varied from 10^–3 to 1.

Considering a family of a given aspect ratio α of three-dimensional oblate spheroidal voids (α < 1) are embedded in the elastic solid matrix with Poisson’s ratio ν_m and shear modulus G_m (i.e.,$G_{m} = {{3\left( {1 - 2\nu_{m} } \right)K_{m} } \mathord{\left/ {\vphantom {{3\left( {1 - 2\nu_{m} } \right)K_{m} } {\left[ {2\left( {1 + \nu_{m} } \right)} \right]}}} \right. \kern-\nulldelimiterspace} {\left[ {2\left( {1 + \nu_{m} } \right)} \right]}}$), the compressibility C_pp may be expressed by (Jaeger et al. 2009):

$$C_{{{\text{pp}}}} = \frac{{ - 2\left( {1 - 2\nu_{m} } \right)\left( {1 + 2R} \right) + \left( {1 + 3R} \right)\left[ {1 - 2\left( {1 - 2\nu_{m} } \right)R + 3\alpha^{2} } \right]}}{{4G_{m} \left[ {\left( {1 + 3R} \right)\alpha^{2} + \left( {1 + R} \right)\left( {\nu_{m} + \nu_{m} R + R} \right)} \right]}},$$

(4)

where C_pp is defined as $C_{{{\text{pp}}}} = \left( {{1 \mathord{\left/ {\vphantom {1 {V_{p}^{i} }}} \right. \kern-\nulldelimiterspace} {V_{p}^{i} }}} \right)\left( {{{\partial V_{p} } \mathord{\left/ {\vphantom {{\partial V_{p} } {\partial P_{p} }}} \right. \kern-\nulldelimiterspace} {\partial P_{p} }}} \right)\left| {_{{P_{c} = {\text{constant}}}} } \right.$ and R is given by $R = \frac{ - 1}{{1 - \alpha^{2} }} + \frac{\alpha }{{\left( {1 - \alpha^{2} } \right)^{{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-\nulldelimiterspace} 2}}} }}\arcsin \sqrt {1 - \alpha^{2} }$.

It is necessary to remark that Eq. (4) is derived without considerations of the interactions between oblate spheroidal voids and, thus, may be valid for porous granular sandstone with isotropic and homogeneous void structures (Kachanov 1994). Because the fluid with the same initial pressure is trapped in the same voids resulting in no exchange of fluid flow between them, the pore pressure within these isolated voids is raised equally caused by an increment in external pressure. In other words, these voids seem to be deformed under isolated conditions (i.e., fluid mass in each pore is constant), and thus the corresponding bulk moduli can be treated as undrained (also equivalent to unrelaxed in this case) bulk moduli of fluid-saturated porous rocks. Following Zimmerman’s notation (Zimmerman 2000; Jaeger et al. 2009), the dynamic (undrained and/or unrelaxed) bulk moduli K_d of fluid-saturated porous rock samples, in this case, can be estimated from the poroelasticity as:

$$\frac{1}{{K_{d} }} = \frac{1}{{K_{{{\text{drained}}}} }} - \frac{{C_{bp} C_{pc} }}{{C_{pp} + C_{f} }},$$

(5)

where the compressibilities C_bp and C_pc are defined as $C_{bp} = \left( {{1 \mathord{\left/ {\vphantom {1 {V_{b}^{i} }}} \right. \kern-\nulldelimiterspace} {V_{b}^{i} }}} \right)\left( {{{\partial V_{b} } \mathord{\left/ {\vphantom {{\partial V_{b} } {\partial P_{p} }}} \right. \kern-\nulldelimiterspace} {\partial P_{p} }}} \right)\left| {_{{P_{c} = {\text{constant}}}} } \right.$ and $C_{pc} = - \left( {{1 \mathord{\left/ {\vphantom {1 {V_{p}^{i} }}} \right. \kern-\nulldelimiterspace} {V_{p}^{i} }}} \right)\left( {{{\partial V_{p} } \mathord{\left/ {\vphantom {{\partial V_{p} } {\partial P_{c} }}} \right. \kern-\nulldelimiterspace} {\partial P_{c} }}} \right)\left| {_{{P_{p} = {\text{constant}}}} } \right.$, respectively.

Furthermore, if the rock matrix is microscopically homogeneous and isotropic, the relations between these compressibilities can be expressed using Betti’s reciprocal theorem as (Zimmerman 2000):

$$C_{bp} = \frac{1}{{K_{{{\text{drained}}}} }} - \frac{1}{{K_{m} }},$$

(6a)

$$C_{pp} = C_{pc} - \frac{1}{{K_{m} }},$$

(6b)

$$C_{bp} = \phi C_{pc} .$$

(6c)

Substituting Eq. (4) and Eq. (6) into Eq. (5), dynamic (undrained and/or unrelaxed) and static (drained) bulk moduli of fluid-saturated sample containing various void geometries with different aspect ratios at a given porosity are expected to be acquired.

The calculated ratios between dynamic and static bulk moduli as a function of aspect ratio are shown in Fig. 9a for four different porosities of 0.1%, 0.5%, 10% and 20%. When the voids are mainly compliant cracks with low aspect ratios, the ratio increases strongly with increasing crack porosities (Fig. 9a). In contrast, the discrepancy between dynamic and static bulk moduli is insensitive to stiff pores with α ≥ 0.1. The experimentally observed initial ratio between dynamic and static bulk moduli for water-saturated Bentheim sandstone (initial crack porosity of about 0.47% with an average aspect ratio of 10^–3) is about 3.6 (Fig. 7), which is similar to the modelled value of ~ 3.3 (Fig. 9a). The experimentally observed ratio beyond crack closure pressure is about 1.5, slightly larger than the corresponding modelled value.

Fig. 9

a Predicted ratios of dynamic bulk moduli to static bulk moduli K_d/K_s of water-saturated sandstone samples containing a population of oblate spheroidal pores with aspect ratio varying from 10^–3 to 1 for porosities of ϕ = 20%, ϕ =10%, ϕ = 0.5% and ϕ = 0.1% using a micromechanical model. b Predicted static and dynamic bulk moduli normalized by bulk modulus of rock matrix (K/K_m) as a function of crack porosity with a fixed aspect ratio α=10^–3 or 10⁻², respectively. K_m ≈ 32 GPa and ν_m ≈ 0.085 for Bentheim sandstone (Wang et al. 2020a) and C_f≈ 0.5 GPa⁻¹ were used in above calculation.

The predicted static and dynamic bulk moduli normalized by bulk modulus of rock matrix as a function of crack porosities with a fixed aspect ratio α = 10^–3 or 10^–2 are shown in Fig. 9b. In the case of aspect ratio α = 10^–3, for instance, the normalized static bulk modulus is increased considerably from ~ 0.28 at crack porosity of 0.5% up to unity at crack-free, whereas the similar amount of crack porosity reduction only results in a minor increase in normalized dynamic bulk modulus from ~ 0.93 to 1. Thus, the static bulk modulus of water-saturated porous sandstone increases more rapidly with decreasing crack porosity than corresponding dynamic bulk modulus. This is supported by our experimentally observation that the crack porosity statistically exerts an about 1.6 times stronger effect on static bulk modulus of water-saturated Bentheim sandstone than on corresponding dynamic bulk modulus based on a linear regression analysis (Fig. 5b).

5.5 Comparison of Measured Dynamic Bulk Moduli of Water-Saturated Porous Sandstone with Theoretical Predictions

Based on the measured dry static and dynamic bulk moduli, overall porosity, crack porosity and stiff porosity, we estimated the pressure-dependent dynamic bulk moduli of water-saturated Bentheim sandstone using Gassmann’s equation (Gassmann 1951), the effective medium theory (Mori–Tanaka scheme) (David and Zimmerman 2012) and the squirt-flow model (Gurevich et al. 2010). In Gassmann’s model, the measured static and dynamic bulk moduli of dry samples were used, while only the dry dynamic bulk modulus was used in applying the effective medium model (Mori–Tanaka scheme) and the squirt-flow model. The undrained bulk modulus K_undrained estimated from Gassmann’s equation (Gassmann 1951) is given by:

$$\frac{1}{{K_{{{\text{undrained}}}} }} = \frac{1}{{K_{m} }} + \frac{{\phi \left( {\frac{1}{{K_{f} }} - \frac{1}{{K_{m} }}} \right)}}{{1 + {{\phi \left( {\frac{1}{{K_{f} }} - \frac{1}{{K_{m} }}} \right)} \mathord{\left/ {\vphantom {{\phi \left( {\frac{1}{{K_{f} }} - \frac{1}{{K_{m} }}} \right)} {\left( {\frac{1}{{K_{{{\text{dry}}}} }} - \frac{1}{{K_{m} }}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {\frac{1}{{K_{{{\text{dry}}}} }} - \frac{1}{{K_{m} }}} \right)}}}},$$

(7)

where K_dry indicates either static or dynamic bulk modulus of dry samples. The details of numerical calculation using the effective medium theory (Mori–Tanaka scheme) and squirt-flow model for fluid-saturated porous media are presented in Appendix C and Appendix D, respectively.

As shown in Fig. 10a, the undrained bulk moduli predicted by Gassmann’s equation using dry static and dynamic bulk moduli are both lower than the measured dynamic bulk moduli of water-saturated Bentheim sandstone samples across the entire range of applied pressure. In contrast, the predictions by the effective medium model (Mori–Tanaka scheme) and by the squirt-flow model both overestimate the dynamic bulk moduli of water-saturated Bentheim sandstone at effective pressures below 15 MPa, but reveal an underestimation at higher pressures (Fig. 10a). Obviously, the unrelaxed bulk moduli predicted by the effective medium theory (Mori–Tanaka scheme) and by the squirt-flow model are similar across the entire range of effective pressures up to 180 MPa. The corresponding predicted unrelaxed bulk moduli using both models are found to be insensitive to applied pressure with an increase of about 5% from low to high pressure. This indicates that the trapped fluid within the compliant cracks has no time to distribute into the surrounding stiff pores, mimicking the behavior of closed cracks (i.e., unrelaxed state).

Fig. 10

a Measured and predicted dynamic bulk moduli as a function of effective pressure for water-saturated Bentheim sandstone samples. b Measured and predicted ratios of dynamic bulk moduli to static bulk moduli of water-saturated Bentheim sandstone samples (K_d/K_s) as a function of pressure. Note that the measured static bulk modulus of water-saturated Bentheim sandstone samples was used to calculate K_d/K_s

At low effective pressures (P_e < 10 MPa), the measured dynamic bulk moduli of water-saturated Bentheim sandstone samples are between the undrained bulk moduli predicted by Gassmann’s equation and the unrelaxed bulk moduli predicted by the effective medium model (Mori–Tanaka scheme) and the squirt-flow model. Similar observations were made for other water-saturated sandstone samples (David and Zimmerman 2012). This suggests that the dynamic bulk moduli of water-saturated Bentheim sandstone samples may represent a partially relaxed state. That is, at low pressures, the fluid communication between compliant cracks and adjacent stiff pores may occur locally at ultrasonic frequency, possibly due to the complex void structures present in Bentheim sandstone (see Fig. 1). Therefore, the real inter-granular morphology of grain-grain contacts (e.g., rough surfaces of grains) may not be oversimplified by the parallel disc-shaped gaps assumed in the squirt-flow model or by the smooth penny-shaped cracks used in the effective medium theory (Mori–Tanaka scheme).

At effective pressures above 10 MPa, the measured dynamic bulk moduli are higher than the corresponding predictions using Gassmann’s equation, effective medium model (Mori–Tanaka scheme) and squirt-flow model (Fig. 10a). The cause for this difference between measured bulk moduli and the predicted values at high pressures is not clear yet. One possible reason is an oversimplification by assuming that all stiff pores have similar aspect ratio (e.g., effective medium model) (de Paula et al. 2012; Sun et al. 2019). Another potential reason for this discrepancy may arise from the implicit assumption of K_ϕ = K_m [K_ϕ is the bulk modulus of pore space, defined by $K_{\phi}=-(V_{P}^{i})(\partial P_{c}/\partial V_{P})\left|_{P_{c}-P_{P}={\text{constant}}}\right.$ (Brown and Korringa 1975; Hart and Wang 2010)] made in Gassmann’s equation, the squirt-flow model and in the effective medium model (Mori–Tanaka scheme). Previous studies reveal that K_ϕ may be related to the pore structures of rocks (Brown and Korringa 1975; Hart and Wang 1995). For instance, when more compliant grains are located in the pore space and are not tightly bound or cemented to the rock matrix (Hart and Wang 2010), this assumption might be invalid. K_ϕ has been experimentally found to be pressure-dependent for several types of sandstones (Brown and Korringa 1975; Hart and Wang 1995; Tarokh et al. 2018) and sometimes even negative (Berge and Berryman 1995). In general, undrained bulk modulus K_undrained can be independently obtained from the poroelastic theory, as given by (Brown and Korringa 1975; Jaeger et al. 2009):

$$\frac{1}{{K_{{{\text{undrained}}}} }} = \frac{1}{{K_{m} }} + \frac{{\phi \left( {\frac{1}{{K_{f} }} - \frac{1}{{K_{\phi } }}} \right)}}{{1 + {{\phi \left( {\frac{1}{{K_{f} }} - \frac{1}{{K_{\phi } }}} \right)} \mathord{\left/ {\vphantom {{\phi \left( {\frac{1}{{K_{f} }} - \frac{1}{{K_{\phi } }}} \right)} {\left( {\frac{1}{{K_{{{\text{dry}}}} }} - \frac{1}{{K_{m} }}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {\frac{1}{{K_{{{\text{dry}}}} }} - \frac{1}{{K_{m} }}} \right)}}}}.$$

(8)

In the case of K_ϕ = K_m, Eq. (8) is reduced to Gassmann’s equation (see Eq. 7). If K_undrained and K_dry in Eq. (8) are substituted by the experimentally measured dynamic bulk modulus of water-saturated Bentheim sandstone samples and the dry dynamic bulk modulus, respectively, K_ϕ is found to be in the range of ~ 3.3–5.0 GPa at P_e > 10 MPa. This is in good agreement with the range reported from direct laboratory measurements for Bentheim sandstone (Blöcher et al. 2014).

As a whole, the predictions of dynamic bulk modulus of water-saturated Bentheim sandstone using the effective medium theory (Mori–Tanaka scheme) and the squirt-flow model roughly match the experimental data at P_e > 5 MPa (Fig. 10).

6 Conclusion

We examined experimentally the continuous evolution of static and dynamic bulk moduli for dry and water-saturated porous Bentheim sandstones deformed at hydrostatic pressure conditions. Static bulk moduli were calculated from pressure-volumetric strain curves and dynamic bulk moduli were obtained from ultrasonic wave velocities measurements (~ 1 MHz). The entire void space present in porous sandstones, to the first order, is divided into equant (stiff) pores and narrow (compliant) cracks. The evolution of crack porosity and stiff porosity with pressure was experimentally estimated from the measured deformation of void space. We found that dynamic and static bulk moduli for dry and water-saturated Bentheim sandstones samples increase considerably until the applied pressure leads to closure of the narrow cracks, and afterwards only a minor increase was observed. The deduced crack porosity from bulk moduli using effective medium models is found to be lower than the real crack porosity, possibly due to the fact that the continuous linear closure law for narrow cracks in response to pressure is assumed. Taking into account also published data of other porous sandstones (Berea, Navajo and Weber sandstones), the ratios between dynamic and static bulk moduli (K_d/K_s) for porous sandstones reduce rapidly from about 1.5−2.0 at ambient pressure to about 1.1 at crack closure pressure under dry conditions and from about 2.0−4.0 to about 1.5 under water-saturated conditions, respectively. The pressure-dependent bulk moduli are explained by the pressure-dependent evolution of void structures, revealing that K_d/K_s is positively correlated with the amount of narrow cracks, but remains almost unaffected by the presence of equant pores. The difference between dynamic and static bulk moduli is higher for fluid-saturated porous rocks, compared to dry samples. This may be attributed to the high pore pressure locally maintained in void space during the measurement of ultrasonic wave velocities. In our experiments, the pressure dependence of dynamic bulk modulus of water-saturated Bentheim sandstone at effective pressures above 5 MPa can be roughly predicted by both the effective medium theory (Mori–Tanaka scheme) and the squirt-flow model. Static bulk moduli are found to be more sensitive to narrow cracks than dynamic bulk moduli for porous sandstones under dry and water-saturated conditions.

Appendices

Appendices

1.1 Appendix A: Experimental Setup

See Fig. 11.

Fig. 11

a Sensor map exhibiting positions of P-wave PZTs (P1 to P12), S-wave PZTs (S_H1, S_H2, S_V1 and S_V2) and two pairs of strain gages. Half of the sensors (marked by orange labels) were used for ultrasonic transmission. b Photograph of a sample assembly with mounted sensors and Neoprene jacket. c Sketch of the experimental apparatus (Wang et al. 2020a)

1.2 Appendix B: Inversion of Pressure-Dependent Bulk Moduli of Dry Porous Rocks for Crack Density

The crack density of dry Bentheim sandstone with pressure is inverted using effective medium theory (Mori–Tanaka scheme) (Mori and Tanaka 1973; Benveniste 1987; David and Zimmerman 2012) which considers a random distribution of oblate spheroidal cracks with varied aspect ratio and a family of stiff pores having the same aspect ratio (inverted from high-pressure experimental data) are embedded in an isotropic solid matrix with bulk modulus K_m and Poisson’s ratio ν_m. The inversion procedures include two main steps.

Step 1: Calculating the aspect ratio of stiff, non-closable pores α^hp by inverting the high-pressure static and dynamic bulk compressibility of dry rock samples after compliant cracks are closed. Since the stiff pores hardly change with pressure, the compliant cracks can be considered to be introduced into a host material composed of mineral phase plus the stiff, non-closable pores. Specifically, according to the Mori–Tanaka method, the effective moduli of dry porous solids containing the stiff spheroidal pores having the same aspect ratio α (0.01 ≤ α ≤ 1) are explicitly described as (Benveniste 1987):

$$\frac{{K_{m} }}{K} = 1 + \frac{\phi }{1 - \phi }P\left( {\alpha ,\nu_{m} } \right),$$

(B.1)

where ϕ is the porosity (taken as stiff porosity at high pressure above crack closure pressure), and P is explicitly written as functions of ν_m and α. For the aspect ratio of spheroid pore α ≤ 0.3, P can be approximated by (David and Zimmerman 2011):

$$P \approx \frac{{P_{ - 1} }}{\alpha } + P_{0} + P_{1} \alpha ,$$

(B.2)

where P_–1, P₀ and P₁ are, respectively, given by:

$$P_{ - 1} = \frac{{4\left( {1 - \nu_{m}^{2} } \right)}}{{3\pi \left( {1 - 2\nu_{m} } \right)}},$$

(B.3)

$$P_{0} = \frac{1}{6}\left( {1 - \nu_{m} } \right)\left( {1 - 2\nu_{m} } \right),$$

(B.4)

$$P_{1} = \frac{{1 - \nu_{m}^{2} }}{{12\left( {1 - 2\nu_{m} } \right)}}\left[ {\pi \left( {1 - 2\nu_{m} } \right)^{2} + \frac{{8\left( {7 - 8\nu_{m} } \right)}}{\pi }} \right].$$

(B.5)

Thus, the parameter α^hp would be found by a least-square regression for the experimentally measured data of dynamic and static bulk moduli at high pressure beyond 70 MPa according to the above equations. The K_m and ν_m for Bentheim sandstone are taken as ~ 32 GPa and ~ 0.085, respectively (Wang et al. 2020a). The ϕ is measured from the drained compression tests.

Clearly, the inverted parameters α^hp from the high-pressure data of static bulk moduli of Bentheim sandstone are consistent with that from dynamic bulk moduli (Table 3). Furthermore, the values of high-pressure static and dynamic bulk moduli obtained from such inversion are referred to as K_s^hp and K_d^hp, respectively, both of which would be used in Step 2.

Table 3 Inversion results for α^hp from dry static and dynamic bulk moduli of Bentheim sandstones at high pressure beyond 70 MPa

Step 2: Crack density Γ is inferred from the experimental measurements of static or dynamic bulk moduli of dry rocks at low pressure according to the following equation (Jaeger et al. 2009; David and Zimmerman 2012):

$$\Gamma = \frac{{9\left( {1 - 2\nu^{{{\text{hp}}}} } \right)\left( {K^{{{\text{hp}}}} - K} \right)}}{{16\left[ {1 - \left( {\nu^{{{\text{hp}}}} } \right)^{2} } \right]K}},$$

(B.6)

where the crack density Γ is defined as $\Gamma = {{N\langle a^{3} \rangle } \mathord{\left/ {\vphantom {{N\langle a^{3} \rangle } V}} \right. \kern-\nulldelimiterspace} V}$ in which N is the number of oblate spheroidal cracks with radius a in a representative elementary volume V, and the angle brackets represents an average. Note that K^hp is taken as the inverted K_s^hp or K_d^hp in Step 1, and meanwhile K is given by the corresponding experimentally determined static or dynamic bulk moduli. For simplicity, the high-pressure Poisson’s ratio ν^hp is inferred from the ultrasonic wave velocity of dry rocks at high pressures using $\nu^{hp} = {{0.5\left( {\gamma^{2} - 2} \right)} \mathord{\left/ {\vphantom {{0.5\left( {\gamma^{2} - 2} \right)} {\left( {\gamma^{2} - 1} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {\gamma^{2} - 1} \right)}}$ where γ is the ratio of P-wave to S-wave velocities.

1.3 Appendix C: Prediction of Unrelaxed Bulk Moduli Using Effective Medium Theory (Mori–Tanaka Scheme)

Unrelaxed bulk moduli of fluid-saturated porous materials can be estimated from effective medium theory when assuming randomly oriented oblate spheroidal cracks with various aspect ratios and a family of stiff pores having same aspect ratio are embedded in an isotropic solid matrix. The fluid-filled compliant cracks and stiff pores are implicitly treated to be completely isolated with regards to fluid flow. The detailed prediction involves two main steps.

Step 1: Calculating the high-pressure unrelaxed bulk moduli $K_{{{\text{unrelaxed}}}}^{{{\text{hp}}}}$ of rocks composed of minerals and saturated non-closable pores. Based on the fact that the pore structure of fully fluid-saturated sample should be same as for a dry tests at the same effective pressure, the unrelaxed bulk moduli $K_{{{\text{unrelaxed}}}}^{{{\text{hp}}}}$ can be obtained from α^hp (inverted from the experimental data of dry dynamic bulk moduli, see Appendix B), total porosity ϕ (stiff porosity in drained tests at high pressure), elastic parameters of rock matrix (K_m and ν_m) and the ratio of solid compressibility and fluid compressibility (ξ) according to the following equation (Benveniste 1987):

$$\frac{{K_{m} }}{{K_{{{\text{unrelaxed}}}} }} = 1 + \frac{\phi }{1 - \phi }\left[ {\frac{{P_{u} \left( {1 - \xi } \right)}}{{1 + \frac{\phi }{1 - \phi }P_{u} \xi }}} \right],$$

(C.1)

where K_unrelaxed is the unrelaxed bulk modulus, P_u is undrained pore compliance coefficient which can be expressed as functions of spheroid’s ratio α, the solid’s Poisson’s ratio ν_m and the ratio of solid and fluid compressiblities ξ (ξ = K_f/K_m). Specifically, P_u is theoretically related to P (i.e., in Eq. (B.2) for the dry rock) by:

$$P_{u} = \frac{P}{1 - \xi }\frac{\delta }{1 + \delta },$$

(C.2)

where the parameter δ characterises the coupling between pore fluid and solid pressures which can be expressed as (Shafiro and Kachanov 1997):

$$\delta = \frac{1 - \xi }{{P\xi }},$$

(C.3)

Because P is correlated with α and ν_m, P_u fundamentally relies on aspect ratio α, ν_m and ξ.

Step 2: Predicting the pressure dependence of unrelaxed bulk moduli prior to crack closure pressure. Unlike the bulk moduli predicted in the dry tests (see Eq. B.6), the expression of unrelaxed bulk moduli depend not only on crack density but also on the aspect ratio of the cracks.

Remembering the compliant crack are distributed in host materials containing minerals and stiff pore with saturated fluid prior to crack closure pressure, the unrelaxed bulk compressibility in Eq. (C.1) can be rewritten as:

$$\frac{{K_{{{\text{unrelaxed}}}}^{{{\text{hp}}}} }}{{K_{{{\text{unrelaxed}}}} (P_{e} )}} = 1 + \frac{{\phi_{c} }}{{1 - \phi_{c} }}\left[ {\frac{{P_{u} \left( {\alpha (P_{e} ),{{K_{f} } \mathord{\left/ {\vphantom {{K_{f} } {K_{{{\text{unrelaxed}}}}^{{{\text{hp}}}} }}} \right. \kern-\nulldelimiterspace} {K_{{{\text{unrelaxed}}}}^{{{\text{hp}}}} }},\nu_{{{\text{unrelaxed}}}}^{{{\text{hp}}}} } \right) \cdot \left( {1 - \xi } \right)}}{{1 + \frac{{\phi_{c} }}{{1 - \phi_{c} }}P_{u} \xi }}} \right].$$

(C.4)

Note that K_m, ϕ, α, ν_m and ξ (= K_f/K_m) in Eq. (C.1) should be replaced by $K_{{{\text{unrelaxed}}}}^{{{\text{hp}}}}$ (inverted unrelaxed bulk modulus at high pressure), ϕ_c (crack porosity), α (P_e) (average aspect ratio of cracks), $\nu_{{{\text{unrelaxed}}}}^{{{\text{hp}}}}$ (unrelaxed Poisson’s ratio at high pressure) and ξ ($= {{K_{f} } \mathord{\left/ {\vphantom {{K_{f} } {\overline{K}^{hp} }}} \right. \kern-\nulldelimiterspace} {\overline{K}^{hp} }}$), respectively. For simplicity, $\nu_{{{\text{unrelaxed}}}}^{{{\text{hp}}}}$ is inferred from the ultrasonic wave velocity of fluid-saturated rocks at high pressures. Moreover, at a given effective pressure P_e, the average aspect ratio of crack α (P_e) can be estimated by $\alpha (P_{e} ) = {{\phi_{c} } \mathord{\left/ {\vphantom {{\phi_{c} } {\left( {{{4\pi \Gamma } \mathord{\left/ {\vphantom {{4\pi \Gamma } 3}} \right. \kern-\nulldelimiterspace} 3}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {{{4\pi \Gamma } \mathord{\left/ {\vphantom {{4\pi \Gamma } 3}} \right. \kern-\nulldelimiterspace} 3}} \right)}}$ where Γ is inverted from dry dynamic bulk moduli (see Eq. (B.6)) and ϕ_cis obtained from drained tests. For additional details, refer to David and Zimmerman (2012).

1.4 Appendix D: Prediction of Unrelaxed Bulk Moduli Using Squirt-flow Model

The squirt-flow provided by Gurevich et al. (2010) considers the local flow between compliant cracks and stiff pores induced by high-frequency wave oscillations. The assumption of parallel disc-shaped gap mimicking the inter-granular crack and toroidal pores approximating the stiff pores surrounded by neighbouring grains is made in squirt-flow model. The frequency-dependent bulk modulus of fluid-saturated rocks K_ur (K_ur is expected to be K_unrelaxed if applied frequency is larger than f₂, see Fig. 8) can be estimated from modified frame modulus K_mf (whereby soft pores are fluid-filled while stiff pores are dry) using Gassmann’s equations as follows (Gurevich et al. 2010):

$$\frac{1}{{K_{{{\text{ur}}}} (P_{e} ,\omega )}} = \frac{1}{{K_{m} }} + \frac{{\phi_{s} (P_{e} )\left( {\frac{1}{{K_{f} }} - \frac{1}{{K_{m} }}} \right)}}{{1 + {{\phi_{s} (P_{e} )\left( {\frac{1}{{K_{f} }} - \frac{1}{{K_{m} }}} \right)} \mathord{\left/ {\vphantom {{\phi_{s} (P_{e} )\left( {\frac{1}{{K_{f} }} - \frac{1}{{K_{m} }}} \right)} {\left( {\frac{1}{{K_{mf} (P_{e} ,\omega )}} - \frac{1}{{K_{m} }}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {\frac{1}{{K_{mf} (P_{e} ,\omega )}} - \frac{1}{{K_{m} }}} \right)}}}},$$

(D.1)

where K_m and K_f are the bulk moduli of the rock matrix and the pore fluid, respectively, ω is the angular frequency (ω = 2πf, f is the frequency), ϕ_s is the stiff porosity. The modified frame modulus K_mf interpreted by adding inclusions of compliant crack porosity (ϕ_c) to the background phase (K_drs) consisting of rock matrix (K_m) and dry stiff porosity (ϕ_s) is estimated by:

$$\frac{1}{{K_{mf} (P_{e} ,\omega )}} = \frac{1}{{K_{{{\text{drs}}}} }} + \frac{1}{{\frac{1}{{\frac{1}{{K_{{{\text{dry}}}} (P_{e} )}} - \frac{1}{{K_{{{\text{drs}}}} }}}} + \frac{1}{{\left( {\frac{1}{{K_{f}^{ * } (P_{e} ,\omega )}} - \frac{1}{{K_{m} }}} \right)\phi_{c} (P_{e} )}}}},$$

(D.2)

where K_drs is taken as the high-pressure dry bulk modulus measured from ultrasonic wave velocities after crack closure (here we adopted the experimental data at maximum confining pressure of 180 MPa), K_dry (P_e) is experimentally derived dry bulk modulus from ultrasonic wave velocities with pressure. The frequency-dependent fluid modulus K_f^* related to the mean aspect ratio of pores and the fluid viscosity η, can be expressed by:

$$K_{f}^{ * } (P_{e} ,\omega ) = \left[ {1 - \frac{{2J_{1} (\beta )}}{{\beta J_{0} (\beta )}}} \right]K_{f} ,$$

(D.3)

where $\beta = {{\left( {{{ - 3i\omega \eta } \mathord{\left/ {\vphantom {{ - 3i\omega \eta } {K_{f} }}} \right. \kern-\nulldelimiterspace} {K_{f} }}} \right)^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}}} } \mathord{\left/ {\vphantom {{\left( {{{ - 3i\omega \eta } \mathord{\left/ {\vphantom {{ - 3i\omega \eta } {K_{f} }}} \right. \kern-\nulldelimiterspace} {K_{f} }}} \right)^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}}} } \alpha }} \right. \kern-\nulldelimiterspace} \alpha }$, α is assumed to be constant at different pressure (≈10^–3 in this study), and J_n is the Bessel function of nth order. The required ϕ_s and ϕ_c are directly taken from measurements in the drained tests. Eq. (D.1) is expected to reduce to the Gassmann’s equation if K_mf is replaced by drained bulk moduli at very low frequency.