Rock Physics: Recent History and Advances

Abstract

This chapter presents the basics of rock physics, the science exploring quantitative relations between various properties (attributes) of the holistic object we call natural rock. This chapter includes several sections, starting with the history and basics; proceeding to the effects of the pore fluid on rock properties; discussing several variables that influence the elastic properties of rocks; presenting selected theories that relate the elastic properties to the porosity, mineralogy, and texture of rocks; and introducing the latest development, digital rock physics. Data examples shown here illustrate qualitative reasoning. Equations are presented as well to mathematically express the conceptual theories discussed. Most importantly, rock physics references are listed to help the reader become willing to delve deeper into the topic and start applying rock physics theories, concepts, and ideas to field data.

1. Introduction: subject of rock physics, background, and brief history

Rock physics is often called a “velocity-porosity” science. The idea behind this name is to predict the elastic-wave velocities in porous rock from its porosity or implement an inverse operation and interpret the velocity measured in a well or using seismic tomography or reflection techniques for the porosity of rock. It is important to mention that the elastic-wave velocities are related to the elastic moduli of rock as follows:

$$V_p=\sqrt{\frac{K+\mathrm{4/3}G}{{\rho }_b}};V_s=\sqrt{\frac{G}{{\rho }_b}},$$E1

where $V_p$ and $V_s$ are the P- and S-wave velocities, respectively; $K$ and $G$ are the bulk and shear moduli, respectively; and ${\rho }_b$ is the bulk density. The latter quantity is related to the total porosity $\varphi$ as

$${\rho }_b=\left(1-\varphi \right){\rho }_s+{\mathrm{\varphi \rho }}_f,$$E2

where ${\rho }_s$ is the density of the mineral matrix also called the solid component of the rock, while ${\rho }_f$ is the density of the pore fluid.

Important elastic constants used in rock physics are the bulk ($K$), shear ($G$), and compressional ($M$) moduli, as well as the P-wave ($I_p$) and S-wave ($I_s$) impedances and Poisson’s ratio ($\nu$):

$$\begin{array}{l}M={\rho }_b{V}_p^2;G={\rho }_b{V}_s^2;K=M-\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.G;\\ I_p={\rho }_b{V}_p;I_s={\rho }_b{V}_s;\nu =\frac{1}{2}\frac{{\left(V_p/V_s\right)}^2-2}{{\left(V_p/V_s\right)}^2-1}.\end{array}$$E3

Most of natural rocks contain more than one mineral. In this situation, ${\rho }_s$ can be computed as the arithmetic average of the densities of the individual components:

$${\rho }_s=\sum _{i=1}^N{f}_i{\rho }_i,$$E4

where $f_i$ is the volume fraction of the i-th mineral component in the mineral matrix and ${\rho }_i$ is its density. These individual densities can be found in handbooks, such as Mavko et al. [1]. They can vary between, e.g., 2.58 g/cc in clay and 4.93 g/cc in pyrite.

The same rule applies to the density of the pore fluid:

$${\rho }_f=S_w{\rho }_w+S_o{\rho }_o+S_g{\rho }_g,$$E5

where $S_w$, $S_o$, and $S_g$ are the water, oil, and gas saturations in the pore space, respectively, and ${\rho }_w$, ${\rho }_o$, and ${\rho }_g$ are the densities of these pore fluid components. Of course, it is required that

$$\sum _{i=1}^N{f}_i=1$$E6

and

$$S_w+S_o+S_g=1.$$E7

Because of the link between the elastic-wave velocities and elastic moduli as given by Eq. (1), it is often instructive to relate these elastic moduli to porosity. Such approach opens an avenue to using the so-called effective medium theories where the elastic moduli are theoretically related to porosity and the geometry of rock, referring to the spatial arrangement of pores and grains, as well as shapes of these pores and grains.

It has been discovered early that the velocity and elastic moduli not only depend on porosity, but also on the properties of the mineral frame. A rule of thumb is that at the same porosity, the softer the mineral frame, the smaller the elastic moduli of rock. For example, at the same porosity, rocks containing soft clays have velocities smaller than rocks dominated by stiffer quartz. Hence, rock physics is not only a “velocity-porosity” science but also a “velocity-porosity-mineralogy” science.

The situation becomes more complex if we consider the effects of the pore fluid on the elastic moduli (and velocities) of a porous composite. It is intuitively clear that the less compressible the pore fluid (water versus gas), the stiffer the entire rock, meaning that its bulk modulus is higher. Now we are talking about “velocity-porosity-mineralogy-fluid.”

The science of rock physics also includes understanding and quantification of other rock properties, such as hydraulic permeability and electrical resistivity, and their relation to other attributes, namely, porosity, rock texture, and mineralogy.

Generally, contemporary rock physics treats natural rock as a holistic object whose various properties (attributes) are extracted from experiments simulating processes, such as elastic-wave propagation, fluid and electrical transport, nuclear magnetic resonance (NMR), and breakage. We seek a theoretical understanding of interrelations between such attributes and their mathematical quantification. Such relations are also called rock physics models (RPM) or transforms. Needless to say that such quantification has to be “as simple as possible but not simpler.”

Finally, the newest branch of rock physics is digital rock physics (DRP) whose mandate is to “image and compute,” image rock at the pore scale and digitally simulate various processes within the digital image. For example, simulations of viscous fluid flow yield permeability, simulations of electrical charge transport yield resistivity, and simulations of deformation under stress yield the elastic moduli.

Let us now review some of historic developments in rock physics.

Arguably, the first rock physics velocity-porosity transform was introduced by Wyllie et al. [2]. It simply states that the total P-wave traveltime through rock with porosity $\varphi$ is the sum of the travel times through the mineral and fluid parts of the rock. This is why it is called the time-average equation. In terms of the P-wave velocities, this formulation is

$$\frac{1}{V_p}=\frac{1-\varphi }{V_{\mathrm{ps}}}+\frac{\varphi }{V_{\mathrm{pf}}},$$E8

where $V_p$ is the P-wave velocity, $V_{\mathrm{ps}}$ is the velocity in the mineral phase, and $V_{\mathrm{pf}}$ is that in the fluid phase. Examples for 100% quartz and 100% dolomite rock are shown in Figure 1. Also shown is an example for rock with mixed 50% quartz and 50% dolomite mineralogy. At the same porosity, $V_p$ is highest in stiffer dolomite, lowest in softer quartz, and falls in between for the mixed mineralogy. The pore fluid was water with $V_{\mathrm{pf}}$ = 1500 m/s.

Figure 1.

V_p versus porosity according to the Wyllie et al. [2] and Raymer et al. [3] transforms for quartz, dolomite, and mixed mineralogy. Legend in the middle refers to all plots.

Equation (8) is purely empirical in spite of its physically meaningful form. Indeed, in real rock, the mineral and fluid parts are not arranged in layers to enable a simple summation of the respective traveltimes. Still, this equation gives a reasonably accurate approximation for $V_p$ in “fast” sediments as discussed in Mavko et al. [1]. Also note that it can only work for rock with liquid since in vacuum dry rock, $V_{\mathrm{pf}}$ = 0. Yet, as have been shown by seismic experiments on the moon, $V_p$ in such sediment is finite.

Equation (8) has dominated petrophysical interpretation of velocity for porosity for a long time. It gave rise to the so-called sonic porosity computed from wireline velocity data as

$$\varphi =\frac{V_p^{-1}-V_{\mathrm{ps}}^{-1}}{V_{\mathrm{pf}}^{-1}-V_{\mathrm{ps}}^{-1}}.$$E9

The next historic equation was introduced by Raymer et al. [3]:

$$V_p={\left(1-\varphi \right)}^2{V}_{\mathit{ps}}+\varphi V_{\mathit{pf}}.$$E10

As Eq. (8), it is purely empirical, derived from wireline data. Still, it is very meaningful as it can be applied to rock with any fluid inside, even where $V_{\mathrm{pf}}$ = 0. As shown in Dvorkin et al. [4], it is more accurate than the Wyllie et al. [2] time average if applied to “fast” consolidated sediments. Velocity-porosity examples according to this equation are also shown in Figure 1.

We conclude this section by presenting equations relating the electrical resistivity to porosity and absolute hydraulic permeability to porosity.

The former transform relates the resistivity $R_t$ of rock fully saturated with conductive fluid (brine) with resistivity $R_w$ as

$$F=\frac{R_t}{R_w}=\frac{1}{{\varphi }^m},$$E11

where $F$ is called the formation factor and $m$ is the cementation exponent. In many sandstones $m$ is approximately 2; however it may be much larger in carbonates [1]. Figure 2 shows experimental data for Fontainebleau sandstone [5] with Eq. (11) curves for $m$ = 1.5, 2.0, and 2.5 superimposed.

Figure 2.

Left: F versus porosity according to Eq. (11) for m = 1.5, 2.0, and 2.5 with Fontainebleau experimental data shown as symbols. Right: R_tS/R_w ratio versus water saturation for ϕ = 0.2 and m = n = 2.0 (Eq. (12)).

At partial brine saturation, $S_w$ < 1, the resistivity of rock $R_{\mathrm{tS}}$ not only depends on porosity but also on saturation $S_w$ as

$$\frac{R_{\mathrm{tS}}}{R_w}=\frac{1}{{\varphi }^m{S}_w^n}=\frac{F}{S_w^n},$$E12

where $n$ is the saturation exponent. This exponent is much more elusive than $m$ since laboratory experiments measuring resistivity at partial saturation are scarce. Generally, $n$ should be larger than 1.0 and approach 2.0. An example of $R_{\mathrm{tS}}/R_w$ versus $S_w$ is shown in Figure 2 for porosity 0.2, $m$ = 2.0, and $n$ = 2.0.

Both Eqs. (11) and (12) were discovered by Archie in 1942 [6] and remain the cornerstone of resistivity interpretation for hydrocarbon saturation in the wellbore. Various modifications of these equations dealing with resistivity interpretation in sediments containing clays and shales are discussed in Mavko et al. [1].

The historic absolute permeability prediction equation is called the Kozeny-Carman [7] formula. It is based on an extremely idealized representation of pores as a set of parallel pipes inclined to the direction of pore pressure gradient at an angle $\alpha$. The tortuosity $\tau$ of these pores is defined as

$$\tau =1/cos\alpha \ge 1.$$E13

The permeability $k$ is also a function of the specific surface area $S$ defined as the ratio of surface of the pore space $S_{\text{Pore}}$ to the total volume $V$ of the rock sample:

$$S=S_{\text{Pore}}/V.$$E14

A variable alternative to $S$ is the grain size (or grain diameter) $d$.

The Kozeny-Carman equation reads [1]

$$k=\frac{1}{2}\frac{{\varphi }^3}{S^2{\tau }^2}=\frac{1}{72}{d}^2\frac{{\varphi }^3}{{\left(1-\varphi \right)}^2{\tau }^2}.$$E15

A modified version of this equation is based on the assumption that $k$ becomes zero not at zero porosity but at a finite and very small porosity value ${\varphi }_p$ called the percolation porosity:

$$k=\frac{1}{72}{d}^2\frac{{\left(\varphi -{\varphi }_p\right)}^3}{{\left[1-\left(\varphi -{\varphi }_p\right)\right]}^2{\tau }^2}.$$E16

It follows from Eq. (15) that the unit of absolute permeability is length squared. However, traditionally, the permeability unit is Darcy (D) or milli-Darcy (mD). One D is 10⁻¹³ m², while one mD = 10⁻¹⁵ m².

Figure 3 shows experimental permeability data for Fountainebleau sandstone and two North Sea sand sets with an Eq. (16) curve superimposed for $d$ = 0.25 mm, $\tau$ = 2.5, and ${\varphi }_p$ = 0.02. This theoretical curve matches the Fountainebleau data, while the permeability from the other two datasets falls below this curves. The reason is the varying grain size as discussed in Mavko et al. [1].

Figure 3.

Permeability versus porosity plots as explained in the text.

2. Effect of pore fluid on elastic properties

Laboratory experiments measuring the elastic-wave velocities in rock often show that the presence of the fluid in the pores strongly affects the elastic properties (Figure 4). Such dramatic results, especially for $V_p$, are in part due to the fact that such experiments are commonly conducted at very high frequencies, on the order of 1 MHz. In this frequency range, the fluid in the pores is “unrelaxed” and acts to strongly reinforce the soft mineral frame, thus increasing the bulk modulus (e.g., [1]).

Figure 4.

V_p (left) and V_s (right) of high-porosity unconsolidated sand versus hydrostatic confining pressure. The pore pressure is constant 0.1 MPa. Squares are data obtained in ultrasonic pulse transmission experiments on the water-saturated sample. Circles are for the room-dry sample (after Zimmer [8]).

Arguably, the most important contribution to rock physics is Gassmann’s fluid substitution theory [9]. This theory allows us to compute the bulk modulus of porous rock filled with Fluid A if this modulus is known (measured) in the same rock but filled with Fluid B. These derivations were conducted under the assumption that the wave-induced pore pressure oscillations equilibrate within the sample over the wave period, meaning that Gassmann’s is a low-frequency theory. Hence, it is applicable at the wireline and seismic frequency ranges. It helps predict the seismic response of rock filled with any hypothetical fluid if it is measured in situ where the pore fluid is known. For example, if the elastic properties of rock are measured in situ in rock 100% filled with water, we can predict these properties in the same rock but filled with oil or gas.

Gassmann’s theory provides the bulk modulus in fluid-saturated rock ($K_{\mathrm{Sat}}$) as a function of the dry rock bulk modulus ($K_{\mathrm{Dry}}$), the bulk modulus of the solid phase ($K_s$), that of the pore fluid ($K_f$), and total porosity ($\varphi$). It assumes that the shear modulus is fluid-independent

$$K_{\mathit{Sat}}=K_s\frac{\varphi K_{\mathit{Dry}}-\left(1+\varphi \right)K_f{K}_{\mathit{Dry}}/K_s+K_f}{\left(1-\varphi \right)K_f+\varphi K_s-K_f{K}_{\mathit{Dry}}/K_s},G_{\mathit{Sat}}=G_{\mathit{Dry}}.$$E17

The latter equation can be rearranged as follows:

$$K_{\mathit{Dry}}=K_s\frac{1-\left(1-\varphi \right)K_{\mathit{Sat}}/K_s-\varphi K_{\mathit{Sat}}/K_f}{1+\varphi -\varphi K_s/K_f-K_{\mathit{Sat}}/K_s},G_{\mathit{Dry}}=G_{\mathit{Sat}}.$$E18

Equations (17) and (18) provide us with a fluid substitution recipe as follows. Assume that we know the bulk modulus $K_{\text{SatA}}$ of rock saturated with Fluid A whose bulk modulus is $K_{\mathrm{fA}}$ and density is ${\rho }_{\mathrm{fA}}$. Then from Eq. (17), we obtain:

$$K_{\mathit{Dry}}=K_s\frac{1-\left(1-\varphi \right)K_{\mathit{\text{SatA}}}/K_s-\varphi K_{\mathit{\text{SatA}}}/K_{\mathit{fA}}}{1+\varphi -\varphi K_s/K_{\mathit{fA}}-K_{\mathit{\text{SatA}}}/K_s}.$$E19

The bulk modulus $K_{\text{SatB}}$ of the same rock saturated with Fluid B is (Eq. (17)):

$$K_{\text{SatB}}=K_s\frac{{\mathit{\varphi K}}_{\mathrm{Dry}}-\left(1+\varphi \right)K_{\mathrm{fB}}{K}_{\mathrm{Dry}}/K_s+K_{\mathrm{fB}}}{\left(1-\varphi \right)K_{\mathrm{fB}}+{\mathit{\varphi K}}_s-K_{\mathrm{fB}}{K}_{\mathrm{Dry}}/K_s},$$E20

where $K_{\mathrm{fB}}$ is the bulk modulus of Fluid B.

Of course, the shear modulus of the rock remains the same, no matter what fluid it is saturated with.

It is important to remember that the bulk density ${\rho }_b$ of the rock is also a function of the pore fluid. It depends on the porosity and density of the fluid (${\rho }_{\mathit{fA}}$ or ${\rho }_{\mathit{fB}}$):

$${\rho }_{\mathit{bB}}={\rho }_{\mathit{bA}}-{\varphi \rho }_{\mathit{fA}}+{\varphi \rho }_{\mathit{fB}},$$E21

where ${\rho }_{\mathit{bA}}$ and ${\rho }_{\mathit{bB}}$ are the bulk densities of the rock with the two pore fluids, respectively.

Finally, we can compute the elastic-wave velocities, as well as other seismic attributes, once we know the elastic moduli:

$$V_{\mathit{pB}}=\sqrt{\frac{K_{\mathit{\text{SatB}}}+\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.G_{\mathit{Dry}}}{{\rho }_{\mathit{bB}}}};V_{\mathit{sB}}=\sqrt{\frac{G_{\mathit{Dry}}}{{\rho }_{\mathit{bB}}}},$$E22

and

$$I_{\mathit{pB}}={\rho }_{\mathit{bB}}{V}_{\mathit{pB}};{\nu }_B=\frac{1}{2}\frac{{\left(V_{\mathit{pB}}/V_{\mathit{sB}}\right)}^2-2}{{\left(V_{\mathit{pB}}/V_{\mathit{sB}}\right)}^2-1},$$E23

where $I_{\mathrm{pB}}$ and ${\nu }_B$ are the P-wave impedance and Poisson’s ratio of the rock filled with Fluid B, respectively. Although the shear modulus $G$ is pore-fluid-independent, $V_s$ is since the bulk density varies with varying fluid.

Let us refer to a later important development in theoretical fluid substitution. It stemmed from the fact that Gassmann’s theory [9] requires the knowledge of the bulk modulus that can only be computed using Eq. (1) if both $V_p$ and $V_s$ (and the bulk density ${\rho }_b$) are known. In practice, the shear wave velocity may not be available. To address this issue, Mavko et al. [10] derived an approximate (but quite accurate) $V_p$—only fluid substitution theory that uses the compressional modulus $M={\rho }_b{V}_p^2$ instead of the bulk modulus $K$. The functional form in this theory is the same as that in Gassmann’s:

$$\begin{array}{l}{M}_{\mathit{Sat}}\approx M_s\frac{\varphi M_{\mathit{Dry}}-\left(1+\varphi \right)K_f{M}_{\mathit{Dry}}/M_s+K_f}{\left(1-\varphi \right)K_f+\varphi M_s-K_f{M}_{\mathit{Dry}}/M_s},\\ M_{\mathit{Dry}}\approx M_s\frac{1-\left(1-\varphi \right)M_{\mathit{Sat}}/M_s-\varphi M_{\mathit{Sat}}/K_f}{1+\varphi -\varphi M_s/K_f-M_{\mathit{Sat}}/M_s}.\end{array}$$E24

Figure 5 shows an example of the results of fluid substitution (pure water) on the elastic properties of high-porosity sand measured in the laboratory [11] at room-dry conditions. Clearly, the pore fluid has a dramatic effect on Poisson’s ratio. Such plots are basis for in situ fluid identification from seismic data.

Figure 5.

Sand experimental data and fluid substitution. Left. The bulk and shear moduli versus confining pressure as measured (dry) and water-substituted using Gassmann’s theory [9]. Middle. V_p and V_s versus confining pressure as measured (dry) and water-substituted. Right. The P-wave impedance versus Poisson’s ratio as measured (dry) and water-substituted, color-coded by the confining pressure.

Let us finally describe the details required in practical fluid substitution, specifically the computation of $K_s$, ${\rho }_s$, $K_f$, and ${\rho }_f$.

The elastic moduli of the multi-mineral rock matrix $K_s$ and $G_s$ can be obtained using Hill’s average (e.g., [1]) as

$$K_s=\frac{K_V+K_R}{2},G_s=\frac{G_V+G_R}{2},$$E25

where

$$\begin{array}{c}{K}_V=\sum _{i=1}^N{f}_i{K}_i,G_V=\sum _{i=1}^N{f}_i{G}_i,\\ K_R^{-1}=\sum _{i=1}^N{f}_i{K}_i^{-1},G_R^{-1}=\sum _{i=1}^N{f}_i{G}_i^{-1},\end{array}$$E26

where $N$ is the number of the mineral components, $f_i$ is the volume fraction of $i$^th mineral, and $K_i$ and $G_i$ are the bulk and shear moduli of the $i$^th component. The pure-mineral elastic moduli, as well as their densities, can be found in various sources, including Mavko et al. [1].

The bulk modulus of the pore fluid is

$$\frac{1}{K_f}=\frac{S_w}{K_w}+\frac{S_o}{K_o}+\frac{S_g}{K_g},$$E27

where $K_w$, $K_o$, and $K_g$ are the bulk moduli of water, oil, and gas, respectively. To estimate these moduli, as well as the densities used in Eq. (5), we refer to [12].

3. Variables influencing the elastic properties of rocks

In addition to the pore fluid, there are two more important variables influencing the elastic properties of rocks, their mineralogy and the differential pressure $P_{\text{Diff}}$ (or stress) defined as the difference between the confining $P_{\text{Confining}}$ (the overburden) and pore pressure $P_{\text{Pore}}$:

$$P_{\text{Diff}}=P_{\text{Confining}}-P_{\text{Pore}}.$$E28

Of course there are other influencing factors, such as rock texture (clastics versus carbonates versus unconventional shale), temperature, and diagenetic history. Here we only concentrate on the abovementioned two.

Mineralogy. As an example, let us examine the Han [13] laboratory dataset obtained on a large suite of sandstones with porosity ranging from zero to 30% and clay content between zero and 50%. Figure 6 shows $V_p$ and $V_s$ versus porosity and color-coded by the clay content.

Figure 6.

Dry rock V_p (top) and V_s (bottom) versus porosity, color-coded by the clay content, at confining pressure 50 MPa (left) and 5 MPa (right) (after Han [13]).

Obviously, the clay content plays a dramatic role acting to reduce both $V_p$ and $V_s$ at the same porosity. Also notice that the velocity-porosity-mineralogy trends are much more pronounced at 50 MPa. This is a commonly observed effect due to much clearer manifestations of key rock properties at high confining stress. The high-porosity data point in Figure 6 at porosity about 0.33 is for unconsolidated Ottawa sand sample. The effect of pressure on its velocities is very strong, similar to what we observe in Figure 5 for a sand of different provenance.

Another striking example of velocity discrimination due to mineralogy comes from unconventional shale with data obtained by wireline logging in a vertical well (Figure 7). The data shown in this figure is for 100% wet rock, obtained by fluid substitution from in situ conditions. The velocity-porosity dataset forms an amorphous cloud (Figure 7, top) with both $V_p$ and $V_s$ varying by almost 1.5 km/s at the same porosity. However, as soon as we introduce a third variable, the sum of the clay and kerogen contents, we observe a clear velocity discrimination with the velocity decreasing as the fraction of this softest component of the solid matrix increasing (Figure 7, bottom).

Figure 7.

100% wet rock V_p (left) and V_s (right) without accounting for mineralogy (top) and color-coded by the sum of clay and kerogen contents (bottom) (adopted from Dvorkin et al. [14]).

The Raymer et al. [3] model also predicts a strong dependence of the velocity on mineralogy (Figure 8), as well as the pore fluid, the latter well pronounced at higher porosity.

Figure 8.

V_p versus porosity according to the Raymer et al. [3] model for dry rock (left) and 100% water-saturated rock (right). The mineralogy is quartz and clay. The upper curves is for zero clay content, while the bottom curve is for 100% clay. The in-between curves are for gradually increasing clay content with increment 10% (top to bottom).

Stress. The effect of the confining pressure on the velocity in sand can be clearly seen in Figure 5 with $V_p$ in dry rock increasing by about 2.0 km/s and $V_s$ by 1.5 km/s as $P_{\text{Diff}}$ varies from almost zero to 50 MPa. Two more examples are shown in Figure 9, one for a sandstone from the North Sea and the other from the Gulf of Mexico.

Figure 9.

V_p (left) and V_s (right) versus pressure for two dry sandstone samples from the North Sea (top) and Gulf of Mexico (bottom).

The velocity in carbonate rocks is often not as affected by stress as it is in clastic samples. The magnitude of this effect is often influenced by the presence of compliant cracks in the rock. Such cracks act to strongly affect the velocity at low pressure while they are open. As the pressure increases, these cracks close acting to increase the velocity (Figure 10, Sample A). In samples where the cracks are absent, the velocity hardly varies as a function of pressure (Figure 10, Sample B).

Figure 10.

Same as Figure 9 but for two chalk samples from Ekofisk field in the North Sea. Velocities in sample A (porosity 0.38) is pressure dependent, while the velocities in sample B (porosity 0.31) hardly vary with pressure.

Notice that both historic velocity-porosity model by Wyllie et al. [2] and Raymer et al. [3] do not account for the dependence of the elastic-wave velocities on the confining stress. Both models are suitable for predicting the elastic properties at high, but not at low stress.

The velocity-stress dependence is important in understanding and predicting the seismic responses during hydrocarbon recovery, a process where the differential pressure may increase during production if the reservoir is depleted and the pore pressure is reduced, while the overburden remains constant. This differential pressure may decrease during enhanced oil recovery where water or gas are injected into the reservoir at high pressure, acting to reduce the difference between the overburden and pore pressure. Plots similar to that shown in Figure 5 (right-hand frame) are useful in simultaneously assessing the effects of the pore fluid and differential pressure on the elastic attributes.

4. Theoretical velocity-porosity models

There are two kinds of elastic moduli versus porosity effective medium models: (a) inclusion models and (b) grain-based models. The first kind models build a rock from the zero-porosity endpoint by placing inclusions into the solid matrix [1]. These models are perhaps relevant to some carbonate rocks where the pores appear as inclusions in calcite or dolomite matrix. The second kind assumes that the rock is formed by solid grains which comprise an uncemented grain pack at the high-porosity endpoint (also called the critical porosity) and, as the porosity is reduced, the original pack is altered either by grain contact cement or by smaller grains deposited in the pore space between the original larger grains, or a combination of these two processes.

As an example of the inclusion models, consider the differential effective medium model (DEM), where spheroidal pores are placed inside the solid matrix. A spheroid is an ellipsoid with two large diameters equal to each other and the third diameter smaller or equal to these two. The ratio of the small to large diameter is called the aspect ratio $\alpha \le 1$. If the spheroid is a sphere, $\alpha$ = 1. The inputs are the bulk and shear moduli of the mineral matrix and those of the inclusions.

Figure 11 (top) shows how the bulk and shear moduli depend on the total porosity for pure calcite rock with the bulk and shear moduli of the mineral 76.8 and 32.0 GPa, respectively, and its density 2.71 g/cc. The pores are empty, meaning the bulk and shear moduli of the inclusions are zero. In the same figure (bottom), we plot the respective $V_p$ and $V_s$. The aspect ratio is different for each of the curves shown. It is 1.00 for the upper curves and gradually decreases to 0.50, 0.20, 0.10, and 0.01 for the curves below. The smaller the aspect ratio, the smaller the elastic moduli and velocities at a fixed porosity.

Figure 11.

Elastic moduli (top) and velocities (bottom) versus porosity computed using DEM model for a pure calcite rock. The aspect ratio corresponding to the top curves is 1.00 and for the bottom curve 0.01. The aspect ratio gradually decreases to 0.50, 0.20, and 0.10 for the curves in between (top to bottom).

Figure 12 is the same as Figure 11 except that we use a single aspect ratio 0.10 and compare the results for empty inclusions with those for water-filled inclusions where the bulk modulus is 2.25 GPa and density is 1.00 g/cc.

Figure 12.

Same as Figure 11 but for a single aspect ratio 0.10 and for empty pores (black) and pores filled with water (blue).

We observe that both the bulk and shear moduli increase for pores filled with water as compared to empty pores. So do $V_p$ and $V_s$. This means that DEM is not consistent with Gassmann’s fluid substitution theory [9] which predicts that the shear modulus is pore-fluid-independent and $V_s$ reduces upon saturation due to increasing bulk density.

Notice that DEM curves connect two endpoints, one at zero porosity where the elastic moduli of rock are those of the mineral matrix and the other at 100% porosity where the elastic moduli are those of the inclusions (fluid in the pores). About three decades ago, Nur observed that most natural rocks simply do not exist in the entire zero to 100% porosity range. The maximum geologically plausible porosity for clastic rocks (sands and sandstones) is about 0.40. It may be higher in carbonates, such as chalks, that can have porosity up to 0.50. This porosity can be even higher for foam-like formations, such as volcanic rock (pumice) or artificially manufactured glass foam. This maximum porosity is called the critical porosity. This concept was formalized in Nur et al. [15].

One implication of the critical porosity concept is that the high-porosity endpoint should be at the critical porosity rather than at 100% porosity. It gave rise to the so-called modified elastic bounds. The simplest example is based on the upper elastic bound (also called the Voigt bound) for a composite made of two elastic components (“1” and “2”) with the compressional and shear moduli $M_1$, $G_1$ and $M_2$, $G_2$, respectively.

Assume that $M_2$ = $G_2$ = 0. Then the respective moduli of this composite ($M$ and $G$) at porosity $\varphi$ cannot exceed

$$M=\left(1-\varphi \right)M_1;G=\left(1-\varphi \right)G_1.$$E29

These two curves are plotted in Figure 13. In the same figure, we plot Han’s [13] data for low-clay-content samples at 50 MPa confining pressure. These data fall way below the upper bound curves for pure quartz with $M_1$ = 96.6 GPa and $G_1$ = 45.0 GPa.

Figure 13.

Upper and modified upper elastic bounds for the compressional (left) and shear (right) moduli versus porosity. The critical porosity is 0.36. Data are from Han’s [13] sandstone dataset for the clay content below 7% and with the elastic-wave velocities measured on dry samples at 50 MPa confining pressure.

The modified bounds use the same equations, but with porosity scaled by the critical porosity ${\varphi }_c$:

$$M=\left(1-\varphi /{\varphi }_c\right)M_1;G=\left(1-\varphi /{\varphi }_c\right)G_1;\varphi \le {\varphi }_c,$$E30

giving modified curves that are much closer to the data (Figure 13).

All grain-based theories exploit the critical porosity concept. We start with the contact-cement theory where it is assumed that the grains are not subjected to any confining stress at ${\varphi }_c$, and, as a result, the elastic moduli are zero, and porosity reduction is due to cement rims enveloping the grains (Figure 14). Such contact cement acts to rapidly increase the elastic moduli of the grain pack due to the dramatically expanding contact areas between the grains as porosity decreases, as explained in Dvorkin et al. [4], where the theoretical equations are given as well. This model is only valid in the very high-porosity range.

Figure 14.

Schematic modes of porosity reduction. From top to bottom: Contact-cement and stiff-sand model; soft-sand model; and constant-cement model (adopted from Dvorkin et al. [4]).

The soft-sand model assumes that at the critical porosity and the elastic properties of the grain pack are given by the Hertz-Mindlin [16] contact theory. This theory assumes that the grain pack is made of identical spherical grains whose elastic properties are those of the mineral (solid) matrix as given by Eq. (25). Combined with the mean field approximation that assumes that all grains are subject to identical local stresses and have the same average number of contacts per grain $n$ (also called the coordination number), the respective dry rock bulk ($K_{\mathrm{HM}}$) and shear ($G_{\mathrm{HM}}$) moduli are

$$K_{\mathit{HM}}={\left[\frac{n^2{\left(1-{\varphi }_c\right)}^2{G}_s^2}{18{\pi }^2{\left(1-{\nu }_s\right)}^2}P\right]}^{\frac{1}{3}},G_{\mathit{HM}}=\frac{5-4{\nu }_s}{5\left(2-{\nu }_s\right)}{\left[\frac{3n^2{\left(1-{\varphi }_c\right)}^2{G}_s^2}{2{\pi }^2{\left(1-{\nu }_s\right)}^2}P\right]}^{\frac{1}{3}},$$E31

where $P$ is the differential pressure (Eq. (28)) and $G_s$ and ${\nu }_s$ are the shear modulus and Poisson’s ratio of the solid matrix, respectively. This model implies that porosity reduction is not due to contact-cement deposition but instead due to smaller particles deposited away from grain contacts (Figure 14).

The coordination number $n$ in an identical grain pack at the critical porosity is about 6.

It is assumed in Eq. (31) that the grains have infinite friction (no slip) at their contacts. If we allow only the fraction $f$ of these contacts to have infinite friction while the rest of the contacts are frictionless and can slip, the equation for $K_{\mathrm{HM}}$ does not change but $G_{\mathrm{HM}}$ becomes now

$$G_{\mathit{HM}}=\frac{2+3f-\nu \left(1+3f\right)}{5\left(2-\nu \right)}{\left[\frac{3n^2{\left(1-{\varphi }_c\right)}^2{G}^2}{2{\pi }^2{\left(1-\nu \right)}^2}P\right]}^{\frac{1}{3}}.$$E32

This parameter $f$ is called the shear stiffness correction factor.

Finally, to obtain the dry rock bulk ($K_{\text{Soft}}$) and shear ($G_{\text{Soft}}$) moduli at any porosity $\varphi$ < ${\varphi }_c$, we use the modified (critical porosity scaled) lower Hashin-Shtrikman bound (e.g., [4]):

$$\begin{array}{l}{K}_{\mathit{\text{Soft}}}={\left[\frac{\varphi /{\varphi }_c}{K_{\mathit{HM}}+\frac{4}{3}{G}_{\mathit{HM}}}+\frac{1-\varphi /{\varphi }_c}{K+\frac{4}{3}{G}_{\mathit{HM}}}\right]}^{-1}-\frac{4}{3}{G}_{\mathit{HM}},\\ G_{\mathit{\text{Soft}}}={\left[\frac{\varphi /{\varphi }_c}{G_{\mathit{HM}}+z_{\mathit{HM}}}+\frac{1-\varphi /{\varphi }_c}{G+z_{\mathit{HM}}}\right]}^{-1}-z_{\mathit{HM}},z_{\mathit{HM}}=\frac{G_{\mathit{HM}}}{6}\left(\frac{9K_{\mathit{HM}}+8G_{\mathit{HM}}}{K_{\mathit{HM}}+2G_{\mathit{HM}}}\right).\end{array}$$E33

It is important to emphasize that the critical porosity endpoints here do not necessarily have to be given by the Hertz-Mindlin contact theory. Alternatively, these values can be selected from experimental data. What is most important in this model is the usage of the “soft” connection between the two porosity endpoints.

An alternative “stiff” connection between the aforementioned endpoints is given by the modified upper Hashin-Shtrikman bound as

$$\begin{array}{l}{K}_{\mathit{\text{Stiff}}}={\left[\frac{\varphi /{\varphi }_c}{K_{\mathit{HM}}+\frac{4}{3}{G}_s}+\frac{1-\varphi /{\varphi }_c}{K+\frac{4}{3}{G}_s}\right]}^{-1}-\frac{4}{3}{G}_s,\\ G_{\mathit{\text{Stiff}}}={\left[\frac{\varphi /{\varphi }_c}{G_{\mathit{HM}}+z}+\frac{1-\varphi /{\varphi }_c}{G_s+z}\right]}^{-1}-z,z=\frac{G_s}{6}\left(\frac{9K_s+8G_s}{K_s+2G_s}\right),\end{array}$$E34

where, once again, $G_s$ and $K_s$ are the shear and bulk moduli of the solid matrix, respectively.

This stiff connection, also called the stiff-sand model, can serve to connect the contact-cement curve with the zero-porosity endpoint.

Yet another model belonging to this family is the constant-cement model. It assumes that the grains have initial contact cementation with further porosity reduction due to the placement of small particles away from grain contacts (Figure 14). The functional form of this model is the same as in the soft-sand model (Eq. (33)) but with artificially high coordination number.

Examples of velocity-porosity curves according to the aforementioned grain-based theories are shown in Figure 15, where we assumed that both the grain and cement materials are pure quartz; $n$ for the soft-sand model is 6, while it is 20 for the constant-cement model; and the differential pressure is 20 MPa. The shear stiffness correction factor is 1.

Figure 15.

Velocity-porosity curves according to the soft-sand, stiff-sand, contact-cement, and constant-cement models as explained in the text.

Figure 16 shows an example of using the constant-cement model to describe the elastic behavior of unconventional gas shale, while Figure 17 is an example of applying the stiff-sand model to carbonate reservoirs. The parameters of the models are provided in the captions. These two examples show that the grain-based theories given here are appropriate not only for clastic sediments but also in very different lithological settings.

Figure 16.

V_p (left) and V_s (right) versus porosity for gas shale from wireline data adjusted for 100% water saturation. The color code is the sum of the clay and kerogen volume fractions (red for high and blue for low). The model curves are computed to bound the data. These curves are from the constant-cement model with the coordination number 12, differential pressure 26 MPa, critical porosity 0.40, and shear stiffness correction factor 1 (adopted from Dvorkin et al. [14]).

Figure 17.

Velocity- (left) and impedance-porosity (middle) plots showing chalk (gray) and lower-porosity carbonate (black) data points from wireline data adjusted for 100% water saturation. Graph on the left is the impedance versus Poison’s ratio plot, also for 100% water saturation conditions. The curves are from the stiff-sand model with the coordination number 6, differential pressure 30 MPa, critical porosity 0.40, and shear stiffness correction factor 1. The two model curves are for the two slightly different properties of the pure calcite end member (adopted from Dvorkin and Alabbad [17]).

Figure 18 shows laboratory data obtained at 30 MPa confining pressure on dry high-porosity, almost pure-quartz sand samples from the North Sea. In this classic example, the higher-velocity dataset is contact-cemented turbidite sand, while the lower-velocity dataset is friable and virtually uncemented sand. The former data can be matched by the contact-cement curves transitioning into the stiff-sand trajectories. The latter data are matched by the soft-sand curves.

Figure 18.

V_p (left) and V_s (right) versus porosity for two high-porosity sand datasets as explained in the text. The model curves marked in the plots are computed for 30 MPa differential pressure, critical porosity 0.40, coordination number 7, shear stiffness correction factor 1, and dry rock.

5. Digital rock physics

Digital rock physics is based on the concept “image and compute,” image rock at the pore scale (Figure 19) and then simulate in the computer various processes in such an image to arrive at a desired rock property. These simulations include viscous fluid flow to arrive at hydraulic permeability, electrical charge flow to arrive at electrical resistivity, as well as elastic deformation to arrive at the elastic moduli and velocities.

Figure 19.

Segmented digital images of loose sand (porosity about 30%), sandstone (porosity about 20%), and carbonate (porosity about 15%) showing the mineral matrix and pores. The images are a few mm across.

The advantage of such digital approach is that the same sample can be reused multiple times, unlike in physical experiments where a sample is altered after every test; the sample can be digitally altered by, e.g., introducing digenetic cementation, which is hardly possible in physical experiments, as well as subsampling of a digital volume to investigate how various rock properties vary within the volume and how relations between rock properties depend on the spatial scale of investigation.

Although the aforementioned concept is simple, its implementation is not. First, the imaging has to be conducted at the appropriate scale and resolution to reveal the salient features of natural rock relevant to the process under examination. Second, the image has to be segmented to separate minerals from pores and segregate various minerals within the solid matrix, as well as fluid phases inside the pores. Third, powerful computational engines have to be utilized and verified to simulate processes relevant to the physical experiment.

In spite of these complexities, during the last decade, DRP has emerged as a powerful technique complementing (if not replacing) physical testing, mostly due to the recent advances in imaging hardware and image processing and computational software, the latter combined with steadily improving computational power. Not only DRP has become a novel research tool in academia and national labs, but is has also been adopted by leading oil and service companies.

There is one more inherent feature of DRP that needs to be accounted for. Pore-scale rock images are only a few mm in size, and the higher the resolution needed to revel the salient features, the smaller the field of view. At the same time, these computational results have to be relevant at much larger spatial scales of feet for wireline measurement interpretations in the well or tens and hundreds of feet in seismic prospecting. Even such basic property as porosity may be different if measured on an inch-sized sample an on mm-sized fragment of the same sample.

One way out of this conundrum is instead of directly comparing data points generated by different methods of measurement, compare trends formed by such data points, such as permeability versus porosity trends. Dvorkin et al. [18] show that such trends are often hidden inside a very small digital sample and can be derived by subsampling it. Moreover, these computational trends often match relevant physical trends and/or theoretical rock physics transforms, hence validating computational results and making them relevant at much coarser spatial scales.

The approach is to subsample a digital volume into 2³, 3³, or 4³ subvolumes (Figure 20) and then compute the desired property pairs (e.g., porosity and permeability) on each of these subvolumes. Very often, the property pairs thus computed form a meaningful trend supported by physical measurements and/or theories (see examples in Figures 21–23). We can call this subsampling approach “to see the rock in a grain of sand.”

Figure 20.

Illustration of the subsampling approach.

Figure 21.

Permeability versus porosity in Fontainebleau sandstone. Left: Laboratory data matched with a Kozeny-Carman theoretical curve. Right: Multiple permeability versus porosity data points computed from a few digital Fontainebleau samples and subsamples thereof (adopted from Dvorkin et al. [18]).

Figure 22.

Formation factor versus porosity computed on carbonate cuttings. The curves are from Archie’s equation with the cementation exponent m 2.0, 2.5, and 3.0 (bottom to top) (adopted from Dvorkin et al. [18]).

Figure 23.

V_p versus porosity for Fontainebleau sandstone as computed from a few digital samples and subsamples thereof (squares). Gray circles are from laboratory measurements of dry samples. The curve is from the stiff-sand model (adopted from Dvorkin et al. [18]).

These results open ways to a meaningful utilization of DRP in research and industry. Publications related to DRP are many and the number is growing. We refer the reader to Kameda and Dvorkin [19], Dvorkin et al. [20], Dvorkin and Derzhi [21], and Andra et al. [22, 23].

6. Conclusion

This chapter presents an overview of rock physics, starting with its history and ending with the most recent development, the digital rock physics. This chapter can be used as a basic reference pointing towards published sources where the topic is developed in-depth and detailed equations, tables, and experimental results are given. One of such comprehensive sources is the third edition of the Rock Physics Handbook [24].

Rock physics remains a key component in interpreting seismic and other remote sensing data for the underlying properties and conditions of the subsurface. A plethora of such practical results has appeared and continues to appear in geophysical journals, such as Geophysics (Society of Exploration Geophysicists), Journal of Geophysical Research (American Geophysical Union), and First Break (European Association of Geoscientists and Engineers), as well as presented at conferences worldwide.

An important topic not addressed in this chapter is a simultaneous interpretation of different remote sensing sources, such as seismic prospecting, electric and electromagnetic sensing, and gravity methods. Once again, such materials can be found in the proceedings and books from the aforementioned professional societies.

We feel that the material presented can serve as a detailed introduction into the extensive field of physics of rocks and be of use to graduate students, as well as advanced professional in earth and environmental sciences.