Exact seismic velocities for TI media and extended Thomsen formulas for stronger anisotropies

THOMSEN'S WEAK ANISOTROPY FORMULATION FOR SEISMIC WAVES

Thomsen's weak anisotropy formulation (Thomsen, 1986), being a collection of approximations designed specifically for use in velocity analysis for exploration geophysics, is clearly not exact. Approximations incorporated into the formulas become most apparent for angles $\theta$ greater than about 15$^o$ from the vertical, especially for compressional and vertically polarized shear wave velocities $v_p(\theta)$ and $v_{sv}(\theta)$, respectively. For VTI media, angle $\theta$ is measured from the $\hat{z}$-vector pointing directly into the earth.

For reference purposes, I include here the exact velocity formulas for; quasi-P, quasi-SV, and SH seismic waves at all angles in a VTI elastic medium. These results are available in many places (Postma, 1955; Musgrave, 1959, 2003; Rüger, 2002; Thomsen, 2002), but were taken directly from Berryman (1979) with only some minor changes of notation; specifically, the $a$,$b$,$c$,$f$,$l$,$m$ notation for stiffnesses has been translated to the Voigt $c_{ij}$ stiffness notation wherein $a \to c_{11}$, $b \to c_{12}$, $c \to c_{33}$, $f \to c_{13}$, $l \to c_{44}$, and $m \to c_{66}$. The results are:

$$v_p^2(\theta) = \frac{1}{2\rho}\left\{\left[\left(c_{11}+c_{... ...t(c_{33}+c_{44}\right)\cos^2\theta\right] + R(\theta)\right\}$$ (A-1)

and

$$v_{sv}^2(\theta) = \frac{1}{2\rho}\left\{\left[\left(c_{11}+... ...(c_{33}+c_{44}\right)\cos^2\theta\right] - R(\theta)\right\},$$ (A-1)

where

$$R(\theta) = \sqrt{\left[\left(c_{11}-c_{44}\right)\sin^2\the... ...t]^2 + 4\left(c_{13}+c_{44}\right)^2\sin^2\theta\cos^2\theta}$$ (A-1)

and, finally,

$$v_{sh}^2(\theta) = \frac{1}{\rho}\left[c_{44} + (c_{66}-c_{44})\sin^2\theta\right].$$ (A-1)

I have purposely written equations 1 and 2 in this way to emphasize the fact that $v_p^2(\theta)$ and $v_{sv}^2(\theta)$ are closely related since they are actually the two solutions of a quadratic equation having the form:

$$\left(v^2\right)^2 - \left(v_p^2 + v_{sv}^2\right)v^2 + v_p^2v_{sv}^2= 0.$$ (A-1)

Any approximations made to one of these two wave speeds should therefore always be reflected in the other for this reason. In particular, any approximation to the square root in $R$ should be made consistently for both $v_p$ and $v_{sv}$.

For VTI symmetry, the stiffness matrix $c_{ij}$ is defined for $i,j = 1,\ldots,6$ by

$$c_{ij} = \left( \begin{array}{cccccc} c_{11} & c_{12} & c_{1... ...\ & & & & c_{44} & \\ & & & & & c_{66} \end{array}\right),$$ (A-1)

where $c_{12} = c_{11}-2c_{66}$. In an isotropic system (which is a more restrictive case than our current interests), $c_{12} = c_{13} = \lambda$, $c_{44} = c_{66} = \mu$, and $c_{11} = c_{33} = \lambda + 2\mu$, where $\lambda$ and $\mu$ are the usual Lamé constants. The definition in equation 6 makes use of the Voigt notation, i.e., $6\times6$ matrix instead of 4th order tensor, wherein Voigt single indices $i,j = 1,2,3,4,5,6$ correspond to the pairs of tensor indices 11,22,33,23,31,12, respectively. And it relates stress $\sigma_{ij}$ to strain $\epsilon_{ij}$ via $\sigma_{23} = c_{44}\epsilon_{23}$, $\sigma_{31} = c_{44}\epsilon_{31}$, $\sigma_{12} = c_{66}\epsilon_{12}$, and $\sigma_{ii} = \Sigma_j c_{ij}\epsilon_{jj}$ (no summation over repeated indices is assumed here) for $i,j = 1,2,3$. For VTI symmetry, we typically take $x_3 = z$ (the vertical) as the axis of symmetry. But, for HTI symmetry, we may choose index direction $x_3$ to be some other physical direction (such as horizontal directions $x$ or $y$, or some linear combination thereof); having done this, equations 2-4 apply strictly only in the vertical plane perpendicular to the fracture plane, while a small amount of vector analysis is then required to obtain the velocity values at all azimuthal angles $\phi \ne \pi/2$ away from the fracture plane.

Expressions for phase velocities in Thomsen's weak anisotropy limit can be found in many places, including Thomsen (1986, 2002) and Rüger (2002). The pertinent expressions for phase velocities in VTI media as a function of angle $\theta$, measured as previously mentioned from the vertical direction, are

$$v_p(\theta) \simeq v_p(0)\left(1 + \epsilon\sin^2\theta - (\epsilon-\delta)\sin^2\theta\cos^2\theta\right),$$ (A-1)

$$v_{sv}(\theta) \simeq v_s(0)\left(1 + \left[v^2_p(0)/v^2_s(0)\right](\epsilon - \delta)\sin^2\theta\cos^2\theta\right),$$ (A-1)

and

$$v_{sh}(\theta) \simeq v_s(0)\left(1 + \gamma\sin^2\theta\right).$$ (A-1)

In our present context, $v_s(0) = \sqrt{c_{44}/\rho_0}$, and $v_p(0) = \sqrt{c_{33}/\rho_0}$, where $c_{33}$, $c_{44}$, and $\rho_0$ are two stiffnesses of the cracked medium and the mass density of the isotropic host elastic medium. [For the specific physical examples that follow involving models of fractured reservoirs, I assume that the cracks contain insufficient volume to affect the overall mass density significantly.] The three Thomsen (1986) seismic parameters appearing in equations 7-9 for weak anisotropy with VTI symmetry are $\gamma = (c_{66}-c_{44})/2c_{44}$, $\epsilon = (c_{11}-c_{33})/2c_{33}$, and

$$\delta = \frac{(c_{13}+c_{44})^2-(c_{33}-c_{44})^2}{2c_{33}(... ...ht) \left(\frac{c_{13}+2c_{44}-c_{33}}{c_{33}-c_{44}}\right).$$ (A-1)

Parameter $\gamma$ is a measure of the shear wave anisotropy and birefringence. Parameter $\epsilon$ is a measure of the quasi-P wave anisotropy. Parameter $\delta$ controls the complexity of the shape of the wave fronts for quasi-P and quasi-SV waves; e.g., when $\delta = \epsilon$ the wave fronts are elliptical in shape, whereas for all TI anisotropic systems having $\epsilon - \delta \ne 0$, the wave front will deviate from being elliptical, and it is in such cases that ray arrival triplications may occur.

All three of these parameters $\gamma$, $\epsilon$, $\delta$ can play important roles in the velocities given by equations 7-9 when the anisotropy is large, as would be the case in fractured reservoirs when the crack densities are high enough. If crack densities are very low, then the SV shear wave will actually have no dependence on angle of wave propagation. Note that the so-called anellipticity parameter (Dellinger et al., 1993; Fomel, 2004; Tsvankin, 2005, p. 253), ${A} = \epsilon - \delta$, vanishes when $\epsilon \equiv \delta$ -- which (as will be shown) does happen to a very good approximation for low crack densities. Then, the results are anisotropic but have the special (elliptical) shape to the wave front mentioned previously.

For each of these phase velocities, the derivation of Thomsen's approximation has included a step that removes the square on the left-hand side of equations 1, 2, and 4 -- obtained by expanding a square root of the right hand side. This step introduces a factor of $\frac{1}{2}$ multiplying the $\sin^2\theta$ terms on the right hand side, and -- for example -- immediately explains how equation 8 is obtained from equation 4. The other two equations for $v_p(\theta)$ and $v_{sv}(\theta)$, i.e., equations 7 and 8, involve additional approximations. More of the details about the nature of these approximations are elucidated by first obtaining an alternative approximate formulation.

Exact seismic velocities for TI media and extended Thomsen formulas for stronger anisotropies