# Rock physics

Rock physics aims to characterize rock properties based on the behavior of seismic waves propagating through them. This requires consideration of how the composition of a rock dictates its stress-strain relationship and thus seismic response. The effect of pore fluids is of particular interest due to its applicability to the hydrocarbon industry. In a standard seismic interpretation workflow rock physics is used to relate impedance and elastic parameters derived from seismic data to specific rock properties. This constrains what seismic data is physically capable of resolving and the non-uniqueness associated with a specific interpretation.

## Mechanics

### Hooke's Law

A material is linearly elastic if its stress-strain relationship is governed by Hooke's Law,

${\displaystyle \sigma _{ij}=c_{ijkl}\varepsilon _{kl}}$

Where ${\displaystyle \sigma _{ij}}$ is the second order stress tensor, ${\displaystyle \varepsilon _{kl}}$ is the second order strain tensor, and ${\displaystyle c_{ijkl}}$ is the fourth order elastic tensor, defining a material's resistance to strain due to a stress in each direction.

The fourth order elasticity or stiffness tensor requires 34 terms to specify the stress strain relationship for an arbitrary homogeneous material. Conservation of angular momentum requires the stress tensor to be symmetric, as any asymmetries would imply rotations, which are not described by stress. Stress describes deformation only. This also requires the strain tensor to be symmetric, implying that the elasticity tensor is symmetric, reducing the number of required components to 36. Energy conservation requires thermal energy and strain energy to balance, further reducing the number of required components to 21.

If one makes the simplifying assumption that the material is isotropic, that is, the strain experienced by a material is independent of the direction in which the stress is applied, the stiffness tensor requires only two parameters, Lamé's constants, λ and μ. Hooke's Law can then be expressed as

${\displaystyle \sigma _{ij}=\lambda \delta _{ij}\varepsilon _{kk}+2\mu \varepsilon _{ij}}$

### The Wave Equation

Equations describing the waves resulting from deformation of elastic materials can be derived by applying Hooke's Law to Cauchy's Equation of motion, given by:

${\displaystyle \rho {\frac {\partial ^{2}u_{i}}{\partial t^{2}}}={\frac {\partial \sigma _{ji}}{\partial x_{j}}}}$

Where ρ is the density of the material and u_{i} is the displacement of the material, related to the infinitesimal strain tensor by

${\displaystyle \varepsilon _{ij}={\tfrac {1}{2}}({\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}})}$

Hooke's Law is used to describe the stress tensor in Cauchy's Equation along with the above definition of the infinitesimal strain tensor to give an equation for material displacement as a function of time and space. Applying the divergence to this equation describes the rate at which volumetric perturbations propagate through the material giving the P-wave equation. For isotropic materials, this simplifies to

${\displaystyle {\frac {\partial ^{2}\Theta }{\partial t^{2}}}={\frac {\lambda +2\mu }{\rho }}\bigtriangledown ^{2}\Theta }$

${\displaystyle V_{p}={\sqrt {\frac {\lambda +2\mu }{\rho }}}}$

Where Θ is the dilatation, ${\displaystyle \Theta =\bigtriangledown \bullet {\vec {u}}}$

Similarly, applying the curl operator describes the rate at which shape perturbations propagate through a material giving the S-wave equation. For isotropic materials this gives

${\displaystyle {\frac {\partial ^{2}\omega }{\partial t^{2}}}={\frac {\mu }{\rho }}\bigtriangledown ^{2}\omega }$

${\displaystyle V_{s}={\sqrt {\frac {\mu }{\rho }}}}$

Where ω is the rotation vector, ${\displaystyle \omega =\bigtriangledown \times {\vec {u}}}$

It is important to note that these velocities assume a homogenous, isotropic, perfectly elastic material. All energy associated with deformation is conserved. No energy is lost to friction etc., so velocities do not depend on wave frequency.

## Elastic Moduli

The following definitions assume a linearly elastic, homogenous isotropic material.

### Young's Modulus

Young's modulus, E (Pa), describes a material's ability to stretch or compress by quantifying how resistant a material is to extensional deformation.

${\displaystyle \sigma _{ii}=E\varepsilon _{ii}}$

In the hydrocarbon industry Young's modulus is used to design hydraulic fracturing mechanisms. Rigid rocks with a high value of Young's modulus will result in narrower fractures than rocks with a lower rigidity.

### Shear Modulus

The shear modulus, µ (Pa) or G describes the rigidity of a material, its resistance to shear strain/angular distortions.

${\displaystyle \sigma _{ij}=2\mu \varepsilon _{ij}}$ for i≠j

Note that fluids have a shear modulus of zero as they cannot withstand any shear stress.

### Bulk Modulus

The bulk modulus, K (Pa) of a material describes its ability to withstand volumetric stress,

${\displaystyle \sigma _{ii}=3K\Theta }$

The inverse of the bulk modulus is defined as the compressibility, β

${\displaystyle K={\frac {1}{\beta }}}$

The bulk modulus is related to the Lamé parameters by

${\displaystyle K=\lambda +{\tfrac {2}{3}}\mu }$

When compared to the isotropic form of Hooke's Law, this expresses that volumetric deformation depends on normal and shear stiffnesses. The bulk modulus quantifies its total volumetric stiffness. The Lamé parameter λ can therefore be thought of as a material's resistance to exclusively normal stresses, whereas µ quantifies the material's resistance to shape deformations, that can result from both normal and non-normal stress regimes.

### Poisson's Ratio

Poisson's ratio, ν defines the transverse strain that results from longitudinal strain. That is, when a material is stretched in one direction, Poisson's ratio defines how much it will contract in the other directions. For an isotropic material deformed in the x direction,

${\displaystyle \nu =-{\frac {d\varepsilon _{y}}{d\varepsilon _{x}}}=-{\frac {d\varepsilon _{z}}{d\varepsilon _{x}}}}$

The value of Poisson's Ratio varies from -1 to 0.5. A material is perfectly incompressible if its Poisson's ratio is 0.5, such as rubber. Materials that do not deform transversally due to longitudinal strain have a Poisson ratio of 0, such as cork. Materials with negative Poisson's ratios expand in transverse directions when stretched longitudinally due to their hinge-like molecular structures and are known as auxetic materials. A material is described as an ideal Poisson solid when λ=µ giving ν=0.25. This is roughly true for many rocks.

Poisson's ratio can also be expressed in terms of the materials P and S wave velocities,

${\displaystyle \nu ={\frac {({\tfrac {V_{p}}{V_{s}}})^{2}-2}{2({\tfrac {V_{p}}{V_{s}}})^{2}-2}}}$

Because µ=0 for fluids, Vs=0, making Poisson's ratio very sensitive to changes in fluid content.

### Moduli Relations

Conversion formulas
Homogeneous isotropic linear elastic materials have their elastic properties uniquely determined by any two moduli among these; thus, given any two, any other of the elastic moduli can be calculated according to these formulas.
${\displaystyle K=\,}$ ${\displaystyle E=\,}$ ${\displaystyle \lambda =\,}$ ${\displaystyle \mu =\,}$ ${\displaystyle \nu =\,}$ ${\displaystyle M=\,}$ Notes
${\displaystyle (K,\,E)}$ ${\displaystyle K}$ ${\displaystyle E}$ ${\displaystyle {\tfrac {3K(3K-E)}{9K-E}}}$ ${\displaystyle {\tfrac {3KE}{9K-E}}}$ ${\displaystyle {\tfrac {3K-E}{6K}}}$ ${\displaystyle {\tfrac {3K(3K+E)}{9K-E}}}$
${\displaystyle (K,\,\lambda )}$ ${\displaystyle K}$ ${\displaystyle {\tfrac {9K(K-\lambda )}{3K-\lambda }}}$ ${\displaystyle \lambda }$ ${\displaystyle {\tfrac {3(K-\lambda )}{2}}}$ ${\displaystyle {\tfrac {\lambda }{3K-\lambda }}}$ ${\displaystyle 3K-2\lambda \,}$
${\displaystyle (K,\,\mu )}$ ${\displaystyle K}$ ${\displaystyle {\tfrac {9K\mu }{3K+\mu }}}$ ${\displaystyle K-{\tfrac {2\mu }{3}}}$ ${\displaystyle \mu }$ ${\displaystyle {\tfrac {3K-2\mu }{2(3K+\mu )}}}$ ${\displaystyle K+{\tfrac {4\mu }{3}}}$
${\displaystyle (K,\,\nu )}$ ${\displaystyle K}$ ${\displaystyle 3K(1-2\nu )\,}$ ${\displaystyle {\tfrac {3K\nu }{1+\nu }}}$ ${\displaystyle {\tfrac {3K(1-2\nu )}{2(1+\nu )}}}$ ${\displaystyle \nu }$ ${\displaystyle {\tfrac {3K(1-\nu )}{1+\nu }}}$
${\displaystyle (K,\,M)}$ ${\displaystyle K}$ ${\displaystyle {\tfrac {9K(M-K)}{3K+M}}}$ ${\displaystyle {\tfrac {3K-M}{2}}}$ ${\displaystyle {\tfrac {3(M-K)}{4}}}$ ${\displaystyle {\tfrac {3K-M}{3K+M}}}$ ${\displaystyle M}$
${\displaystyle (E,\,\lambda )}$ ${\displaystyle {\tfrac {E+3\lambda +R}{6}}}$ ${\displaystyle E}$ ${\displaystyle \lambda }$ ${\displaystyle {\tfrac {E-3\lambda +R}{4}}}$ ${\displaystyle {\tfrac {2\lambda }{E+\lambda +R}}}$ ${\displaystyle {\tfrac {E-\lambda +R}{2}}}$ ${\displaystyle R={\sqrt {E^{2}+9\lambda ^{2}+2E\lambda }}}$
${\displaystyle (E,\,\mu )}$ ${\displaystyle {\tfrac {E\mu }{3(3\mu -E)}}}$ ${\displaystyle E}$ ${\displaystyle {\tfrac {\mu (E-2\mu )}{3\mu -E}}}$ ${\displaystyle \mu }$ ${\displaystyle {\tfrac {E}{2\mu }}-1}$ ${\displaystyle {\tfrac {\mu (4\mu -E)}{3\mu -E}}}$
${\displaystyle (E,\,\nu )}$ ${\displaystyle {\tfrac {E}{3(1-2\nu )}}}$ ${\displaystyle E}$ ${\displaystyle {\tfrac {E\nu }{(1+\nu )(1-2\nu )}}}$ ${\displaystyle {\tfrac {E}{2(1+\nu )}}}$ ${\displaystyle \nu }$ ${\displaystyle {\tfrac {E(1-\nu )}{(1+\nu )(1-2\nu )}}}$
${\displaystyle (E,\,M)}$ ${\displaystyle {\tfrac {3M-E+S}{6}}}$ ${\displaystyle E}$ ${\displaystyle {\tfrac {M-E+S}{4}}}$ ${\displaystyle {\tfrac {3M+E-S}{8}}}$ ${\displaystyle {\tfrac {E-M+S}{4M}}}$ ${\displaystyle M}$

${\displaystyle S=\pm {\sqrt {E^{2}+9M^{2}-10EM}}}$

There are two valid solutions.
The plus sign leads to ${\displaystyle \nu \geq 0}$.
The minus sign leads to ${\displaystyle \nu \leq 0}$.

${\displaystyle (\lambda ,\,\mu )}$ ${\displaystyle \lambda +{\tfrac {2\mu }{3}}}$ ${\displaystyle {\tfrac {\mu (3\lambda +2\mu )}{\lambda +\mu }}}$ ${\displaystyle \lambda }$ ${\displaystyle \mu }$ ${\displaystyle {\tfrac {\lambda }{2(\lambda +\mu )}}}$ ${\displaystyle \lambda +2\mu \,}$
${\displaystyle (\lambda ,\,\nu )}$ ${\displaystyle {\tfrac {\lambda (1+\nu )}{3\nu }}}$ ${\displaystyle {\tfrac {\lambda (1+\nu )(1-2\nu )}{\nu }}}$ ${\displaystyle \lambda }$ ${\displaystyle {\tfrac {\lambda (1-2\nu )}{2\nu }}}$ ${\displaystyle \nu }$ ${\displaystyle {\tfrac {\lambda (1-\nu )}{\nu }}}$ Cannot be used when ${\displaystyle \nu =0\Leftrightarrow \lambda =0}$
${\displaystyle (\lambda ,\,M)}$ ${\displaystyle {\tfrac {M+2\lambda }{3}}}$ ${\displaystyle {\tfrac {(M-\lambda )(M+2\lambda )}{M+\lambda }}}$ ${\displaystyle \lambda }$ ${\displaystyle {\tfrac {M-\lambda }{2}}}$ ${\displaystyle {\tfrac {\lambda }{M+\lambda }}}$ ${\displaystyle M}$
${\displaystyle (\mu ,\,\nu )}$ ${\displaystyle {\tfrac {2\mu (1+\nu )}{3(1-2\nu )}}}$ ${\displaystyle 2\mu (1+\nu )\,}$ ${\displaystyle {\tfrac {2\mu \nu }{1-2\nu }}}$ ${\displaystyle \mu }$ ${\displaystyle \nu }$ ${\displaystyle {\tfrac {2\mu (1-\nu )}{1-2\nu }}}$
${\displaystyle (\mu ,\,M)}$ ${\displaystyle M-{\tfrac {4\mu }{3}}}$ ${\displaystyle {\tfrac {\mu (3M-4\mu )}{M-\mu }}}$ ${\displaystyle M-2\mu \,}$ ${\displaystyle \mu }$ ${\displaystyle {\tfrac {M-2\mu }{2M-2\mu }}}$ ${\displaystyle M}$
${\displaystyle (\nu ,\,M)}$ ${\displaystyle {\tfrac {M(1+\nu )}{3(1-\nu )}}}$ ${\displaystyle {\tfrac {M(1+\nu )(1-2\nu )}{1-\nu }}}$ ${\displaystyle {\tfrac {M\nu }{1-\nu }}}$ ${\displaystyle {\tfrac {M(1-2\nu )}{2(1-\nu )}}}$ ${\displaystyle \nu }$ ${\displaystyle M}$

The stiffness matrix (9 by 9, or 6 by 6 in Voigt notation) in Hooke's law (in 3D) can be parametrized by only two components for homogeneous and isotropic materials. One may choose whichever pair one prefers among the elastic moduli given below. Some of the possible conversions are listed in the table.

## Effective Moduli

Most geophysical theories make the simplifying assumption that inherently heterogenous Earth materials obtain some effective, averaged homogenous properties over the dimension and scale of measurement. It is therefore essential to have an understanding of how the components of a material are arranged to define its overall physical properties. The properties of the whole can be significantly different than the properties of its parts.

The process of defining the effective properties of a heterogenous material is referred to as homogenization. Given a material with small scale, rapidly varying properties, homogenization attempts to redefine the material over a larger volume such that it becomes homogenous in a statistical sense at a particular scale. These homogenized volumes are referred to as representative volume elements (RVE’s). The parameters that describe the behavior of these RVE’s are defined as its effective parameters. A material can be considered homogenous if its constituents are smaller than the wavelength of measurement.

Seismic studies therefore require effective elastic moduli to be defined. This depends on the elastic properties of the constituents and their geometric relationships. Because it is essentially impossible to obtain sufficient information to uniquely define the effective moduli of a natural substance, the theoretical maximum and minimum moduli, known as the Voigt and Reuss bounds are often used instead.[1].

More refined effective moduli can be defined by accounting for the arrangement of the material's components. In the case of a rock, the mineral matrix components are similar enough that a simple Voigt-Reuss-Hill average is sufficiently accurate. Pores, pore filling fluids, clays and cements must be considered separately as they can significantly change the effective moduli. There are two general methods of doing so, depending on the relative volume fraction of the component being considered. For large relative volumes a granular approach is usually taken, where the effective moduli is determined by considering how the contact forces on grains change the overall moduli. Such methods are used for high porosity, more unconsolidated to conventional sands. For low relative volumes, the components are considered as isolated deviations from a background medium, suitable for low porosity, unconventional isotropic rocks.

### Voigt Bounds

The Voigt Bound is obtained by assuming the strain field remains constant throughout the material when subjected to an arbitrary average stress field. For a composite composed of N components, the effective Voigt moduli, Mv, are given by

${\displaystyle M_{v}=\sum _{i}^{N}f_{i}M_{i}}$

Where fi is the volume fraction of the ith component with a modulus Mi. It was proven by Hill [2] using energy conservation principles that this gives the maximum possible moduli of the composite.

### Reuss Bounds

The Reuss Bound is obtained by assuming the stress field remains constant throughout the material in an arbitrary average strain field. Energy conservation principles imply that gives the minimum possible moduli for the material. For a composite composed of N components, the effective Reuss moduli, MR, are given by

${\displaystyle {\frac {1}{M_{R}}}=\sum _{i}^{N}f_{i}{\frac {1}{M_{i}}}}$

Where fi is the volume fraction of the ith component with a modulus Mi.

### Voigt-Reuss-Hill average

A common simplistic approach to estimating the effective moduli is to calculate the average of the Voigt and Reuss bounds, referred to as the Voigt-Reuss-Hill (VRH) average.

${\displaystyle M_{VRH}={\frac {M_{V}+M_{R}}{2}}}$

This is relatively accurate when the moduli of the components have similar values and is commonly used to estimate the matrix moduli of multi-mineralic rocks.

### Hashin-Shtrikman Bounds

Hashin and Shtrikman[3] use a more robust approach using the calculus of variations. They derive the effective elastic moduli that minimize the potential elastic energy of the composite. The strain energy of the homogenous, isotropic reference material, U0 is defined then altered by introducing an inclusion with different elastic properties, and the resultant change in the elastic strain energy, δU=U-U0 is considered. An arbitrary number of inclusions can be added as long as the assumptions required for homogenization remain valid. The change in elastic strain energy due to the inclusion give maximum bounds when the reference material is the phase with the highest moduli and minimum bounds when the reference material is the phase with the lowest moduli. This is analogous to the method used to derive Voigt and Reuss bounds, except that instead of averaging stress and strain fields, Hashin minimizes the strain energy. Obtaining expressions for the effective moduli requires and assumption of spherical symmetry. Thus Hashin-Shtrikman bounds assume spherical inclusions distributed symmetrically about the matrix. For a two phase medium the resultant effective moduli are given by

${\displaystyle K_{HS}=K_{1}+{\frac {f_{2}}{(K_{2}-K_{1})^{-1}+f_{1}(K_{1}+{\tfrac {4}{3}}\mu _{1})^{-1}}}}$
${\displaystyle \mu _{HS}=\mu _{1}+{\frac {f_{2}}{(\mu _{2}-\mu _{1})^{-1}+{\tfrac {2f_{1}(K_{1}+2\mu _{1})}{5\mu _{1}(K_{1}+{\tfrac {4}{3}}\mu _{1})}}}}}$

Where the subscript 1 refers to the matrix material, and the subscript 2 refers to the inclusion material. The upper bound is obtained when the matrix has the maximum moduli, and the lower bound is obtained when the matrix has the minimum moduli. Note that if one of the components is a fluid, the lower bound for the shear modulus is zero, as the lower bound represents a suspension, which is effectively a fluid with a shear modulus of zero.

### Modified Hashin-Shtrikman Bounds

The Modified Hashin-Shtrikman Bounds separates consolidated and unconsolidated elastic moduli trends according to the rock's critical porosity. Consider the rock's diagenetic history beginning with an unconsolidated, highly porous sediment that gradually gets compacted until a critical porosity, ϕc (0.36-0.4 for sandstones) is reached. After this critical porosity is reached, cementation takes over as the primary diagenetic process that decreases porosity causing a faster increase in elastic moduli.

Empirically, it has been found that the cementation trend can be modelled using the Upper Hashin-Shtrikman bound by replacing the porosity, or fractional fluid content, ϕ with ${\displaystyle {\tfrac {\phi }{\phi _{c}}}}$, this is known as the Modified Upper Hashin-Shtrikman bound.

For porosities greater than the critical porosity, the compaction trend can be modelled using the Reuss bound, or the Modified Lower Hashin-Shtriman bound obtained by replacing the porosity ϕ with ${\displaystyle {\tfrac {\phi }{\phi _{c}}}}$ in the Lower Hashin-Shtrikman bound.

### Kuster and Toksöz formulation

Kuster and Toksöz [4]derive expressions for effective moduli of a two phase isotropic material using first order scattering theory. This formulation emphasizes the effect of inclusion shape by the incorporation of an aspect ratio term,α.

It is assumed that the concentration, f, of inclusions is small, so that interactions between inclusions can be neglected. The concentration of inclusions can be taken to be small if f/α<1, otherwise the inclusions overlap. This allows multiple scattering to be neglected in that the wave field incident on each inclusion is the original, undisturbed wave field, unaffected by other scatterers. It is also assumed that inclusions are randomly oriented, making the effective material isotropic.

To obtain equations for the effective moduli, an expression for the displacement, u(x) of an arbitrary material within an infinite matrix due to a wave, u0 incident from infinity is obtained using an integral equation of Green's function. The effective properties of a volume of material with N inclusions is taken to be the properties that equate the wave displacement of this entire volume, u*(x) to the sum of the displacement of each inclusion, ui(x), that is.

${\displaystyle u(x)=u_{0}(x)+\sum _{i}^{N}u_{i}(x_{i})=u_{0}(x)+u^{*}(x)}$

For a solid matrix with properties K μ and ρ, embedded with a fraction f of inclusions with properties K' μ' and ρ' and an aspect ratio α, the effective properties K* μ* and ρ* are then found to be given by:

${\displaystyle {\frac {K^{*}-K}{3K^{*}+4\mu }}=f{\frac {K'-K}{3K+4\mu }}{\frac {1}{3}}T_{iijj}}$

${\displaystyle {\frac {(\mu ^{*}-\mu )}{6\mu ^{*}(K+2\mu )+\mu (9K+8\mu )}}={\frac {f(\mu '-\mu )}{25\mu (3K+4\mu )}}(T_{ijij}-{\tfrac {1}{3}}T_{iijj})}$

${\displaystyle \rho ^{*}=\rho (1-f)+\rho 'f}$

Tijkl is Wu's inclusion tensor[5] which for ellipsoidal inclusions is given by

${\displaystyle T_{iijj}=-{\frac {3F_{1}}{F_{2}}}}$

${\displaystyle T_{ijij}-{\tfrac {1}{3}}T_{iijj}={\frac {2}{F_{3}}}+{\frac {1}{F_{4}}}+{\frac {F_{4}F_{5}+F_{6}F_{7}-F_{8}F_{9}}{F_{2}F_{4}}}}$

${\displaystyle F_{1}=1+A[{\tfrac {3}{2}}(g+h)-R({\tfrac {3}{2}}g+{\tfrac {5}{2}}h-{\tfrac {4}{3}})]}$

${\displaystyle F_{2}=1+A[1+{\tfrac {3}{2}}(g+h)-{\frac {R}{2}}(3g+5h)]+B(3-4R)+{\frac {A}{2}}(A+3B)(3-4R)[g+h-R(g-h+2h^{2})]}$

${\displaystyle F_{3}=1+{\frac {A}{2}}[R(2-h)+{\frac {(1+\alpha ^{2})}{\alpha ^{2}}}h(R-1)]}$

${\displaystyle F_{4}=1+{\frac {A}{4}}[3h+g-R(g-h)]}$

${\displaystyle F_{5}=A[R(g+h-{\tfrac {4}{3}})-g]+Bh(3-4R)}$

${\displaystyle F_{6}=1+A[1+g-R(g+h)]+B(1-h)(3-4R)}$

${\displaystyle F_{7}=2+{\frac {A}{4}}[9h+3g-R(5h+3g)]+Bh(3-4R)}$

${\displaystyle F_{8}=A[1-2R+{\tfrac {g}{2}}(R-1)+{\tfrac {h}{2}}(5R-3)]+B(1-h)(3-4R)}$

${\displaystyle F_{9}=A[g(R-1)-Rh]+Bh(3-4R)}$

${\displaystyle A={\frac {\mu '}{\mu }}-1,B={\tfrac {1}{3}}({\frac {K'}{K}}-{\frac {\mu '}{\mu }}),R={\frac {3\mu }{3K+4\mu }}}$

${\displaystyle h={\frac {\alpha }{(1-\alpha ^{2})^{3/2}}}[cos^{-1}\alpha -\alpha (1-\alpha ^{2})^{1/2}]}$

${\displaystyle g={\frac {\alpha ^{2}}{1-\alpha {2}}}(3h-2)}$

For materials with M different aspect ratios,

${\displaystyle {\frac {K^{*}-K}{3K^{*}+4\mu }}={\frac {K'-K}{3K+4\mu }}\sum _{i}^{M}f(\alpha _{i}){\tfrac {1}{3}}T_{iijj}(\alpha _{i})}$

${\displaystyle {\frac {(\mu ^{*}-\mu )}{6\mu ^{*}(K+2\mu )+\mu (9K+8\mu )}}={\frac {(\mu '-\mu )}{25\mu (3K+4\mu )}}\sum _{i}^{M}f(\alpha _{i})[T_{ijij}(\alpha _{i})-{\tfrac {1}{3}}T_{iijj}(\alpha _{i})]}$

More specific pore shapes the two phase Kuster Toksöz formulation reduces to [6]

${\displaystyle (K^{*}-K){\frac {K+{\tfrac {4}{3}}\mu }{K^{*}+{\tfrac {4}{3}}\mu }}=\sum _{i}^{M}f_{i}(K_{i}-K)P_{i}}$
${\displaystyle (\mu ^{*}-K){\frac {\mu +\zeta }{\mu ^{*}+\zeta }}=\sum _{i}^{M}f_{i}(\mu _{i}-\mu )Q_{i}}$
Where

${\displaystyle \zeta ={\frac {\mu (9K+8\mu )}{6(K+2\mu )}}}$

Where Pi and Qi for the ith inclusion with aspect ratio α given by,

Shape Pi Qi
Sphere ${\displaystyle {\frac {K+{\tfrac {4}{3}}\mu }{K_{i}+{\tfrac {4}{3}}\mu }}}$ ${\displaystyle {\frac {\mu +\zeta }{\mu _{i}+\zeta }}}$
Needle ${\displaystyle {\frac {K+\mu +{\tfrac {1}{3}}\mu _{i}}{K_{i}+\mu +{\tfrac {1}{3}}\mu _{i}}}}$ ${\displaystyle {\frac {1}{5}}({\frac {4\mu }{\mu +\mu _{i}}}+2{\frac {\mu +\gamma }{\mu _{i}+\gamma }}+{\frac {K_{i}+{\tfrac {4}{3}}\mu }{K_{i}+\mu +{\tfrac {1}{3}}\mu _{i}}})}$
Disk ${\displaystyle {\frac {K+{\tfrac {4}{3}}\mu _{i}}{K_{i}+{\tfrac {4}{3}}\mu _{i}}}}$ ${\displaystyle {\frac {\mu +\zeta _{i}}{\mu _{i}+\zeta _{i}}}}$
Penny Crack ${\displaystyle {\frac {K+{\tfrac {4}{3}}\mu _{i}}{K_{i}+{\tfrac {4}{3}}\mu _{i}+\pi \alpha \beta }}}$ ${\displaystyle {\frac {1}{5}}(1+{\frac {8\mu }{4\mu _{i}+\pi \alpha (\mu +2\beta )}}+2{\frac {K_{i}+{\tfrac {2}{3}}(\mu _{i}+\mu )}{K_{i}+{\tfrac {4}{3}}\mu _{i}+\pi \alpha \beta }})}$

Where ${\textstyle \beta =\mu {\frac {3K+\mu }{3K+4\mu }},\gamma =\mu {\frac {3K+\mu }{3K+7\mu }},\zeta ={\frac {\mu (9K+8\mu )}{6(K+2\mu )}}}$

It's important to note that for a solid matrix material, this formulations assumes pores are small and isolated and therefore there is no fluid flow due to wave propagation. This is equivalent to assuming the material is observed by high frequency waves. For use at low frequencies, Mavko et. al [6] recommend Toksöz formulation to obtain dry frame moduli for the material, then obtain the fluid saturated moduli using the Gassmann Equation, which is valid only at low frequencies.

If the matrix is a fluid, the problem becomes dynamic, as wave disturbances will cause relative motion between the fluid and the inclusions accounted for in a dynamic density term.

${\displaystyle {\frac {\rho -\rho ^{*}}{\rho +2\rho ^{*}}}=f{\frac {\rho -\rho '}{\rho +2\rho '}}}$

This is not the same as the usual mass balance formula as it accounts for mass associated with the coupled solid-fluid motion that occurs in the dynamic case. µ*=0 as a material requires a solid continuum to sustain shear stress, consequently the effective properties have no dependence on the shear modulus of the inclusions. K* can be obtained as above, setting µ to zero.

This formulation emphasizes the effect pore shape has on observed seismic properties. Flatter, lower aspect ratio pores are shown to have a much larger effect effective moduli than more spherical pores. Consequently observed fluid effects are strongly dependent on the pore shape. Because spherical inclusions have such a small effect on the effective moduli, transitioning from water filled to gas filled pores will only significantly change the density, causing seismic velocities to increase. However if the pores are flat, the effect of the inclusions becomes significant, and therefore the decrease in fluid bulk modulus when transitioning from water to gas is greater than the density decrease, so seismic velocities decrease.

### Self-consistent approximations

Self-consistent approximations attempt to expand effective moduli formulations beyond small inclusion constraints by considering the inclusion of a material within a background matrix with the same properties as the effective medium. The moduli of background medium are free to vary until they match the moduli of the composite, meaning that the inclusion is invisible to incoming waves and thus homogenization of the material has been achieved.

### Differential Effective Medium Model (DEM)

In the Differential Effective Medium model, the effective moduli of a composite are obtained by incrementally replacing small volumes of a background medium with small volumes of inclusions. This is done by considering Hooke's Law, averaged over a representative volume, Vo. An infinitesimal volume of background matrix, ∂v is removed and replaced with an infinitesimal volume of inclusion material, ∂vi[7]

${\displaystyle V_{0}\mathbf {\bar {\sigma }} =V_{0}(\mathbf {c} +\Delta \mathbf {c} )\mathbf {\bar {\varepsilon }} =\int _{V_{0}-\Delta v}\mathbf {c\varepsilon } dv+\int _{\Delta v_{i}}\mathbf {c} _{i}\mathbf {\varepsilon } dv}$
${\displaystyle \Delta \mathbf {c{\bar {\varepsilon }}} ={\frac {1}{V_{0}}}\int _{\Delta v_{i}}(\mathbf {c} _{i}-\mathbf {c} )\mathbf {\varepsilon } dv}$
Eshelby[8] found, using energy methods that the strain associated with an arbitrary inclusion is given by ${\textstyle \varepsilon =T\varepsilon _{0}}$where T is Wu's inclusion tensor. Thus in the infinitesimal limit,

${\displaystyle d\mathbf {c} (f)=(\mathbf {c_{i}} -\mathbf {c} )\mathbf {T} {\frac {df}{1-f}}}$
${\displaystyle \mathbf {c} (0)=\mathbf {c} _{matrix}}$

For an isotropic material with ellipsoidal inclusions this leads to

${\displaystyle {\frac {d}{df}}K(f)={\frac {(K'-K^{*})P'}{1-f}}}$
${\displaystyle {\frac {d}{df}}\mu (f)={\frac {(\mu '-\mu ^{*})Q'}{1-f}}}$
${\displaystyle K(0)=K_{matrix},\mu (0)=\mu _{matrix}}$

with P and Q as defined in the Kuster Toksöz section.

Note that this is a not a symmetric solution, depending on the order in which components are added. Solutions for multiple inclusions exist. It differs from self-consistent approximations in that a specific host material is identified and inclusions are added to it, whereas self-consistent methods continually modify the properties of the host to "absorb" the effect of adding inclusions. DEM moduli tend to be softer than Kuster Toksöz moduli but stiffer than self-consistent moduli.

## Important papers

• Ojala, I.O. Using rock physics for constructing synthetic sonic logs, ROCKENG09: Proceedings of the 3rd CANUS Rock Mechanics Symposium, Toronto, May 2009 (Ed: M.Diederichs and G.Grasselli), Paper 4016, 10p.
• Castagna, J, M Batzle, and T Kan (1993), Rock Physics - The Link Between Rock Properties and AVO Response. In Castagna, J and M Backus, eds. Offset-Dependent Reflectivity—Theory and Practice of AVO Analysis, 1993, p 135–171, DOI: 10.1190/1.9781560802624.
• Han, D, A Nur, and D Morgan (1986), Effects of porosity and clay content on wave velocities in sandstones. Geophysics, November 1986, p 2093–2107, DOI: 10.1190/1.1442062.
• Hashin, Z., and Shtrikman, S., 1963, A variational approach to the elastic behavior of multiphase minerals: J. Mech. Phys. Solids, 11, 127-140,DOI:10.1016/0022-5096(63)90060-7.

## References

1. Watt, P. J., G. F. Davies and R. J. O’Connell, The elastic properties of composite materials., Rev. Geophys. Space Phys. 14, 541-563, 1976. http://onlinelibrary.wiley.com/doi/10.1029/RG014i004p00541/full
2. Hill, R., Elastic properties of reinforced solids: Some theoretical principles, J. Mech. Phys. Solids, 11, 357-372, 1963. http://www.sciencedirect.com/science/article/pii/002250966390036X
3. Hashin, Z., and S. Shtrikman, A variational approach to the elastic behavior of multiphase materials, J. Mech. Phys. Solids, 11. 127-140, 1962.
4. Kuster, G.T., and Toksöz, M.N., 1974, Velocity and attenuation of seismic waves in two-phase media: Part I. theoretical formulations: Geophysics, 39, 587-606, DOI:10.1190/1.1440450
5. Wu, T. T., The Effect of Inclusion Shape on the Elastic Moduli of a Two-Phase Material., Int. J. Solid Structures. 2, 1-8, 1966. http://www.sciencedirect.com/science/article/pii/0020768366900023
6. Mavko, G, T. Mukerji and J. Dvorkin (2009), The Rock Physics Handbook: Cambridge University Press.
7. Norris, A.N., A Differential Scheme for the Effective Moduli of Composities, Mechanics of Materials, 4, 1-16, 1985. DOI: 10.1016/0167-6636(85)90002-X
8. Eshelby, J.D., The Determination of the Elastic Field of an Ellipsoidal Inclusion and Related Problems, Proc. Roy. Soc. London, Ser. A241, 376, 1957, DOI: 10.1098/rspa.1957.0133