Orthorhombic anisotropy: A physical seismic modeling study

Scott P. Cheadle∗, R. James Brown‡ , and Don C. Lawton‡

ABSTRACT and +40 percent for "sand-, shale-, and carbonate-prone"

units in the Prudhoe Bay area.
An industrial laminate, Phenolic CE, is shown to possess
seismic anisotropy. This material is composed of Physical seismic modeling can be extremely useful in
laminated sheets of canvas fabric, with an approximately bridging the gap between theory and the complexities
orthogonal weave of fibers, bonded with phenolic resin. It observed in field seismic data, and this seems particularly
is currently being used in scaled physical modeling studies true in the context of seismic anisotropy. Many theoretical
of anisotropic media at The University of Calgary. predictions of wave-propagation phenomena can be tested
Ultrasonic transmission experiments using this material in sealed laboratory experiments. Ultrasonic modeling
show a directional variation of compressional- and shear- using phenolic laminate is ideally suited to the study of
wave velocities and distinct shear-wave birefringence, or velocity anisotropy because the ambiguities inherent in
splitting. Analysis of group-velocity measurements taken field data are absent. Tatham et al. (1987), Sayers (1988),
for specific directions of propagation through the material Ebrom et al. (1990) and Rathore et al. (1990) have
demonstrates that the observed anisotropy is characteristic described physical modeling experiments simulating
of orthorhombic symmetry, i.e., that the material has three fracture-induced anisotropy. This type of work is being
mutually orthogonal axes of two-fold symmetry. For P carried further within the CREWES Project (Consortium
waves, the observed anisotropy in symmetry planes varies for Research in Elastic Wave Exploration Seismology) at
from 6.3 to 22.4 percent, while for S waves it is observed to The University of Calgary.
vary from 3.5 to 9.6 percent.
This paper describes the results of experiments to
From the Kelvin-Christoffel equations, which yield phase determine the anisotropic elastic properties of Phenolic CE.
velocities given a set of stiffness values, expressions are We first determine the variation of the body-wave group
elaborated that yield the stiffnesses of a material given a velocities with direction by measuring traveltimes over a
specified set of group-velocity observations, at least three selection of paths. The theory of wave propagation in
of which must be for off-symmetry directions. anisotropic media (Appendix) is then used to relate these
observed group velocities to the nine (in the orthorhombic
INTRODUCTION
case) elastic coefficients or stiffnesses, permitting us to
Studies of anisotropy and shear-wave splitting are gaining compute the details of elastic-wave propagation in any
importance as part of the ongoing effort to enhance seismic direction through the phenolic and, in particular, to
data interpretation and reservoir exploitation. Several compute the variation of quantities such as Thomsen’s
authors have dealt with the relationships among anisotropy, (1986) anisotropy parameters and the phase velocity.
shear-wave polarization and fracture patterns (e.g. Keith
PHYSICAL MODEL EXPERIMENTS
and Crampin, 1977; Crampin, 1981, 1984, 1985; Lewis et
al., 1989; Yale and Sprunt, 1989). Liu et al. (1989) used Laboratory set-up
numerical modeling results to outline the potential and We are using piezoelectric P-wave and S-wave transducers
limitations of shear-wave splitting analysis for the as both sources and receivers in our multicomponent
crosswell configuration. Both compressional- and shear- physical modeling. Both types are flat-faced cylindrical
wave anisotropy impact on velocity analysis for contact transducers with an active element 12.6 mm in
multicomponent seismic imaging and on methods of diameter. With reference to a horizontal profile, the
estimating subsurface stress based on the VP/VS ratio compressional or P-wave transducer (Panametrics V103) is
(Thomsen, 1986, 1988). Banik (1984) reported errors in vertically polarized, with the maximum sensitivity normal
depth estimates of between 150 and 300 m in areas of the to the contact face; the shear-wave transducer (Panametrics
North Sea basin due to anisotropy within some shaly units. V153) is horizontally polarized, with the maximum
Ensley (1989) described anisotropy values of between -40 sensitivity parallel to a line across the contact face. During
operation, these contact faces are coupled to a selected flat

Presented at the 60th Annual International Meeting, Society of Exploration Geophysicists. Manuscript received by the Editor August 24,
1990; revised manuscript received March 21, 1991.

∗ Formerly Department of Geology & Geophysics, The University of Calgary; presently Veritas Seismic Ltd., 200, 615-3rd Ave. SW
Calgary, Alberta, Canada T2P OG6.
‡ Department of Geology & Geophysics, The University of Calgary, Calgary, Alberta, Canada T2N ]N4.

1
 Orthorhombic Physical Seismic Modeling

surface of the phenolic and, for a particular experiment, a Experimental results

profile direction and sagittal plane are established. To Shear-wave splitting experiments were conducted using
record the radial component, the shear-receiver transducer cubes of the phenolic as described above. Figures 3, 4, and
is used with its polarization parallel to the direction of the 5 show the transmission records through Faces 1, 2, and 3,
profile (inline), whereas for the transverse component, it is respectively, of an approximately 9.6 cm cube of phenolic.
rotated so that the polarization is perpendicular to the Each trace records the signal transmitted through the cube
azimuth of the profile and to the sagittal plane (cross-line). at 5-degree intervals of rotation with respect to the
polarization direction of the shear-wave transducers. The
The source transducer is driven with a 28-volt square wave
0-degree direction was chosen to correspond to the
tuned to produce a broadband wavelet with a central
polarization azimuth of the amplitude maximum of the
frequency of 600 kHz. Amplified data are sampled using a
faster of the two shear-wave arrivals. The sample interval
Nicolet digital oscilloscope connected through an IBM-XT,
used in this study was 50 nanoseconds, and the arrival
which controls the experiments, to a Perkin-Elmer 3240
times are shown in microseconds. The faster shear arrival
seismic processing system for storage. Traces of up to
is designated S1 and the slower mode S2. While it is more
4096 samples are recorded sequentially and stored on tape
correct to refer to the split shear waves and the
or disk in SEG-Y format.
compressional waves under most conditions as quasishear
The CE-grade phenolic laminate is composed of layers of a and quasicompressional modes, except for special cases,
woven canvas fabric saturated and bonded with a phenolic such as propagation in one of the principal directions, that
resin, and has a density of 1364 kg/m3. In one direction of prefix will usually be implied rather than explicitly stated.
the fabric the fibers of the warp run more or less straight, [Crampin (1989) and Winterstein (1990) have provided
like the fixed threads on a loom; in the orthogonal direction authoritative manuals of terminology for seismic
the fibers of the woof run back and forth across the warp. anisotropy.] In Figures 3, 4, and 5, the weakly coupled P-
Initial tests with the material showed a directional wave arrival is barely visible. The compressional velocities
dependence of the velocity for both P and S waves, were determined separately using the P-wave transducers.
suggesting its suitability for physical modeling of an
In Figures 3 and 4, for propagation in the 1- and 2-
anisotropic medium. Shear-wave splitting was observed
directions, respectively, the polarization (particle motion) at
during transmission tests when the sample was rotated
the S1 amplitude maximum is, in each case, parallel to the
between two shear-wave transducers. The polarizations of
"bedding plane" of the canvas layers, whereas for the S2
the split shear waves were approximately parallel to the
amplitude maxima, the polarizations are perpendicular to
orientations of the orthogonal weave of fibers in the canvas
this plane. In Figure 5, for propagation in the 3-direction,
fabric. For this reason, subsequent experiments were
the polarization of the S1 amplitude maximum is parallel to
conducted on pieces of phenolic that were cut with faces
the 1-direction (the woof), while the S2 amplitude
parallel or orthogonal to the observed fiber directions as
maximum is parallel to the 2-direction (the warp). A plot
well as to the plane of the canvas layers. A sample of the
of amplitude versus polarization direction for a record
phenolic with the faces labeled with the convention used in
this study is shown in Figure 1. The factory-machined
surface of the laminate sheet parallel to the fabric layers
was designated Face 3, consistent with the conventional
choice of x3 as the vertical direction and with a horizontal
attitude for the layering of the medium. Since the 3-
direction turned out to be slowest for P-wave propagation,
the other two principal (or symmetry) directions were
labeled such that the 1-direction (parallel to the woof) is
fastest and the 2-direction (parallel to the warp)
intermediate for P-wave propagation.
The apparatus used for studying split shear waves is shown
in Figure 2. The cube of material is placed between two
fixed shear-wave transducers which are aligned with
parallel polarizations. The cube is rotated between the
transducers. and a pointer on the cube is used to determine
the azimuth of the sample with respect to a fixed circular
protractor. A similar experimental procedure was described
by Tatham et al. (1987) for a study of fracture-induced FIG. 1. Phenolic CE laminate is composed of layers of a canvas
shear-wave splitting. weave fabric bonded with phenolic resin. The faces are labeled as
used in this study.

2
FIG. 2. The apparatus used for shear-wave splitting experiments clamps the cube of phenolic between two shear-wave transducers with mutually
parallel polarization. The cube is rotated while the transducers remain fixed. A circular protractor scale is used to determine the azimuth of
rotation.

FIG. 3. The record through Face I of a 9.6 cm cube of phenolic, FIG. 4. The record through Face 2, showing the faster S I arrival at
showing the faster S1 arrival at 1665 m/s and the slower S2 arrival 1658 m/s and the slower S, arrival at 1506 m/s. The compressional
at 1602 n-l/s. The compressional velocity in the 1-direction is velocity in the 2-direction is 3365 m/s. The polarization directions
3576 m/s. The polarization direction of the S 1 amplitude of the S1 and S2 amplitude maxima are parallel and perpendicular
maximum is parallel to the "bedding plane" of the canvas layers, respectively to the canvas layering, as in Figure 3.
while that of the S2 amplitude maximum is perpendicular to that
plane.

3
 Orthorhombic Physical Seismic Modeling

through Face 2 is shown in Figure 6. This and other defined on the basis of the P-wave velocities (Figure 7).
transmission records through the phenolic show that the S1 The values quoted are group velocities based on the transit
mode generally has a greater maximum amplitude than the time measured with respect to the onset of the pulse. The
S2 arrival, indicating greater attenuation for the S2 mode. velocities are the averages of values measured through 10-
The ratios of the amplitudes of the S1 arrivals to those of and 8-cm cubes. The measured velocities for the phenolic
the S2 arrivals, measured at their maxima, have ranged from cubes were repeatable to within ±15 m/s (≈0.5 percent) for
1.1 to 1.4 for the samples tested. P-waves and ±4 m/s (≈0.25 percent) for shear waves. The
variations are likely related to small inconsistencies in the
thickness of the coupling agent used to bond the
transducers to the phenolic. Velocity variations between
different samples of phenolic ranged up to 2 percent. The
time picks used to calculate the velocities were made
directly on the digital oscilloscope for maximum accuracy.
For the following discussion, the velocities will be labeled
with 2 subscripts indicating, respectively, the directions of
propagation and polarization with respect to the three
symmetry axes (Figure 7). For example, V11 is the group
velocity for propagation and particle motion in the 1-
direction (a P wave) while V12 indicates propagation in the
1-direction with polarization in the 2-direction (an S wave).
The six shear-wave velocities measured in the principal
directions were paired as follows: V23≈V32; V31≈V13;
V12≈V21; indicating, with very small error (Table 1), only
three independent values.

FIG. 5. The record through Face 3, showing the faster S1 (1610 PRINCIPAL AXES
m/s) and slower S2 (1525 m/s) shear waves. The compressional
velocity in the 3-direction, determined separately with P-wave RAY (GROUP) VELOCITIES
transducers, is 2925 m/s. The traces for the records of Figures 3 to
5 were recorded at 5-degree intervals of rotation.

FIG. 6. The plot of amplitude versus azimuth with respect to the

polarization direction of the shear-wave transducers for a record FIG. 7. The P- and S-wave velocities measured along the principal
through bv Face 2. The scatter of the measured amplitudes from axes are summarized, with the heavy arrow designating the
the sinusoidal variation with azimuth is due to variable coupling of direction of propagation and the lighter arrow the direction of
the transducers to the sample during rotation. particle motion of the shear waves. The subscripts correspond to
the directions of propagation and particle motion, respectively. Of
the six shear-wave velocities, three distinct pairs of values are
recognized.
The P, S1, and S2 velocities measured along the principal
axes are summarized in Figure 7 and those along the 45-
degree diagonals in Figure 8. Fast, medium, and slow For the cases of diagonal raypaths in symmetry planes we
directions through the cube (1, 2, and 3, respectively) were adopt, for the purpose of this paper alone, a special index

4
 Orthorhombic Physical Seismic Modeling

convention (Figure 8). For the 23-plane, the direction at 45 C11 C12 C13 
degrees to the 2- and 3-directions is denoted by the index 4.  C22 C23 
The group velocitv of the quasi-P (qP) wave in this   (1)
 C33 
direction is thus designated V44. Polarization quasi-normal cmn = 
to this 4-direction but still within the 23-plane is denoted by  C44 
the index 4 . Thus the quasi-SV (qSV) velocity is designated  C55 
 
V4 4 . The velocity of the corresponding SH wave, with  C66 
particle motion in the 1-direction, is labelled V41.
Similarly, we use the indices 5 and 6 to denote propagation of nine independent coefficients (e.g. Nye, 1985). Using
in the 31 - and 12-planes, respectively, at 45 degrees. The the elastic equations of motion the stiffnesses Cmn may be
P-, SV-, and SH-wave group velocities are thus labeled V55, estimated from the observed body-wave velocities and the
density of the phenolic (see the Appendix). The computed
V5 5 , and V52, in the 31-plane, and V66, V , and V63, in the
66 stiffnesses are summarized in Table 1. Along the principal
12-plane (Figure 8), in each instance only for the special axes the phase and group velocities are equal and the
cases of rays at 45 degrees to the symmetry directions. stiffnesses were computed directly using equations (A-45)
Each of the velocities along the diagonal raypaths is the and (A-46). Along the diagonal raypaths, the direction of
average of two measurements (between the two pairs of the wavefront normal (i.e. the slowness direction) is not, in
opposing edges of the cube) which had equivalent raypaths general, the same as the 45-degree direction of the raypath
relative to the principal axes within each of the three (i.e., of energy transport). The procedure used to compute
principal planes. The two traveltimes for each of the the slowness directions, the phase velocities, and the related
diagonal raypath pairs were virtually identical, differing by stiffnesses for the diagonal raypaths is described in the
two sample points (100 ms) or less in all cases. Four Appendix.
measurements were also recorded for raypaths from corner Nine independent velocity values, are required to enable
to corner of the cube, with similarly small differences in the complete determination of the stiffness matrix for the case
velocities. This symmetry confirmed that the presumed of orthorhombic anisotropy. These could include the three
principal planes, chosen to correspond to the planar P-wave velocities along the principal axes, three shear-
layering of the canvas fabric and the orthogonal weave of wave velocities (one of each pair, or their average) also
fibers in the phenolic, are indeed the seismic anisotropic along the principal axes, and three P-wave or SV-wave
symmetry planes. velocities, each for a raypath perpendicular to one and at 45
45° AXES degrees to the other two principal axes. In principal,
RAY (GROUP) VELOCITIES measurements at other orientations could be used but these
would require considerably more complex solutions.
Since we actually observe more than nine velocities, the
internal consistency of the orthorhombic symmetry model
can be checked. In addition to observations in this context
already mentioned above, equation (A-44) was used with
the shear-wave velocities observed along the principal axes
(Figure 7) to calculate additional independent values for the
SH-wave velocities along the diagonal raypaths (Figure 8).
For example,

(
V41′ = 2V13V12 / V132 + V122 )
12
(2)
= 1633m / s,
where V41, the observed value, is 1636 m/s, differing by
0.18 percent from V’41 the calculated value. Similarly, V’52
FIG. 8. The results of transmission measurements between
opposing edges of the phenolic cube are summarized. The =1583m/s, equal toV52, and V’63 = 1559m/s, differing by
propagation directions were at 45 degrees to two of the principal 0.l9 percent from the V63 value of 1556 m/s. Clearly, the
axes and perpendicular to the third. SH-mode velocities observed along the diagonal raypaths
conform very well to the assumed orthorhombic symmetry
ORTHORHOMBIC ANISOTROPY model.
For the orthorhombic symmetry system. the 3 x 3 x 3 x 3 The off-diagonal stiffnesses have been computed (Table 1)
stiffness tensor Cijkl (see the Appendix) may be reduced to a using the velocities measured on diagonal raypaths
6 x 6 symmetric matrix, namely: [equations (A-47)] both for P (giving Cij) and for SV
(giving cij ). The differences between pairs (Cij, cij ) seem

5
 Orthorhombic Physical Seismic Modeling

rather large at first; however, the stiffness values are quite The results of these experiments indicate that the
sensitive to small velocity changes. For each wave type we orthorhombic symmetry system is appropriate to describe
compare the two velocities (one measured and one the anisotropy of this material. In particular, three distinct
calculated) corresponding to a stiffness pair (Cij, cij ). The compressional velocities were measured in mutually
orthogonal directions consistent with the visible structural
following are the relative deviations from the mean of these symmetry of the material; the observed shear-wave
velocity pairs, with the corresponding deviation for the velocities along the principal axes are equal in pairs (Vij =
stiffness in square brackets (all given as percentages) -21- Vji); the observed diagonal SH velocities agree well with
plane: ±0.3 (P), ±1.0 (SV), [±2.5]; 13-plane: ±0.6 (P), ±2.4 those calculated independently [equation (2); Table 1]; and
(SV),[ ±5.7]; 23-plane: ±0.5 (P), ±2.0 (SV), [±4.4]. the two independent determinations of each off-diagonal
stiffness (C23, C31, C12) , from a qP and a qS measurement,
are in good agreement (Table 1).

Table 1. Body-wave velocities and phenolic stiffnesses

Mode Group Velocity (m/s) Phase Velocity (m/s) Phase Angle (deg) Stiffness (x109 N/m2)
Raypaths in Principal Directions
P v33 2925 C33 11.670

P V 22 3365 C 22 15.445

P V 11 3576 C 11 17.443

S V 21 1658 1662 C 66 3.768

S V 12 1665 (avg)

S V 31 1610 1606 C 55 3.518

S V 13 1602 (avg)

S V 32 1525 1516 C 44 3.135

S V 23 1506 (avg)

Raypaths at 45 Degrees to Principal Directions

P V 66 3378 3373 41.6 C21 7.104 7.283

SV V66 1810 1809 47.1 C21 7.462 (avg)

SH V 63 1556 (1559)∗ 1556 48.3

P V 55 3219 3155 33.4 C13 6.258 6.633

SV V55 1620 1620 45.0 C13 7.008 (avg)

SH V 52 1583 (1583) * 1577 39.8

P V 44 3094 3066 37.1 C23 6.097 6.379

SV V 44 1569 1569 45.7 C 23 6.662 (avg)

SH V41 1636 (1633)* 1632 47.0

∗
See equation (2)

6
 Orthorhombic Physical Seismic Modeling

DISCUSSION macroscopic, including preferred orientation of mineral

grains, pores or fractures (Crampin, 1981, 1984, 1985),
Degree of anisotropy thin-layer lamination (Helbig, 1983) and regional stress
The conventional measures of anisotropy for the transverse (Nikitin and Chesnokov, 1984). Anisotropy has been
isotropy case are given by Thomsen (1986) as recognized in many rocks (Banik, 1984; Thomsen, 1986;
ε ≈ [Vp (90 degrees) – Vp (0 degrees]/Vp (0 degrees) (3) Lewis et al., 1989; Ensley, 1989), but the physical cause
and symmetry systems of specific cases of anisotropic
and media are seldom unambiguously identified. Transverse
γ ≈ [VS (90 degrees) – VS (0 degrees]/VS (0 degrees) (4) isotropy can be invoked for horizontal thin-bed layering,
for example, in shale sequences, while azimuthal
anisotropy may arise in the idealized case of aligned
At least in the case of transverse isotropy, the ε parameter vertical fractures. Both of these examples would be
is not always useful in the context of the limited ray angles degenerate cases of the more general orthorhombic system.
typical of surface seismic gathers. Thomsen (1986) also Two or more sources of anisotropy superimposed
defined the parameter and discussed its use in conjunction orthogonally within the same lithologic unit, such as
with moveout-velocity and stress analysis. aligned vertical fracturing of a horizontally laminated
sequence, could result in orthorhombic anisotropy. The
Table 2. Measured anisotropy phenolic laminate is being used to simulate media with
similar velocity properties regardless of the different
δ ε γ
physical causes of the anisotropy.
21-plane -0.047 0.063 0.059
CONCLUSIONS
31-plane 0.183 0.224 0.096
Ultrasonic modeling with Phenolic CE laminate has
32-plane 0.081 0.150 0.035 demonstrated the anisotropic elastic properties of the
δ ≈ [Vp (45 degrees)/ Vp (0 degrees) – 1] material. The patterns of shear-wave splitting observed
– VP (90 degrees]/VP(0 degrees) – 1] (5) through each face of a cube of the phenolic, along with the
measured compressional-wave velocities, were used to
The measurements of these velocity ratios determined in define orthogonal principal (or symmetry) axes related to
the principal planes of the phenolic are shown in Table 2, the fast, medium, and slow directions through the material.
and fall within the range of the values reported by Thomsen Shear- and compressional-wave velocities were also
(1986) for a variety of rocks. The P-wave anisotropy ranges measured in directions between opposing edges of the cube
form 6.3 percent in the 21-plane to 22.4 percent in the 31- to support the determination of the orientations of the
plane, the SH-wave anisotropy from 3.5 percent in the 32- planes of symmetry. Within a principal plane, the SV wave
plane to 9.6 percent in the 31-plane. The plane of weakest has equal velocities for propagation in either of the axial
P anisotropy is not the same as the plane of weakest SH directions. Velocities computed (using equations of
anisotropy, but the planes of strongest anisotropy do propagation based on velocities from other directions and
correspond. Anisotropy of the SV mode is observed along assuming the orthorhombic model, closely matched the
the 45 degrees raypaths. In the 12-plane, V is 1810 m/s, observed values. Analysis of the data supports the
66

8.9 percent higher than V21. In the 32 plane, V is 1569 interpretation that the anisotropy conforms very closely to a
44 system of orthorhombic symmetry.
m/s, 3.5 percent higher than V32. SV anisotropy in the 31-
plane is lowest (0.9 percent) despite this plane exhibiting Physical modeling is currently proceeding with the
the strongest P and SH anisotropy. Although some of these phenolic and involves the recording of shot gathers as well
observations may seem surprising intuitively, they are all as simulated VSP and crosswell experiments. Observations
quite reasonable since the material has nine independent of the effect of orthorhombic anisotropy on moveout
stiffnesses. In the 23-plane for instance, P and SH velocities are being reported by the present authors (Brown
anisotropies depend on C22 and C33 (P) and C55 and C66 et al., 1991) and tomographic reconstruction will be
(SH) whereas the SV anisotropy depends not only on C22, described in future publications. Physical model data using
C33, and C44, but also on C23. phenolic laminate is proving to be a valuable adjunct to
numerical studies of the increasingly important topic of
Origin of the anisotropy seismic anisotropy.
The cause of the anisotropy in the phenolic laminate
appears to be related to the layering and the weave of the ACKNOWLEDGMENTS
canvas fabric. The material behaves like a stack of nets set This research is supported by funding provided by the
in a gel, with different fiber densities and orientations in the corporate sponsors of the CREWES Project. The authors
directions of the three principal axes. The many causes of wish to acknowledge the technical contributions of Mr.
anisotropy in rocks range from the microscopic to the Malcolm Bertram of the Department of Geology and

7
 Orthorhombic Physical Seismic Modeling

Geophysics and Mr. Eric Gallant of the CREWES Project, Liu, E., Crampin, S., and Booth, D. C., 1989, Shear-wave
both at The University of Calgary. Finally, we would like splitting in cross-hole surveys: Modeling: Geophysics, 54,
to thank Dan Ebrom and Linda Zimmerman for their 57-65.
knowledgeable reviews of this paper and for their
Musgrave, M. J. P., 1970, Crystal acoustics: Holden-Day.
constructive suggestions, which have led to an improved
final version. Nikitin, L. V., and Chesnokov. E. M., 1984. Wave
propagation in elastic media with stress-induced
REFERENCES anisotropy: Geophys. J. Roy. Astr,. Soc., 76, 129-133.
Aki, K., and Richards, P. G., 1980, Quantitative Nye, J. F., 1985, Physical properties of crystals: Oxford
seismology: Theory and methods: W. H. Freeman and Co. University Press.
Banik, N. C., 1984, Velocity anisotropy of shales and depth Rathore, J. S., Fjaer E., Holt, R. M., and Renlie, L., 1990.
estimation in the North Sea basin: Geophysics, 49, 1411- Acoustic anisotropy of synthetics with controlled crack
1419. geometries: Presented at the 4th Internat. Workshop on
Brown, R. J., Lawton, D. C., and Cheadle, S. P., 1991, Seismic Anisotropy.
Scaled physical modelling of anisotropic wave propagation: Sayers, C. M., 1988, Stress-induced ultrasonic; wave
multioffset profiles over an orthorhombic medium: velocity anisotropy in fracture rock: Ultrasonics, 26, 311-
Geophys. J. Internat., in press. 317.
Bullen, K. E., 1963, An introduction to the theory of Tatham, R. H., Matthews, M. D., Sekharan, K. K., Wade,
seismology: Cambridge University Press. C. J., and Liro, L. M., 1987, A physical model study of
Crampin, S., 1981, A review of wave motion in anisotropic shear-wave splitting and fracture intensity: 57th Ann.
and cracked elastic-media: Wave Motion, 3, 343-391. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts,
-- 1984, An introduction to wave propagation in anisotropic 642-645.
media: Geophys. J. Roy. Astr. Soc., 76, 17-28. Thomsen, L., 1986, Weak elastic anisotropy: Geophysics,
-- 1985, Evaluation of anisotropy by shear-wave splitting: 51, 1954-1966.
Geophysics, 50, 142-152.
-- 1989, Suggestions for a consistent terminology for ----- 1988, Reflection seismology over azimuthally
seismic anisotropy: Geophys. Prosp., 37, 753-770. anisotropic media: Geophysics, 53, 304-313.

Ebrom, D. A., Tatham, R. ll., Sekharan. K. K., McDonald, Vlaar, N. J., 1968. Ray theory for an anisotropic
J. A., and Gardner, G. H. F., 1990. Hyperbolic traveltime inhomogeneous elastic medium: Bull. Seis. Soc. Am., 58,
analysis of first arrivals in azimuthally anisotropic rnedium: 2053-2072.
A physical modeling study: Geophysics, 55, 185-191 Winterstein, D. F., 1990, Velocity anisotropy terminology
Ensley, R. A., 1989, Analysis of compressional- and shear- for geophysicists: Geophysics, 55, 1070-1088.
wave seismic data from the Prudhoe Bav Field: The leading Yale, D. P., and Sprunt, E. S., 1989, Prediction of fracture
Edge, 8, no. 11, l@13. direction using shear acoustic anisotropy: The Log Analyst,
Fedorov, F. I., 1968, Theory of elastic waves in crystals: 30, 65-70.
Plenum Press.
APPENDIX
Helbig. K., 1983, Elliptical anisotropy - Its significance and
meaning: Geophysics, 48, 825-832. RELATIONSHIPS AMONG VARIOUS ELASTIC-
WAVE PARAMETERS IN AN ANISOTROPIC
Keith, C. M., and Crampin, S., 1977, Seismic body waves MEDIUM OF ORTHORHOMBIC SYMMETRY
in anisotropic media: Reflection and refraction at a plane
interface: Geophys. J. Roy. Astr. Soc., 49, 181-208. Basic theory and the Kelvin-Christoffel equations
The equations of motion governing wave propagation in a
Kendall, J.M., and Thomson, C. J., 1989, A comment on
generally isotropic elastic medium are given by many
the form of the geometrical spreading equations, with some
authors (e.g. Bullen, 1963; Fedorov, 1968; Musgrave,
numerical examples of seismic ray tracing in
1970- Aki and Richards, 1980; Crampin, 1981, 1984; Nye,
inhomogeneous, anisotropic media: Geophys. J. Internat.,
1985) For infinitesimal displacements ui, Cartesian
99, 401-413.
coordinates xi, density ρ, stress tensor σij and body forces
Lewis, C., Davis, T. L., Vuillermoz, C., Gurch, M., per unit mass gi:
Iverson, W., and Schipperijn, A. P., 1989,Three-
dimensional imaging of reservoir heterogeneity, Silo Field, ρu&&i = σ ij , j + ρg i (A-1)
Wyoming: 59th Ann. Internat. Mtg., Soc. Expl. Geophys.,
Expanded Abstracts, 763-765.

8
 Orthorhombic Physical Seismic Modeling

where ",j" denotes the partial derivative with respect to xj will be real. (Throughout this appendix vertical bars denote
and where the Einstein summation convention (for repeated
determinant.)
indices) applies.
A further consequence of the symmetries embodied in (A-
The stress tensor, in terms of the strain tensor εkl and the 8) is that there are only 21 independent stiffnesses, Cijkl.
stiffness tensor Cijkl, is given in accordance with Hooke’s Following, e.g., Musgrave (1970), Nye (1985), and
law by: Thomsen (1986), the fourth-order stiffness tensor may be
σ ij = C ijkl + ε kl , (A-2) written as a second-order (6 x 6) symmetric matrix:

where C ijkl → C mn

ε kl = 1
2
(u l ,k + u k ,l ) . (A-3) where

Substituting (A-2) and (A-3) into (A-1), neglecting any

m =1 if i = j,  (A-10)

body forces, yields: m = 9 − (i + j ) if i ≠ j 
C ijkl u k , lj − ρ u&&i = 0 . (A-4) and n and kl are analogous to m and ij.

These equations of motion, and their solution for By introducing the so called Kelvin-Christoffel stiffnesses,
monochromatic plane-wave motion, are considered by given by Musgrave (1970) as:
many authors (e.g., Fedorov, 1968; Musgrave, 1970; Keith Γ ik = C ijkl n j n l (A-11)
and Crampin, 1977; Aki and Richards, 1980; Crampin,
1981, 1984) but here we follow Musgrave’s treatment most
closely. We assume harmonic plane-wave displacement,
expressed as equations (A-7) and (A-9) may be rewritten as:

u k = Ap k exp [i ω (s r x r − t )] , (A-5)
 Γ11 − ρv 2

Γ12 Γ13   p1 
 (A-12)
 Γ21 Γ22 − ρv 2 Γ23   p 2  = 0
where A is the amplitude factor, pk is the unit polarization  Γ31 Γ32 2
Γ33 − ρv   p 3 

(or particle displacement) vector, ω is angular frequency, sr
is the slowness vector, and in this equation only, i = − 1 . and
The slowness vector gives the direction of the wavefront  Γ11 − ρ v 2 
Γ12 Γ13
normal and may further be written:   (A-13)
 Γ21 Γ22 − ρ v 2 Γ23  = 0
s r = v −1 n r , (A-6)  Γ31
 Γ32 Γ33 − ρ v 2 

where v is phase velocity and nr is the unit slowness (or Equations (A-12) and (A-13) are known as the Kelvin-
wavefront-normal) vector. From equation (A-4), (A-5), Christoffel equations.
and (A-6) one obtains:
In the case of orthorhombic symmetry only nine of the, in
(C ijkl n j nl − ρv δ ik p k = 0
2
) . (A-7) general, 21 indepent stiffnesses, Cmn, are nonzero. The six
independent Kelvin-Christoffel stiffnesses are then
Thus, the determination of the details of the wave motion Γ 11 = n 12 C 11 + n 22 C 66 + n 32 C 55
has been cast as an eigenvalue problem in which, having Γ 22 = n 12 C 66 + n 22 C 22 + n 32 C 44
specified Cijkl (the stiffnesses of the medium) and nr (the (A-14)
direction of phase propagation), one can solve for pk (the Γ 33 = n 12 C 55 + n 22 C 44 + n 32 C 33

particle motion or polarization vector) and three values (in Γ 23 = n 2 n 3 (C 23 + C 44 )

general) for v (phase velocity). Γ 31 = n 3 n 1 (C 31 + C 55 )
Due to the well known symmetries involved (see e.g., Γ 12 = n 1 n 2 (C 12 + C 66 ).
Musgrave, 1970; Nye, 1985):
Propagation along a principal direction
C ijkl = C ijlk = C jikl = C klji (A-8)
For slowness vector in the 1-direction,
and therefore the matrix (Cijkl nj nl - ρv σik) is symmetric. = (1 , 0 , 0 )
2
n j
(A-15)
This implies in turn that the three eigenvalues obtained for
ρv2 by setting and equations (A-14) reduce to:
C ijkl n j n l − ρ v 2 σ ik = 0 (A-9)

9
 Orthorhombic Physical Seismic Modeling

Γ11 = C 11  (phase) with V (group) in equations (A-21), (A-22), and (A-

Γ 22 = C 66  (A-16) 23).


Γ 33 = C 55  Propagation of 45 degrees to two principal axes or
Γ 23 = Γ 31 = Γ12 = 0 . “edge to edge”
Equation (A-12) then becomes: Equation (A-21) to (A-23) and their cyclically varied
analogs allow one to determine the six stiffnesses along the
Γ11 − ρv 2
0 0   p1  diagonal of the Cmn matrix from velocities measured along
   (A-17)
 0 Γ22 − ρv 2 0  p2  = 0 principal directions. In order to determine the three
 0 0 2
Γ33 − ρv   p 3  independent off diagonal stiffnesses, one must measure

velocities for raypaths along different directions. The next
For this rather simple case, that of propagation along a simplest directions to consider would seem to be those in
principal direction, there are three obvious eigenvalues principal planes at 45 degrees to each of two principal
which will zero the determinant of the 3 x 3 matrix. For directions. We have measured velocities along each such
each of these, the associated eigenvector pk is the raypath for the three different polarization.
polarization (or unit-particle-displacement) vector.
Unfortunately, the raypath or group-velocity direction is
The P wave. – Choosing the eigenvalue solution: not, in general, the same as the wavefront-normal or phase-
velocity direction (Figure A-1). So we cannot make simple
Γ 11 − ρ v = 0 (A-18)
( )
2
substitutions for ni = 0 or 2 / 2 in equations (A-14).
reduces the three equation of (A-17) to two, namely: We need additional equations that will allow determination
 C 66 − C 11 0  p2  of ni and v from knowledge of ξi ( the unit vector in the
=0
. (A-19)
 C 55 − C 11   p 3 
group-velocity direction; Figure A-1)
 0
Such theory has been dealt with in several works (e.g.
The only permissable solution to (A-19) is Vlaar, 1968; Musgrave, 1970; Kendall and Thomson,
1989). Here we take the result from Musgrave (1970) and
p = p3 = 0 , (A-20)
2
refer the reader to the works cited for details. Starting from
since otherwise at least two of the six independent the geometrical relationship (Figure A-1):
stiffnesses would have to be equal, violating the
v V = cos ∆ = n i ξ i , (A-24)
assumption of orthorhombic symmetry. It follows from
equations (A-16), (A-18), and (A-20) that Musgrave (1970, p.89) gives:
p k = (1,0,0 ) v11 = (C11 / ρ )
12
and (A-21) 1  p2   2 ∂α 2 ∂A  ,
vi =  2  ρ v − A +α 2 
(A-25)
where v11 denotes that v which applies for propagation 2 ρv α  k  ∂ni ∂ni  k
(slowness) in the 1-direction with particle motion
(polarization) in the 1-direction, that is, the P-wave Where
velocity. 12
 Γik Γ jk  i≠ j≠k
The S waves.-Choosing each of the other two eigenvalue αk =   , (A-26)
solutions leads to the two solutions:  Γij 

p k = (0,1,0) v12 = (C 66 / ρ ) and

12
and (A-22)

and
Ak = Γkk (no summation) . (A-27)

pk = (0,0,1) v13 = (C55 / ρ ) , In equation (A-25), p, α, and A inside brackets should be

12
and (A-23)
represented by their kth components and the products of the
these representing S waves polarized in the 2- and 3- brackets are summed. This notation follows Musgrave
directions, respectively. (1970) except that, whereas Musgrave use ξi as the group-
velocity vector, we use Vi for group velocity and ξi as its
The corresponding velocities and polarizations for the
unit vector.
propagation in the 2- and 3- directions are obtained from
equations (A-21), (A-22), and(A-23) by cyclic variation of For a medium of orthorhombic symmetry, we get from
the indices (1,2,3) and the indices (4,5,6). Since, for these equations (A-26), (A-27), and (A-14):
axial propagation directions, the wavefront normal and the
raypath have the same direction, one could replace v

10
 Orthorhombic Physical Seismic Modeling

(C12 + C66 )(C31 + C55 ) p1 = 0 (A-33)

α 12 = n12 (A-28)
(C 23 + C 44 ) or
and ρ v 2 = n 22 C 66 + n 32 C 55 . (A-34)
A1 = n C11 + n C 66 + n C 55 ,
2
1
2
2
2
3
(A-29) The quasi-P and –SV waves. – Equation (A-33) implies
polarization entirely within the 23-place (the sagittal plane),
and similarly for k=2 and 3. Substitution into (A-25) results i.e., P-SV types. From the analogs to equation (A-30) for
in: i=2 and 3 we then get:
ρ vV 2 = p 22 n 2− 1 (ρ v 2 − n 32 C 44 ) + p 32 n 2 C 44
and (A-35)

ρ vV 3 = p 22 n 3 C 44 + p 32 n 3− 1 (ρ v 2 − n 22 C 44 ) .
Further, from equation (A-12) for this case we get:

p2 n n (C + C 44 ) ρv 2 − n 22 C 44 − n32 C 33
= 2 2 3 2 23 =
p 3 ρv − n2 C 22 − n32 C 44 n2 n3 (C 23 + C 44 )

(A-36)
from which

p 22 ρ v 2 − n 32 C 33 − n 22 C 44
= . (A-37)
p 32 ρ v 2 − n 22 C 22 − n 32 C 44
since the raypath or group-velocity direction is at 45
degrees to the 2- and 3-axes, V2 = V3 and the two right-hand
sides of equations (A-35) are equal. From this and equation
(A-37) one can eliminate p2 and p3 and obtain:

ρ 2 v 4 (n3 − n 2 ) + ρv 2 [n23 (C 22 + C 44 ) − n33 (C 33 + C 44 )]

+ C 33 C 44 n34 (n2 + n3 ) − C 22 C 44 n24 (n2 + n3 ) = 0
(A-38)
It is clear from (A-36) that there are two solutions for v2
Figure A-1. Schematic diagrams of a wavefront in an anisotropic and thus for p2/p3. One of these solutions is the quasi-P or
medium at times t and δt, showing the phase and group velocities, qP-wave case, the other the qSV-wave case. For qP we
v and V respectively, their respective unit vectors (and directions), denotes the phase velocity v44 and the group velocity V44.
ni and ξi, and the angle ∆ between ni and ξi. For qSV these are v and V , respectively. There is, it is
44 44
true, a fundamental incongruity between the single- and
ρvV1 = p12 [(ρv 2 − n12 C11 − n22 C 66 − n32 C55 )n1−1 + n1C11 ] double-subscript notations for V. However, we do not try
(
+ n1 p 22 C 66 + p32 C 55 ) to combine the two and thus no problem ever arises here.

(A-30) Defining θ=cos-1n3, we can write (figure A-1) that ∆=θ-45

degrees. Therefore, from equation (A-24) for qP:
( )
and similarly for i=2 and 3. If we consider propagation in
the 23-plane of symmetry: v ij = V ij cos ∆ = 2 2 (n 2 + n 3 )V ij , (A-39)

n 1 = 0 and V 1 = 0 . (A-31) where ij=44, 44 or 41, for qP,qSV, or SH, respectively.

Using equations (A-21) to (A-23) and their cyclic variants
Thus, from equation (A-30)
to eliminate stiffnesses, and (A-39) to eliminate v from (A-
( )
p 12 n 1− 1 ρ v 2 − n 22 C 66 − n 32 C 55 = 0 , (A-32) 38), we obtain:

and therefore either

11
 Orthorhombic Physical Seismic Modeling

( ) [ (
V444 (n 2 + n3 ) n32 − n 22 + 2V442 (n 2 + n3 ) n 23 V222 + V232 ,
2
) Expression for stiffnesses in terms of group velocities
( )] (
− n33 V332 + V232 + 4 n34V332V232 − n 24V222 V232 = 0 ) For completeness, expressions for the nine stiffnesses, for
the case of orthorhombic symmetry, are here summarized.
(A-40) These equations follows directly from (A-21) to A-23), (A-
36), and (A-39), as well as their cyclic variant.
in which all of the Vij have been measured experimentally.
Now, since n 2
2
+ n32 = 1 , equation (A-40) can, in C 11 = ρ V 112 

principle, be solved for n2 and n3. In practice, we determine C 22 = ρ V 222  (A-45)
n2 and n3 by iterative substitution. Knowing n2 and n3 and C 33 = ρ V 332 
v44 from equation (A-39), equation (A-36) and (A-37) can
be solved for C23 and the polarization, p2/p3. Similarly,
C 44 = ρ V 232 = ρ V 322 
using V44 (qSV ) in (A-40), we will get different values 
C 55 = ρ V 312 = ρ V 132  (A-46)
(in general) for n2, n3 and v 44 ; but assuming the C 66 = ρ V 122 = ρ V 212 
orthorhombic model to be a reasonable one, we should get
the same result for C23.
C 23 =
ρ
n 2 n3
{[
1
2
(n2 + n3 )2 V442 − n32V232 − n32V332 ] (A-47a)
The SH wave. – Choosing equation (A-34) instead of (A-
33) we have, from (A-14) and (A-31):
[1
2
(n2 + n3 )2 V442 − n22V222 − n32V232 ]}
12
− ρ V232
ρv 2
= n C 66 + n C 55 = Γ11
2 2
(A-41)
ρ
{[ (n ]
41 2 3

C 31 = + n1 ) V552 − n32V312 − n12V112

1 2
so that the Kelvin-Christoffel equations (A-12) result in n3 n1 2 3 (A-47b)
p k = (1 , 0 , 0 ) . (A-42)
[ (n
1
+ n1 ) V552 − n32V332 − n12V312
2
]}
12
− ρ V312
2 3
Incorporating (A-31) and (A-42) into (A-30) and its cyclic
variants yields: C12 =
ρ
n1n 2
{[ (n
1
2 1 + n 2 ) V662 − n12V122 − n 22V222
2
] . (A-47c)
ρv 41V2 = n2 C 66 ρ v 41V 3 = n 3 C 55 .
[ (n + n ) V ]}
and (A-43) 12
− n22V112 − n32V122 − ρ V122
1 2 2
2 1 2 66
Again applying V2=V3 (for ray at 45 degrees) and equation
(A-39), we obtain: an in (A-47) Vii(i=4,5,6) may be replaced by Vii .
n2 C55 2V31V12
= and V41 = . (A-44)
n3 C 66 (V 2
31 + V122 )
12

12