Rock Physics: Recent History and Advances

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.


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: 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 ρ b is the bulk density. The latter quantity is related to the total porosity ϕ as where ρ s is the density of the mineral matrix also called the solid component of the rock, while ρ 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 (ν): Most of natural rocks contain more than one mineral. In this situation, ρ s can be computed as the arithmetic average of the densities of the individual components: where f i is the volume fraction of the i-th mineral component in the mineral matrix and ρ 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: where S w , S o , and S g are the water, oil, and gas saturations in the pore space, respectively, and ρ w , ρ o , and ρ g are the densities of these pore fluid components. Of course, it is required that X N i¼1 f i ¼ 1 (6) and S w þ S o þ S g ¼ 1: 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 "velocityporosity-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 ϕ 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 where V p is the P-wave velocity, V ps is the velocity in the mineral phase, and V 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 pf = 1500 m/s. 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 pf = 0. Yet, as have been shown by seismic experiments on the moon, V p in such sediment is finite.  [2] and Raymer et al. [3] transforms for quartz, dolomite, and mixed mineralogy. Legend in the middle refers to all plots. 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 The next historic equation was introduced by Raymer et al. [3]: 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 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 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.
At partial brine saturation, S w < 1, the resistivity of rock R tS not only depends on porosity but also on saturation S w as 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 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 α. The tortuosity τ of these pores is defined as The permeability k is also a function of the specific surface area S defined as the ratio of surface of the pore space S Pore to the total volume V of the rock sample: A variable alternative to S is the grain size (or grain diameter) d.
The Kozeny-Carman equation reads [1] k ¼ 1 2 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 ϕ p called the percolation porosity: 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 À13 m 2 , while one mD = 10 À15 m 2 . 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, τ = 2.5, and ϕ 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].

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]).
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 Sat )asa function of the dry rock bulk modulus (K Dry ), the bulk modulus of the solid phase (K s ), that of the pore fluid (K f ), and total porosity (ϕ). It assumes that the shear modulus is fluid-independent The latter equation can be rearranged as follows: Equations (17) and (18) provide us with a fluid substitution recipe as follows. Assume that we know the bulk modulus K SatA of rock saturated with Fluid A whose bulk modulus is K fA and density is ρ fA . Then from Eq. (17), we obtain: The bulk modulus K SatB of the same rock saturated with Fluid B is (Eq. (17)): where K 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 ρ b of the rock is also a function of the pore fluid. It depends on the porosity and density of the fluid (ρ fA or ρ fB ): where ρ bA and ρ 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: and where I pB and ν 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-fluidindependent, 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 ρ 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 p instead of the bulk modulus K. The functional form in this theory is the same as that in Gassmann's: 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 roomdry 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.
Let us finally describe the details required in practical fluid substitution, specifically the computation of K s , ρ s , K f , and ρ 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 where 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 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.
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].

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 Diff (or stress) defined as the difference between the confining P Confining (the overburden) and pore pressure P Pore : 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. 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 highporosity 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).
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.
Stress. The effect of the confining pressure on the velocity in sand can be clearly seen in Figure 5    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  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.

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 highporosity 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 α ≤ 1. If the spheroid is a sphere, α = 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 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.
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 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 ϕ 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.
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  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 grainbased 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 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 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.

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.
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. 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]). 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. Porescale 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 ,3 3 ,or4 3 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." 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].

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.