We demonstrate several aspects of processing 3-C data to obtain anisotropy information from shear-wave splitting. This includes: analysis for unknown or uncertain field orientations of the receivers; a new method of splitting analysis which uses both radial and transverse component data to estimate the S1 direction, which is compared with a transverse only method; and some analysis of the difference between P-S1 and P-S2 imaging versus anisotropy-corrected “radial-prime” (R’) imaging. Results of this splitting analysis to characterize changes in the reservoir of an active in-situ combustion project are provided as well. These results highlight the value of 3C data in thermal recovery projects.


Lately there has been increasing interest in using PS data from multicomponent surveys over heavy oil or oil sands reservoirs for characterizing stress – particularly in the overburden and as it relates to caprock integrity analysis. An example of time-lapse stress analysis from PS data for the overburden is described in some detail in Wikel et al. (2012). The basic principle is that the increases in stress associated with production activity (in that case the Toe-to-Heel-Air- Injection or “THAI” method), give rise to anomalies in the amount of shear-wave splitting as measured by the time delay between S1 and S2 (fast and slow shear) modes. Hence, the accurate measurement of shear-wave splitting attributes, such as S1 azimuth and time delay, is of paramount importance for monitoring thermal recovery methods. Shear-wave splitting measurement is also important elsewhere, such as fracture analysis for shale gas plays.

In this paper we will look at several aspects of this process, including analysis to verify correct orientation of the receivers, comparison of two shear-wave splitting analysis methods and some observations from imaging of the S1 and S2 modes directly. This is in contrast to the alternative combined result from correction and rotation back to radial, sometimes referred to as “radial-prime” (R’).

Different aspects of these processing methods are demonstrated using two heavy oil datasets. One of these is a medium sized survey (13,000 shots) which was shot with relatively sparse receiver line spacing. The second (Kerrobert) is a small survey (<1000 shots) which had tighter acquisition parameters and was shot with a focus on stress analysis for a THAI production field. Results from the Kerrobert survey are then analysed in more detail in the final section of this paper.

Orientation Analysis

It is desirable when acquiring a 3-C survey to consider the information required to properly process the data. From a geometry point of view, this means not only accurate field positioning of source and receivers, but also accurate orientation of the inline direction of the 3-C receivers. Without reliable orientation information, estimation of azimuthally dependent quantities required for splitting analysis is impossible, and even simple radial PS output is adversely affected. A typical practice is to orient the receivers towards magnetic north, and apply the magnetic declination correction in processing. Sometimes practical considerations (e.g. difficult ground conditions) can compromise this ideal. In these cases, it may be necessary to derive an estimate of the receiver orientation directly from the recorded data.

Various methods for orientation analysis have been published (e.g. Li and Yuan, 1999; Olofsson et al., 2007), but these generally are addressing the more challenging (and less well constrained) problem of estimating three Euler or equivalent angles for a geophone of unknown tilt and bearing. Here we assume the tilt is vertical, and only the bearing is in question.

The geometry for orientation analysis is illustrated in Figure 1. Suppose we assume that the orientation of H1 is known to be Ψ, and that H2 is 90° clockwise from H1. Then for shots j=1, …, N with shot-receiver azimuth of θj , we can calculate the corresponding radial and transverse data as

Formula 1a

where H1j(t), H2j(t) are the H1 and H2 component records from shot j, and A is a rotation matrix

Formula 1b
Fig. 01
Figure 1. Geometry of orientations.

First breaks on R are assumed to be P-wave arrivals which should correlate maximally with the vertical component, V. This should apply for all shots into the receiver and we can construct a correlation function:

Formula 2

Where Vj is the vertical phone record for shot j. Ideally, the correlation is only performed over a window which contains first arrival data, to avoid contamination from reflections which may have been influenced by anisotropy. Now, regarding the receiver orientation Ψ as an unknown in equation (2) we search for a value of Ψ which maximizes F(t).

Fig. 02
Figure 2. Comparison of field orientations (a) and estimated orientations (b) using first break radial-vertical match. The absent field orientations are locations where receivers had been removed prior to field measurement. In (c) the estimated errors are shown, based on standard deviation of individual shot-based orientation estimates. Circles in (b) and (c) indicate the location of the gathers shown in Figure 3.

In the Kerrobert survey illustrated in figure 2, and later shown as our case study, the orientations of some of the receivers were measured after the survey was acquired (the receivers with no line segments in the figure had already been picked up). There was enough variance from the nominal orientation (magnetic North) that we decided to estimate orientations of all the receivers from the data, providing an opportunity to compare measured and estimated values where available.

The method described above for orientation analysis is based on the assumption that first breaks measured on both vertical and horizontal components are radially oriented P-wave arrivals. This assumption warrants some scrutiny, as it neglects near surface effects such as anisotropy and scattering. However, we found the analysis to be fairly robust when applied over many shots for each receiver.

We can also perform the maximization of F(t) in equation (2) for each shot separately. The resulting estimates of Ψ have some scatter and the standard deviation of these is used to obtain an error estimate. Because these single shot measurements are quite noisy, the resulting standard deviation is generally larger than our actual errors, but is nevertheless useful to identify suspect locations. Figure 2(c) shows the error measurements for the orientations estimate in figure 2(b). We then visually compared the radial and transverse rotated data using the orientations which were recorded with those estimated by our method, at a number of locations including both high and low error spots. An example is indicated by the circles in figure 2(b) and (c), where the error QC was about 43°. Figure 3 shows the shallow data for this location, for offsets from 400m to 1100m, after flattening on first breaks and bulk shifting by 200ms. The analysis was performed on first arrival energy above 250ms on the shifted data, and over offsets up to 1100m. At this location, the estimated orientation deviates from the recorded orientation by 47°, giving rise to significant differences between the resulting rotations. The rotation to radial and transverse indicates a more plausible distribution of first break energy (higher on radial, lower on transverse) when the estimated orientation is used. We found that in general the orientations estimated from first breaks were at least as plausible as the recorded ones. Ultimately, we processed the entire dataset with the estimated values.

Fig. 03a Fig. 03b
Figure 3. Comparison of data rotated with field orientation and orientation estimated from first break analysis. The gathers have been flattened on first-breaks, and then bulk shifted downwards by 200ms, and are all scaled with the same single scalar. The recorded orientation is approximately 17° clockwise from North, whereas the estimated orientation is 330° clockwise from North.

Shear-wave Splitting Analysis

Shear-wave splitting analysis can be viewed as a two-step process: first, estimate the azimuth of the S1 direction, and; second, estimate the time-delay between S1 and S2. Typically the transverse component is used to estimate the S1 azimuth, based on the azimuthal position of polarity flips or using an amplitude fitting approach (see for example Li, 1998). We refer to this as the T-only method.

However, amplitude information from the radial can also be incorporated in the process. We describe an algorithm to do this, based upon the theoretical form of the radial and transverse amplitude variations for an S1 direction ϕ, and N traces with radial directions θi . Assuming there are two time functions A1(t) and A2(t) which correspond to the S1 and S2 reflected amplitudes, then the radial and transverse amplitudes R(t, θi ) and T(t, θi ) can be written as

Formula 3

for i = 1, ..., N, where A(t) = ½ [A1(t) + A2(t)]

and ΔA(t) = ½ [A1(t) – A2(t)] .

Equation (3) can be written in matrix form as

Formula 4a


Formula 4b

which can be inverted, provided N ≥ 2, to give the usual leastsquares estimate

Formula 5

The second and third terms of the solution vector X can then be combined to determine ϕ. We refer to this as the R+T method.

What is the advantage in using both radial and transverse as opposed to transverse alone? A common assumption is that, since the signal-to-noise ratio on the radial is higher than on the transverse, then an estimate using radial data in addition to transverse will be more robust. However, equation (3) makes it clear that the signal of interest is the deviation ΔA(t) from the background amplitude A(t). This deviation signal is no stronger (on average) on radial than it is on transverse. However, it does have a different sensitivity to the azimuth distribution, and so may help in areas where the azimuth distribution is poor.

Formula 4
Figure 4. Comparison (from left to right) of S1 azimuth estimates with superbinned input, using 5x5 T-only, 5x5 R+T, 11x11 T-only and 11x11 R+T.

Figure 4 shows a comparison of S1 azimuth estimates measured on a heavy oil dataset at the Devonian level (about 1 second PS time) using 5x5 and 11x11 superbin sizes between T-only and R+T based methods. The receiver line spacing for this survey is 150m, which causes a noticeable acquisition footprint that contaminates the analysis results for shallow targets. As is seen in figure 4, for 5x5 superbin sizes, the footprint is somewhat reduced using the R+T method compared to the T-only method. However, if the superbin size is increased to 11x11, which is necessary to remove the footprint satisfactorily, then both methods perform equally well.

For the Kerrobert data, which had tighter acquisition parameters so that footprint was not a significant problem, we utilized analysis based on transverse component only.

P-S1 and P-S2 versus R’

There are two alternative ways to image the PS data after shearwave splitting analysis.

Perhaps the most obvious is to rotate data to S1 and S2 directions and image these two datasets (P-S1 and P-S2) separately. This is not as straightforward as it first sounds, as both P-S1 and P-S2 have azimuthal amplitude variations related to the sourcereceiver azimuth, including polarity changes. These must be dealt with appropriately, for example by the weighted stacking method of Bale et al. (2000). In this approach the P-S2 can be left in delayed time, or shifted to match P-S1.

Fig. 5a Fig. 5b Fig. 5c
Figure 5. Illustration of process from original Radial and Transverse through splitting correction to derive Radial and Transverse “prime”. The input data (a) is shown as common asymptotic conversion point (ACP ) gathers where data are indexed first by offset range, and then by shot-receiver azimuth, (see profiles above the gathers). The gathers are superbinned using 7x7 bins to improve S/N. The blue dashed box shows the range used for shear-wave splitting analysis. Using S1 azimuth derived from the analysis, the data are rotated to the coordinate system defined by S1 and S2 direction (b). Using a time delay estimated from the analysis (about 10ms), the S2 data is shifted upwards (right panel of b) to align with the S1 data. Finally in (c), the data are rotated back to the R-T coordinate system where the removal of the anisotropy effect is evident. These are referred to as R’ and T’.

The alternative approach is to shift the P-S2 data to match P-S1 first and rotate back to the radial-transverse coordinate system. The resulting data may be referred to as radial-prime and transverse- prime (R’ and T’). This procedure is illustrated in figure 5, using the Kerrobert dataset shown later as our case study.

In order to layer strip through multiple layers, it is necessary to compute R’ and T’ for each layer in succession, in order to perform analysis for the subsequent layer. Often the R’ dataset is the most suitable for the final goal of imaging the reservoir.

Fig. 6a Fig. 6b Fig. 6c Fig. 6d
Figure 6 (right). Comparison of amplitudes for different PS imaging methods. The three stacks in (a) show P-S1, PS2 and R’ results, with the S1-S2 time delay profiled above them. The amplitude spectra in (b) are for the window from 700-1350ms. Note that P-S1 is richer in high frequencies and P-S2 is richer in low frequencies. However, at a different location (c) the reverse behaviour is observed, as seen in (d).

Nevertheless, there can be advantages in analysis of P-S1 and PS2 images separately, in particular for amplitude subtleties or spectral response. Figure 6 shows a comparison of data from the same heavy oil survey used in figure 4, processed for P-S1 and P-S2 images, and compared with a R’ image. Figure 6(a) shows a portion of the image from one line, while figure 6(b) compares the amplitude spectra for the three images. What is notable is that the P-S1 image is richer in the higher frequencies while the P-S2 is biased towards lower frequencies. This suggests (assuming equivalent reflectivity responses) that there is some degree of differential attenuation being observed between the two modes, with the P-S2 mode having greater frequency dependent losses than the PS1 mode. We have observed this to be fairly prevalent within this dataset. However, there are examples elsewhere of the reverse behavior, as shown in figure 6(c). In both cases, the R’ spectrum benefits from the combined contributions of P-S1 and P-S2.

Case Study: Splitting Analysis for Recovery Monitoring over a THAI Reservoir

The Kerrobert field, near the town of the same name in Saskatchewan, is the location of the first commercial THAI operation in Canada (two THAI pilots have been operated in Canada since 2005). The reservoir of interest is ~750m deep and consists of well consolidated sands that range in thickness from 15-30m. The caprock in this area consists of over 50m of marine shales in the Mannville group, providing thick continuous containment of the process. The reservoir has been developed using cold production methods, however, as in most heavy oil reservoirs, primary recovery is quite low at 1.5% of original oil in place. The pilot phase for THAI utilized two wells (KP/KA 1-2) from 2009-2011 with another 10 wells drilled for commercial development in early 2011 (KP/KA 3-12). The Kerrobert expansion wells were brought online in July-August 2011, while the first 3C seismic shoot occurred in October 2011.

The THAI process utilizes one air injector/horizontal producer pair for insitu combustion. A detailed overview of the THAI process can be found in Kendall (2009) and Wikel et al. (2012). We know from reservoir modeling and previous piloting experience that the THAI front and its advancement change the pore pressure and fluid content within the reservoir as combustion progresses. It follows that a change in fluid state, gas content, and an increase in pore pressure with increased air injection should result in a change in seismic response over time and should produce 3C anomalies at the time of the shoot. It is believed that stress is the dominant factor in shear wave splitting anomalies within the reservoir, however, small contributions in splitting from fluid phase changes cannot be ruled out.

4D-3C seismic was chosen as the primary method for monitoring these variations through time, as the changes throughout the reservoir and caprock will be heterogeneous, meaning that simple wellbore monitoring alone is insufficient. The examples below are the baseline for future 3C shoots.

Shear-wave splitting analysis at Kerrobert was performed using four layers, selected using a combination of stress analysis objectives and signal quality considerations. These layers correspond to overburden, two layers above the caprock (the second layer bounded at the bottom by the caprock) and a final layer over the reservoir. Figure 7 shows splitting analysis results from the overburden and the layer bounded by the caprock. Figures 8 and 9 are focused on the reservoir analysis.

Fig. 7
Figure 7. Overburden (a) and caprock (b) analysis for shear-wave splitting. The colour underlay represents the time delay between S1 and S2 modes, with a range from 0 to 15ms, as shown in the histograms. The line segments represent both time delay (length) and S1 azimuth.

As has been observed numerous times in heavy-oil PS data (Whale et al., 2009; and Cary et al., 2010) the most significant Swave anisotropy effect, at least in a raw time-delay sense, comes from the near surface, as seen in figure 7a. The measured overburden anisotropy has a mode around 10ms which, for analysis over 750ms from surface, corresponds to an average S-wave anisotropy of about 1.3%. Though this is not high in percentage terms, it can have a noticeable impact on imaging if not corrected for. Of course, layer stripping of this layer is also critical for obtaining meaningful analysis results for subsequent layers. The area between overburden and caprock was divided into two layers. Results from the second of these, after correction for the first, are shown in figure 7b. This might be characterized as largely isotropic with isolated pockets of measureable splitting.

We now turn our attention to the final layer, which was analyzed after layer stripping overburden and caprock layers.

Fig. 8
Figure 8. PS1/PS2 time lag in colour with the PS1 direction overlain as a vector scaled by the time lag in October 2011. White wells are THAI development air (directional vertical wells) and production (horizontal wells). KP1 and KP2 were the initial THAI pilot wells that have been injecting air since November 2009. KP3-12 were brought online with air injection in July- August 2011 and are in the startup stages of THAI. Red wells are legacy primary wells.

Figure 8 shows the PS1/PS2 time lag with the PS1 direction referenced to true north overlaying the horizon as a vector (vector size scaled by time lag). In this display we normalize the time lag (~7ms) by the time interval from the previous layer to this one (75ms), to get the percent anisotropy over the reservoir. Figure 9 shows a zoom of anomaly in figure 8, as raw time lags (not normalized), with a different colour scale to highlight the PS1 vectors. As shown in figure 8 and 9 there is a large anomaly which shows 10% anisotropy over the window. This is a large anomaly, given that layer stripping was performed on three layers overlying this reservoir. In addition, this large PS time lag anomaly is located at the bottomhole location of the two air injectors in the field that have injected the largest volumes of air and have been implementing THAI for the longest period in the field (KA/KP 1-2). These two air injectors and horizontal producers were the location of the THAI pilot in this field, as noted.

Fig. 9
Figure 9. Zoom of the anomaly in figure 8 with a different colour scale to highlight the PS1 vector attribute (wells are not shown).


PS data is proving useful for identification of stress variation in the near surface over heavy oil reservoirs, in the absence of natural fracturing. To enable this analysis, it is important to accurately measure receiver orientation, or, if necessary, to estimate it as well as possible from the data.

The first step in splitting analysis is to estimate S1 azimuth, which is typically done using the transverse component alone (T-only). We found some marginal benefit to a method using both radial and transverse (R+T), which reduced the footprint on azimuth estimates with large receiver line intervals, but not enough to obviate the need for substantial superbinning.

A comparison of P-S1 and P-S2 data after weighted stacking indicates possible differential attenuation between these modes and that the R’ spectrum benefits from the combination of the two modes.

The results of this type of analysis, when processed correctly, are crucial for monitoring reservoir changes that may not be apparent on PP data only. This includes monitoring overburden and reservoir changes via converted wave splitting. Monitoring in this fashion provides areal coverage that cannot be achieved with wellbore monitoring alone. 3C methods are well suited for most thermal and EOR processes where appreciable reservoir changes occur, such as CO2 floods, cyclic steam stimulation (CSS), and steam assisted gravity drainage (SAGD).



We thank Key Seismic Solutions and Petrobank Energy and Resources for permission to publish this work. We also would like to thank Jeff Deere for his careful review.

About the Author(s)

Richard Bale is currently Director of R&D for Key Seismic Solutions. He is a graduate of the University of Cambridge in the UK, obtaining both a BA and an MA in mathematics. He began his geophysical career with GECO in the UK in 1984, and continued with Schlumberger Geco-Prakla, in both the UK and Norway, as a senior research geophysicist. From 1997 to 2001 he was the project manager for multicomponent processing development for Geco-Prakla/WesternGeco, responsible for technology applied on 4-C data in the North Sea and Gulf of Mexico. Richard moved to Calgary in 2001 to pursue research with the CREWES consortium and obtained a Ph.D. in geophysics from the University of Calgary in 2006. Richard worked for Veritas/CGGVeritas in Calgary from 2004 to 2011, where he was most recently Manager for Multicomponent Research and Technology, before joining Key Seismic in 2011. He has authored or co-authored a number of SEG, EAGE and CSEG papers on aspects of multicomponent processing, including converted-wave binning, shear-wave splitting analysis and elastic migration.

Tobin Marchand graduated with a B.Sc. Joint Honours in Mathematics and Physics in 2002 and a M.Sc. in Applied Mathematics in 2005, both from McGill University. Tobin started his career as a junior seismic processor with CGGVeritas; then Veritas DGC. Over the span of 5.5 years with CGGVeritas, Tobin’s exposure widened from processing 2D data to 3D data, then onto 3D multicomponent processing to finally becoming an instrumental member of the multicomponent production support team. In October 2011, Tobin joined Key Seismic Solutions Ltd. as the multicomponent team leader and principle multicomponent processor.

Keith Wilkinson graduated in 1983 from the University of Saskatchewan with a B.Eng. degree in Geological Engineering – Geophysics Option. He spent 12 years with Shell Canada Limited primarily in seismic acquisition and processing roles. Keith moved into the service sector in 1996 and held technical management positions at Enertec, Veritas ,and GSI. Keith was a founding partner of Key Seismic Solutions Ltd. in 2004 and is currently Vice President – Technology with Key Seismic.

Kurtis Wikel graduated with a B.Sc. in Geology from the University of Montana. He received his M.Sc. from the Colorado School of Mines in Geophysics with a minor in Petroleum Engineering in 2008. A graduate of the Reservoir Characterization Project, he worked with Schlumberger DCS Denver on geomechanics applied to time lapse seismic data. He worked for Talisman Energy for 2.5 years in exploration geophysics and was the International Exploration Geomechanics Specialist until July 2010. Kurt is currently working on Subsurface Geophysics and Geomechanics for Petrobank Energy and Resources in Calgary, Alberta, Canada.

Robert Kendall is currently President of Kendall Geophysics Ltd., specializing in full-cycle multicomponent and timelapse seismic analysis. From 2006-2012, he was the Subsurface Manager at Petrobank Energy and Resources Ltd where he managed the geology, geophysics and reservoir engineering departments. From 1999 until 2006 he worked for Veritas GeoServices where he started the Multicomponent Processing Group. In 1992 he began work for Amoco Exploration and Production in Houston, he then moved to the Amoco Tulsa Research Center and eventually the Offshore Business Unit in New Orleans. While at Amoco he worked with shear data (both pure- and converted-mode) from the western United States, the Netherlands, Norway, the Gulf of Mexico (US and Mexico) and many of the fields of western Canada in both the borehole and surface seismic environments. From 1988 until 1990 he worked in seismic data acquisition with Solid State Exploration mostly in western Canada and northwestern U.S.A. He received his M.S. from the Colorado School of Mines in 1992 and his B.Sc. from the University of Calgary in 1988.


Bale, R., G. Dumitru, and T. Probert, 2000, Analysis and stacking of 3D converted wave data in the presence of azimuthal anisotropy: 70th Annual International Meeting, SEG, Expanded Abstracts, 1189–1192.

Cary, P., X. Li, G. Popov, and C. Zhang, 2010, Shear wave splitting in compliant rocks: The Leading Edge, 29, 1278–1285.

Kendall, R. [2009] Using time lapse seismic to monitor the THAI heavy oil production process. SEG Expanded Abstracts, 2009, 3954-3958.

Li, X. -Y., 1998, Fracture detection using P-P and P-S waves in multicomponent sea-floor data: 68th Annual International Meeting, SEG, Expanded Abstracts, 2056–2059.

Li, X. Y., and J. Yuan, 1999, Geophone orientation and coupling in three-component sea-floor data: A case study: Geophysical Prospecting, 47, 995–1014.

Olofsson, B., C. Massacand, and M. Millington, 2007, Determining 3-C geophone orientation from a single seismic shot: SEG Expanded Abstracts 26, 970-974.

Whale, R., R. Bale, K. Poplavskii, K. Douglas, X. Li, and C. Slind, 2009, Estimating and compensating for anisotropy observed in PS data for a heavy oil reservoir: 79th Annual International Meeting, SEG, Expanded Abstracts, 28, 1212–1216.

Wikel, K., R. Kendall, R. Bale, J. Grossman, K. DeMeersman, 2012, 4D-3C Geomechanical Study of In-Situ Bitumen Recovery Utilizing Toe to Heel Air Injection THAI®), Alberta, Canada: First Break, 30, 43-53.


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