Elastic anisotropy is a seismic attribute defined as the variation of wave speed as a function of its propagation direction within a medium. With the increase in interest of unconventional reservoirs, the incorporation of anisotropic factors is a necessity to provide more accurate seismological models. Without consideration of anisotropic effects, sub-surface imaging (especially from wide-azimuth and far-offset data) can produce erroneous anomalies and lead to misinterpretation (Wild, 2011). Although reflection seismic data can provide insight into the anisotropic velocity fields (Alkhalifah & Tsvankin, 1995), laboratory studies at the core scale are ideal for characterizing the magnitude of anisotropy as it allows for greater control and precision of the testing conditions.

These measurements are important to study the anisotropy of the rock and the relationship to its petrophysical properties where they may be observed at the seismic scale. In this study, ultrasonic velocity measurements are conducted on Duvernay shale samples to gain insight into that formation’s elastic properties.


Shales tend to make up a substantial component of the average sedimentary basin by volume (Johnston & Christensen, 1995; Jones & Wang, 1981) and are presently utilized as reservoir rocks. They also tend to exhibit a relatively large degree of anisotropy as the orientation of microcracks, preferential clay platelet alignment, and layering all contribute significantly to this phenomenon. Such features can be seen in Figure 1.

Fig. 01
Figure 1. SEM of sample ECA HZ 5-13 showing alignment of clay minerals and crack orientation. Way up direction to the right.

We model the overall elastic compliance of the rock as

Equation 01

where Equation is the compliance of an anisotropic medium void of compliance cracks and Equation is the excess compliance due to the clay platelet contact regions (Sayers & Kachanov, 1995; Mavko, et al., 2009). In addition, Sayers (1999) showed that these cracks are more compliant in shear than in compressional stress. Therefore, at large values of hydrostatic stress where these contact regions of the clay platelets close, we assume the rock to have an elastic compliance of Equation and inversely, an elastic stiffness Equation. Schijns, et al. (2012), for example, have recently used this theory to show that cracks must be considered in the overall anisotropy of a thick formation from comparisons of laboratory and walk-a-way vertical seismic profile measurements.

Here, we make the assumption that the samples are homogeneously transversely isotropic with a vertical axis of symmetry (VTI). For this case, the elastic stiffness tensor is of the form (King, 1970)

Equation 02

where C is the Voigt elastic stiffness matrix with 5 independent elastic constants. Equivalently, Equation 2 can be expressed as the elastic compliance tensor with elastic moduli more familiar to the engineering community (Meléndez-Martínez, 2014)

Equation 03

where E is the Young’s modulus, ν is Poisson’s ratio, and μ is the shear where the subscripts denote the 3 orthogonal axes X1, X2 (perpendicular to symmetry axis) and X3(parallel to symmetry axis). More specifically, the Young’s modulus E1= E2 is the ratio of the uniaxial stress applied to the isotropy plane to its strain along the same axis. E3 would be the same ratio, but with stress applied and strain measured perpendicular to the isotropy plane. The shear modulus μ23 is defined as the ratio of shear stress to shear strain in the X2 direction on the plane with a normal in the X3 direction. Poisson’s ratio is defined as the negative ratio of transverse to axial strain where ν31 corresponds to the stress applied perpendicular to the isotropy plane ( X3 ) with strain measured parallel (X1 = X2) . Similarly, ν12 corresponds to stress applied parallel to the isotropy plane with strain measured parallel as well (Sayers, 2013).

From the elastic constants, Thomsen’s parameters of anisotropy (Thomsen, 1986) can be calculated for a VTI medium

Equations 04-06

where and are dimensionless parameters used to qualify the P-wave and S-wave anisotropy respectively. The parameter can be viewed as the anellipticity of the P-wavefront.

For estimations of the elastic constants needed to fully describe transversely isotropic media, Daley & Hron (1977) derived expressions of the 3 phase velocities and as a function of the phase angle and from these expressions, we have

Equations 07-11

where ρ is the bulk density of the material, and the subscripts P, SH, and SV refer to the P-wave, horizontally polarized S-wave, and vertically polarized S-wave respectively. These equations allow us to estimate the elastic constants through the ultrasonic wavespeeds with the variation of the propagation direction. Examination of Equations 7-10 show that most of the elastic moduli can be determined relatively accurately by simple measurements of wave speeds along principle axes of the material. Determination of C13, however, is not easily achieved in practice as propagation of errors through Equation 11 results in large uncertainties.

An alternative equation for C13 can be derived by revisiting the phase velocity equations from Daley & Hron (1977). Hemsing (2007) proposed Equation 12 to greatly reduce the uncertainty by removing the influence of several elastic constants. The expression, instead, makes use of the 45° vertically polarized shear velocity in addition to the 45° P-wave velocity.

Equation 12


In this contribution, the degree of anisotropy is investigated experimentally using the Pulse Transmission Method using piezoelectric ceramic transducers (Bakhorji, 2010; Hemsing, 2007; Kebaili, 1997). These transducers produce an impulse with the application of an electric potential, which when mounted on the surface of a sample will propagate through as an elastic wave. Upon reaching the other side, another transducer mounted to the surface converts this mechanical disturbance back to a voltage due to its piezoelectric properties. Taking the ratio transducer offset to the travel time of the signal through the sample gives us an estimation of the phase velocity at ultrasonic frequencies.

These elastic constants can be estimated from bulk density of the material and the P-wave and S-wave phase velocities measured perpendicular (VP (0°)), parallel (VP (90°), VSH(90°)), and 45° parallel (VP(45°), VSV(45°)) to the foliation plane. The most common method of measuring these velocities is through the production of 3 cylindrical cores obtained from the main sample (Figure 2) as carried out for example by Meléndez-Martínez & Schmitt (2013). For this study, we will instead adapt a single prismatic sample configuration after Wong, et al. (2008) and Meléndez-Martínez (2014).

Fig. 02
Figure 2. Common coring configuration for anisotropy measurements. Meléndez-Martínez (2014)

The sample configuration consists of cutting the material into a prism-like geometry (Figure 3) as opposed to the traditional 3 cylindrical cores, where each opposing face is cut to provide the propagation directions needed to estimate the elastic constants of a VTI medium (0°, 45°, and 90° to bedding). This configuration allows us to overcome the heterogeneity issues of measuring the velocities through multiple cores and reduces issues associated with attempting to core weak materials in set directions.

Fig. 03
Figure 3. Model of the multi-face core and proposed ultrasonic pulse orientations from Meléndez-Martínez (2014).

The setup consists of attaching copper pads directly to the sample to act as common grounds for the transducers. The transducers themselves are then glued onto the copper pads using a silver conductive epoxy with another copper pad on top acting as the positive ground. For static measurements, linear foil strain gauges are mounted parallel and perpendicular to sample bedding to measure deformation. Wires are soldered onto each electrode and dried before the application of a non-conductive urethane compound around the sample. The compound is used as a membrane to protect the transducers and sample material from the confining fluid of the pressure vessel as the samples are measured under dry and undrained conditions. Figure 4 shows a picture of the completed sample before the application of the compound. The entire sample is then placed into the pressure vessel with wires fed out through the top and connected to a data acquisition setup consisting of an oscilloscope, power supply, pulser / receiver, and digitizer.

Bulk density values of the samples were obtained from Mercury Intrusion Porosimetry.

Fig. 04
Figure 4. Picture of completed sample before application of urethane compound


In this section, the phase velocities as well as estimations of the static and dynamic moduli of one sample are shown as results of the experiment. Figure 5 shows the P-wave waveforms for the propagation direction perpendicular to bedding. Changes in arrival time with increasing pressure can be observed directly from this figure.

Fig. 05
Figure 5. Normalized P-wave waveforms for the propagation direction perpendicular to bedding.

Travel times are picked from the first extremum of each waveform and are marked by arrows in Figure 6. These times are then shifted to the first break through a calibration to achieve proper timing. For additional details on this calibration, refer to Meléndez-Martínez (2014).

Fig. 06
Figure 6. P-waveform at a confining pressure of 50 MPa. Arrows indicate point at which travel time is picked.

In addition to the ultrasonic velocities, deformation is measured parallel and perpendicular to the bedding plane (Figure 7). These strain values will be used to calculate the static elastic moduli for comparison with their dynamic equivalent.

Fig. 07
Figure 7. Deformation of sample measured by linear foil strain gauges in X1 (90° from symmetry axis) and X3 (0° from symmetry axis) directions


Figure 8 plots the phase velocities as a function of the confining pressure with arrows showing the progression of measurements in the pressurization and depressurization cycles.

Fig. 08
Figure 8. Ultrasonic phase velocities plotted as a function of confining pressure. Arrows show progression of measurements.

The observations in the velocities are consistent with what would be expected for a medium with vertical transverse isotropy where VP (0°) < VP (45°) < VP (90°), and VSH (0°) < VSV (45°) < VSH (90°). All velocities increase with confining pressure and the sample becomes stiffer due to the closure of microcracks and fractures in the rock. As a result, the anisotropy also decreases with confining pressure (Table 1), though at high pressures, the large discrepancies between the velocities show that these microcracks and fractures only account for a small source of anisotropy.

In addition to these trends, a hysteresis effect is observed in all measurements due to the microcracks and fractures closing at a faster rate during pressurization than opening during depressurization (Hemsing, 2007).

Table 01
Table 1. Showing estimated Thomsen parameters

Table 1 shows the Thomsen (1986) parameters of anisotropy calculated from the estimated elastic constants where ε and γ are the P-wave and S-wave anisotropy parameters respectively.

From the elastic constants, dynamic Young’s modulus and Poisson’s ratio are calculated with both single crystal and polycrystalline aggregate assumptions (Figure 9).

Fig. 09
Figure 9. Dynamic Young’s modulus and Poisson’s ratio plotted as a function of confining pressure. Arrows show progression of measurements.

The numerical subscripts denote the orthogonal directions X1, X2, and X3 of variables calculated with the single anisotropic crystal assumption. The subscript VRH refers to the Voigt-Reuss-Hill arithmetic average of the upper Voigt and lower Reuss bounds (Hill, 1952) to provide an estimation of the effective isotropic elastic moduli when considering a polycrystalline aggregate. These values are shown as a reference for the moduli calculated with single crystal assumptions.

As expected for shales, the discrepancy in values of Young’s modulus shows that the sample is stiffer in the X1 direction (parallel to bedding) than in X3 (perpendicular to bedding). However, Sayers (2013) examined the relationship between ν31 and ν12 and found that their ratio can be greater than, equal to, or less than 1. An inequality of ν31 > ν12 is expected for a partial alignment of clay particles where as ν31 < ν12 may indicate the presence of kerogen inclusions, microcracks, or low aspect ratio pores parallel to the bedding plane.

Dynamic and static tangent bulk modulus values were also calculated and plotted as a function of confining pressure in Figure 10. As expected, the static measurements are lower than the dynamic measurements where the physical causes can be attributed to the nonlinearity of the stress-strain relation where the ratio differs between measurement over large and small strains. These discrepancies can also be attributed to the sample’s inelasticity or the presence of microcracks (Mavko, et al., 2009; Cheng & Johnston, 1981).

Fig. 10
Figure 10. Dynamic and static tangent bulk moduli as a function of confining pressure. Arrows show progression of measurements.


Ultrasonic velocities experiments were conducted on Duvernay shale samples to investigate their elastic properties with the assumption that the material is transversely isotropic with a vertical axis of symmetry; velocities increase with confining pressure and are consistent with this assumption. Calculated Thomsen’s parameters of anisotropy showed the P-wave anisotropy ε decreased from 0.60 at low pressures to 0.33 at high pressures. Similarly, S-wave anisotropy γ decreased from 0.41 to 0.32.

Plots of the dynamic Young’s modulus and Poisson’s ratio showed that with increasing confining pressure, the closure of microcracks and low-aspect ratio pores result in the sample becoming stiffer. This effect may be indicated by the ratio ν31 / ν12 decreasing from greater than 1 to less than 1.



We would like the University of Alberta, Institute for Geophysical Research, and the Natural Sciences and Engineering Research Council of Canada (NSERC) for their generous support. We would also like to thank the Alberta Geological Survey for providing the samples needed for this study.


About the Author(s)

Oliver N. Ong obtained his B.Sc. in Geophysics from the University of Alberta in 2014 and has since worked as a Research Assistant in the Experimental Geophysics Group at the University of Alberta. Starting this fall, he will be enter a thesis based M.Sc. program in Rock Physics. His research interests have focused on seismic anisotropy and on the pressure and temperature sensitivity of seismic wave speeds in bitumen saturated carbonates.

Jaime Meléndez-Martínez recently completed his Ph.D. in Geophysics from the University of Alberta and he is now a Researcher at the Mexican Institute of Petroleum in Mexico City. His research is focused on studying physical properties of shales, particularly anisotropy. He graduated from UNAM, National Autonomous University of Mexico (BSc in Physics, 2003 and M.Sc. in Geophysics, 2006). During his M.Sc. program he was awarded a Mexico's National Science and Technology Council Scholarship by CONACYT. Prior to starting his PhD program, he worked at IMP, Mexican Institute of Petroleum, developing illumination maps to help understand the influence of the acquisition geometry on the final quality of processed seismic images.

Douglas R. Schmitt is the Canada Research Chair in Rock Physics at the University of Alberta. His team, the Experimental Geophysics Group, is currently carrying out a variety of laboratory rock physics tests to measure fluid properties, rock anisotropy and strength, wave propagation and electrical properties of a variety of rocks under conditions of pressure and temperature. This laboratory work is balanced by field studies particularly in scientific drilling with current projects in New Zealand, Sweden, and India.

Randy Kofman is currently a Research Professional in the Experimental Geophysics Group at they University of Alberta. He holds a B.Sc. in Physics and a B.Sc. and M.Sc. in Geology from the University of Alberta. He collaborates with students and visitors to EGG and is involved with most of the ongoing laboratory and field studies.


Alkhalifah, T. & Tsvankin, I., 1995. Velocity analysis for transversely isotropic media. Geophysics, 60(5), pp. 1550 - 1566.

Bakhorji, A. M., 2010. Laboratory Measurements of Static and Dynamic Elastic Properties in Carbonate, Edmonton, AB: University of Alberta.

Cheng, C. H. & Johnston, D. H., 1981. Dynamic and static moduli. Geophysical Research Letters, 8(1), pp. 39-42.

Daley, P. F. & Hron, F., 1977. Reflection and transmission coefficients for transversely isotropic media. Bulletin of the Seismological Society of America, 67(3), pp. 661-675.

Hemsing, D. B., 2007. Laboratory determination of seismic anisotropy in sedimentary rock from the Western Canadian Sedimentary Basin, Edmonton, AB: University of Alberta.

Hill, R., 1952. Elastic behaviour of a crystalline aggregate. Physical Society of London. Proceedings, 65(6), pp. 349-354.

Johnston, J. E. & Christensen, N. I., 1995. Seismic anisotropy of shales. Journal of Geophysical Research, 100(B4), pp. 5991-6003.

Jones, L. E. A. & Wang, H. F., 1981. Ultrasonic velocities in Cretaceous shales from the Williston basin. Geophysics, Volume 46, pp. 288-297

Kebaili, A., 1997. Velocity anisotropy determination in tau-p space, Edmonton, AB: Geophysics.

King, M., 1970. Static and dynamic elastic moduli of rocks under pressure. Berkeley, CA, Soc. Mining Engineers.

Mavko, G., Mukerji, T. & Dvorkin, J., 2009. Rock Physics Handbook: Tool for Seismic Analysis in Porous Media. 2nd ed. ed. Cambridge, New York: Cambridge University Press.

Meléndez-Martínez, J., 2014. Elastic Properties of Sedimentary Rocks, Edmonton, AB: University of Alberta.

Meléndez-Martínez, J. & Schmitt, D. R., 2013. Anisotropic elastic moduli of carbonates and evaporites from the Weyburn-Midale reservoir and seal rocks. Geophysical Prospecting, Volume 61, pp. 363-379.

Sayers, C. M., 1999. Stress-dependent seismic anisotropy of shales. Geophysics, 64(1), pp. 93-98.

Sayers, C. M., 2013. The effect of anisotropy on the Young’s moduli and Poisson’s ratios of shales. Geophysical Prospecting, 61(2), pp. 416 - 426.

Sayers, C. M. & Kachanov, M., 1995. Microcrack-induced elastic wave anisotropy of brittle rocks. Journal of Geophysical Research, 100(B3), pp. 4149-4156.

Schijns, H., Schmitt, D. R., Heikkinen, P. J. & Kukkonen, I. T., 2012. Seismic anisotropy in the crystalline upper crust: observations and modelling from the Outokumpu scientific borehole, Finland. Geophysical Journal International, 189(1), pp. 541-553.

Thomsen, L., 1986. Weak elastic anisotropy. Geophysics, 51(10), pp. 1954-1966.

Wild, P., 2011. Practical applications of seismic anisotropy. First Break, Volume 29, pp. 117 - 124.

Wong, R. C. K., Schmitt, D. R., Collis, D. & Gautam, R., 2008 . Inherent transversely isotropic elastic parameters of over-consolidated shale measured by ultrasonic waves and their comparison with static and acoustic in situ log measurements. Journal of Geophysics and Engineering, 5(1), pp. 103-117.


Join the Conversation

Interested in starting, or contributing to a conversation about an article or issue of the RECORDER? Join our CSEG LinkedIn Group.

Share This Article