A 3D seismic data set was acquired near Alder Flats, Alberta, to image several coal zones of the Ardley Coals. Deterministic inversion of the data shows a band-limited acoustic impedance result that over estimates the very low impedance of the coals as well as overestimating the thickness of the coal zones. Additionally, the thin sub-zones of the Ardley Coals are not differentiated in the inversion. A geostatistical inversion is conducted in an attempt to improve the inversion resolution and assess model uncertainty. The inversion provides multiple model realizations each of which honors the seismic data, the well data and the geostatistics. The mean of the realizations gives a highly resolved estimate of the acoustic impedance. As well, the multiple model realizations provide a range of acoustic impedance and net coal-thickness estimates that can be used to quantify uncertainty in the results. However, while the inversion appears to be accurate at the level of the shallow-most coals, at the deeper coal subzones a thin but significant sub-zone is not detected in the inversion, even when the inversion is constrained by the well data. This phenomena is possibly explained by the results of synthetic modelling that investigates the effects of short-path multiple reflections. The modelling shows that the effect of short-path multiple reflections, built up in the seismic pulse from the stacked coal sub-zones, could cause a phase delay in the reflected energy of the deepest coals and which would explain the results of the inversion.


A 3D vertical component seismic survey was conducted at the Alder Flats location in June 2007. Details of the acquisition design, processing, deterministic inversion via several methods and final interpretation are discussed in McCrank and Lawton (2009) and McCrank (2008). This paper investigates a geostatistical inversion of the data using the Fugro Jason – Jason Geoscience Workbench (JGW) StatMod algorithm. Discussion of the method can be found in Torres-Verdin et al. (1999) and van der Laan and Pendrel (2001).

Fig. 01
Figure 1. The field layout of the 3D survey.

The 3D field layout (Figure 1) was a square of approximately 560 x 560 m. The survey was intended to image the stratigraphy of the Ardley Coal Zone which is illustrated in the well correlation cross-section shown in Figure 2. The Ardley Coal Zone is informally broken into the Upper and Lower Ardley coal zones. The Upper Ardley is an approximately 10 m gross thickness coal zone and the Lower Ardley coal zone is represented by two thinner zones, the Silkstone coals and the deeper Mynheer coals that have gross coal thicknesses of approximately 4 m and 8 m, respectively.

Fig. 02
Figure 2. An approximately 1 km long north-south geological cross-section of the Ardley Coal Zone traversing the 3D seismic volume.

Well #2 (Figure 1 and Figure 2) is inside the survey bounds and has high quality sonic and density logs that were used for wavelet extraction and petrophysical analysis (discussed below). Well #3, although within the survey bounds, had a poor quality density log and no sonic log that was not used in the analysis.

The 3D volume of seismic data (Figure 3a) contained a wavelet (Figure 3b) that was significantly tuned in the thin but highly reflective beds of the various coal zones.

Fig. 03
Figure 3. (left) the 3D seismic volume and (right) the seismic wavelet.

The results from deterministic Zp inversion

After carefully extracting a wavelet, the post-stack seismic data were deterministically inverted to estimate acoustic impedance. Figure 4 shows a cross-section of the inverted Zp volume. This west-east section that intersects the Well #2 shown in Figure 4 will be the section shown in all subsequent cross-section figures in this paper.

Figure 5 shows a cross-plot of Zp versus gamma ray response in the well logs from Well #2. Coals can be identified with an upper limit Zp of 6x106 kg/m3*m/s and the lithologies with higher Zp are defined as "clastics" (representing shales, siltstones and sandstones). The data for the cross-plot are generated at "logging" resolution as recorded in the field (resolution of approximately 15 cm to 1 m).

Fig. 04
Figure 4. A cross-section (EW) through the acoustic impedance volume that intersects Well #2. The acoustic impedance well log is superimposed on the section where the well log data have been high cut filtered at 60-70 Hz in order to match the resolution of the band-limited inverted data.

It is known that the Zp estimate from seismic inversion will have lower band-width and less resolution than the well logs, so the cross-plot is revisited in Figure 6 where the well log data have been logged in time (1 ms sampling) and high-cut filtered at 60- 70 Hz. The lithology cut-offs are unchanged, but the resolution is diminished. The net effect can be seen in the Zp logs and the lithology logs of Figure 7. The Zp log shows that the vertical detail is reduced after high-cut filtering and that the extreme amplitudes of the low Zp coals are smoothed and increased when high cut filtered. The lithology logs resulting from these two Zp curves are also shown in Figure 7. It can be seen that the thin coal sub-zones in each of the Lower Ardley coal zones are "smeared" into one coal zone when the lithology is established based on the high-cut filtered Zp log. Also, the net thicknesses of the zones in the high-cut filtered lithology log are incorrect.

Fig. 05
Figure 5. Cross-plot of the acoustic impedance versus gamma ray. The coal facies is defined with acoustic impedance less than 6x106 kg/m3*m/s.
Fig. 06
Figure 6. Cross-plot of the acoustic impedance versus gamma ray after the well logs have been high-cut filtered at 60 Hz.
Fig. 07
Figure 7. (a) The Zp log at logging resolution and after 60 Hz high-cut filtering, (b) the lithology log generated from the Zp log at logging resolution and (c) the lithology log based on the high-cut filtered Zp logs.

The need for geostatistcal inversion

Deterministic seismic inversions yield a result that has the maximum likelihood of being correct given the seismic data alone. However, because the Zp estimate inverted from the seismic data is band-limited in approximately the same way that the well log data just discussed, the results will suffer from the same shortcomings. Specifically:

1. The absolute values of the Zp estimate will be higher than the actual very low Zp values of coal facies.

2. The thickness of the coal zones will not be accurate because at the frequency resolution of the seismic data, the low impedance of several coal zones will be smeared together and result in an "average" low Zp inversion zone with an overestimated thickness.

These limitations are illustrated in Figure 8 that shows the Zp log of Well #2 at "logging" resolution superimposed on the deterministic Zp inversion result. It is obvious that the inversion result is not capturing the detail available in the well log.

3. Interpretation becomes difficult because there are multiple models that could result in the same seismic response. Therefore, interpretations are not unique and the uncertainty of any interpretation can not be gauged.

While none of these issues can be solved using deterministic means as there are an infinite number of models that could produce the measured data, it can not be said that each of the infinite possibilities is equally probable. The well data at high resolution can also constrain the inversion.

While the log data is known to be accurate at the well location, it provides no accuracy away from the well. However, geostatistics can provide probabilistic constraints away from the well. Geostatistical methods can not provide a unique model of the subsurface away from hard data points, but can provide probabilistic models or multiple possible realizations (i.e. from methods such as sequential Gaussian simulation). Thus geostatistical inversion proposes to solve for a set of models each of which are constrained by the available seismic data, the well data and the known geostatistics of the study area. Non-uniqueness is not eliminated but the range of possibilities and their probabilities can be estimated.

Fig. 08
Figure 8. The deterministic Zp inversion results with the Zp well log superimposed at “logging” resolution.

Geostatistical inversion

The JGW StatMod MCTM algorithm uses the wavelet that was used in the deterministic inversion. The algorithm uses a Markov Chain Monte Carlo method to sample from the complex probabilistic density function that results from the product of the probabilities generated from the given seismic data, the well data and the geostatistical probabilities. The result of the inversion is a set of high resolution model realizations each of which match the seismic and the geostatistical constraints that are input to the inversion. The results are a representative sampling of the Bayesian posteriori probabilistic density function.

Histograms of the Zp values for each of the two facies (coal and clastics) were generated from the well logs (Figure 9). Also, the vertical variograms for the facies and Zp values were modeled using the well data (Figure 10). The Zp variance was modelled by an exponential variogram and the facies variance was modelled as a Gaussian variogram. The range of all vertical variograms was approximately 2 ms.

Fig. 09
Figure 9. Zp histograms for each facies.

A lack of lateral data made it impossible to empirically model the lateral variograms, so the range and model were estimated by running the inversion repeatedly until the results matched an intuitive belief about the variation of the facies and the Zp values given the regional geology and demonstration cross-sections such as that shown in Figure 2. It was found that when the lateral facies variance was again modelled with a Gaussian variogram and the Zp variance modelled with an exponential variogram, the results were robust as long as the range was set between 200 – 800 m. For the final inversion the range was set to 400 m.

Fig. 10
Figure 10. The vertical variograms for Zp and facies for each of the two facies.

Inversion unconstrainted by the well

The inversion can be run with or without the Zp well log used as a constraint. Initially, the inversion was conducted without well constraint in order to understand how much detail was being provided by the seismic data and the geostatistics alone. The inversion was run to find 10 model realizations. Figure 11 shows cross-sections of the resulting Zp estimates. The panel on the upper left is the average Zp estimate after generating 10 model realizations from the Markov Chain Monte Carlo algorithm. Five of these ten realizations (selected arbitrarily) are illustrated in the other five panels of Figure 11. Superimposed on each of the panels is the Zp well log from Well #2 that has not been filtered and is displayed at "logging" resolution. The inversion results show that the level of detail in the Zp estimates is comparable to the well log resolutions. Also, the absolute Zp values of the inversion are comparable to the absolute Zp values of the unfiltered well log.

Fig. 11
Figure 11. The unconstrained Zp inversion results. The cross-section in the upper left is the mean Zp of 10 inversion realizations. The other 5 cross-sections are an arbitrary sampling of 5 of the 10 realizations.

Each of the 10 realizations should be viewed as possible model solution to the inversion. Figure 12 illustrates the residual difference between the seismic data and the synthetic seismogram generated from one of the 10 model realizations. The figure shows that the residual difference is very small. Each of the 10 realizations shows an equally minimal residual meaning that each is an inversion that is compatible with the seismic data.

Fig. 12
Figure 12. Analysis of one arbitrary realization. (a) The Zp realization, (b) the original seismic data, (c) the residual difference between the seismic data and the synthetic seismogram generated from the Zp model realization and (d) the synthetic seismogram generated from the Zp model realization.

The most likely facies estimate is shown in the panel on the upper left of Figure 13, and the other five panels show the facies associated with the five realizations from Figure 11. The figure shows that the statistical form of the coal zones appears to be modelled in each realization. Each of the individual realizations shows very small coal occurrences in locations that are outside of the expected coal zones, however, when these realizations are averaged to produce the most probable facies, coals are predicted in the known zones only. In the most probable facies panel, continuous but distinct thin coal zones are resolved. However, while the Upper Ardley coal zone appears to be accurately resolved, the Lower Ardley coal zone is not. In the Lower Ardley coal zone, the most probable facies results appear to show two sub-zones but these are together delayed in time and the Silkstone coal zone is not resolved at all.

Fig. 13
Figure 13. The unconstrained inversion facies results. The cross-section in the upper left is the most likely facies of 10 inversion realizations. The other 5 cross-sections are an arbitrary sampling of 5 of the 10 realizations.

Figure 14 summarizes the results of the unconstrained inversion. It shows the deterministic inversion results (for comparison) as well as the most probable facies, the likelihood of coal facies, and the minimum, mean and maximum Zp from the set of model realizations. The probability of coal is high in each of the coal zones and this data can be used to give a quantitative understanding of the likelihood of encountering different gross thicknesses of coal. Again the results are robust for the Upper Ardley coal zone but indicate greater uncertainty in the Lower Ardley coal zone. The data seems to suggest that the existence of the Silkstone coal zone is improbable, since there is a very low probability of coal facies at that level.

Fig. 14
Figure 14. A summary of the unconstrained inversion results. (a) The results of the deterministic inversion (for comparison), (b) the most likely facies, (c) the probability of coal facies, (d) the minimum Zp estimate of the 10 realizations, (e) the mean Zp of the 10 realizations and (f) the maximum Zp estimate of the 10 realizations.

The Inversion constrainted by the well

The inversion was then run with the Zp values and the facies log of Well #2 included as a constraint on the geostatistical realizations. Again, the inversion found 10 possible model realizations that are each constrained by the well data, the geostatistics and the seismic data. Figure 15 shows the residual difference between the seismic data and the synthetic seismogram created from one of the realizations. It shows that the residual is again very small. It is noteworthy that the Silkstone coal zone is a part of the model solution at least at and near the Well #2 location.

Fig. 15
Figure 15. Analysis of one arbitrary realization with the Zp of Well #2 (shown) as a constraint. (a) The Zp realization, (b) the original seismic data, (c) the residual difference between the seismic data and the synthetic seismogram generated from the Zp model realization and (d) the synthetic seismogram generated from the Zp model realization.

Figure 16 shows nine of the ten realizations along with the facies log from Well #2. As in the unconstrained results, the absolute Zp values of the inversion match those of the well log with coal zones showing Zp values as low as 3.0x106 kg/m3*m/s. Each of the realizations shows evidence of the Silkstone coal zone at or near the Well #2 location, but in each realization the low Zp zone of these coals disappears away from the well.

Fig. 16
Figure 16. Nine of the ten Zp model realizations constrained by the Well #2 Zp well log (shown as a facies log for visualization purposes).

Figure 17 shows the mean Zp estimate of the 10 realizations and Figure 18 shows the associated most probable facies generated from the inversion with well constraint. Again, the expected continuity and thickness of the coal zones is produced.

Fig. 17
Figure 17. The mean Zp estimate of the 10 constrained realizations with the Zp well log from Well #2 superimposed.
Fig. 18
Figure 18. The most probably facies realization from the 10 realizations constrained by the Well #2 facies (shown).

Figure 19 shows a summary of the constrained inversion results. The panel in the upper right illustrates the probability of coal facies existence and shows that the probability of coal is highest closest to the well location and Figure 20 illustrates the comparison of the mean Zp and most probable facies results for the constrained and unconstrained inversions. The figure shows that the constrained inversion clearly honors the well data where the unconstrained results do not. However, both are comparable in terms of the coal zone thicknesses and lateral continuity.

Fig. 19
Figure 19. A summary of the constrained inversion results. (a) The results of the deterministic inversion (for comparison), (b) the most likely facies, (c) the probability of coal facies, (d) the minimum Zp estimate of the 10 realizations,(e) the mean Zp of the 10 realizations and (f) the maximum Zp estimate of the 10 realizations.

The Upper Ardley coal zone is well characterized in the constrained inversion; however, the Lower Ardley coal zone character is more problematic. While the Silkstone coal zone is predicted at and near the Well #2, it disappears within 50 – 100 m away from the well. Thus it appears that the seismic data is overriding the likely geostatistical prediction at a distance of less than a quarter of the lateral variogram range (400 m). The likely explanations for this are there is a fundamental flaw with the seismic trace waveform, the wavelet may not be accurate or the time-to-depth table for the well log may be incorrect. However, given that the inversion works very well for the Upper Ardley coal zone, it may be that the problem is localized to the seismic event of the Lower Ardley coals.

Fig. 20
Figure 20. A comparison of the constrained and the unconstrained inversion results showing the most probable facies (left) and mean Zp (right).

A possible explanation comes from Schoenberger and Levin (1974 and 1978) who investigated the effect of short-path multiples on the seismic response of cyclically bedded coal seams. They showed an example (replicated in Figure 21) where the energy of the primary reflection from coals deep in a cyclic series of coal cyclothems was overpowered by the energy from interand intra-bed short-path multiples. The net effect is a delay in the arrival time of the peak energy of the reflections.

Fig. 21
Figure 21. An example illustrating the energy of a primary pulse at time zero and the energy of short-path multiple reflections for a signal that has passed through a series of coal beds. The examples shows how the initial pulse is dominated by the delayed arrivals of pulses from the short-path multiples (after Schoenberger and Levin, 1974).

To investigate whether short-path multiples play a role in the trace wave form of the Alder Flats data, offset gathers were synthesized with the reflectivity method. The method models the full waveform response of the earth and can be calculated with or without including multiples and mode conversions. The results are illustrated in Figure 22 which shows the offset gathers with NMO removed and the stacks of the modelled gathers. In the stack that includes all multiple and mode conversion reflections, there is a slight phase delay (3-5 ms) in the onset of the reflection event that ties to the base of the Lower Ardley coals. This delay corresponds to the scale of delay that would be required to explain why the inversion result is unable to predict the location of the Silkstone coals.

Fig. 22
Figure 22. Synthetic seismograms showing the effect of short-path multiples and mode conversions by comparing gathers and stacks with and without including multiples in the model. All models were created with a zero-phase 5-10-50-60 Hz Ormsby wavelet. The dotted lines are constant time lines for comparing the phase of the different models.


The geostatistical inversion successfully improved the Zp inversion result in the Upper Ardley coals. The absolute estimated Zp values were equivalent to the values from the well logs at logging resolution and bed thicknesses were accurately predicted. In addition, since each of realizations is consistent with the seismic data and is a possible representation of the earth model, the range of Zp values and coal zone net-thicknesses that are captured by the realizations allow the uncertainty of the coal characteristics to be assessed. For example, the probability of achieving a certain coal net-thickness could be calculated from the multiple realizations.



The authors would like to thank Fugro-Jason for donating the use of the Jason Geoscience Workbench to conduct this work. Additionally, the assistance of John Pendrel, Jawwad Ahmad, Brian Hoffe and Olga Rehkopf was very much appreciated. Of critical importance was the support of the CREWES sponsors, the volunteers and workers who contributed to the field acquisition, and the faculty, staff, and students of CREWES and the University of Calgary Geoscience Department. Finally, Hampson Russell Software was used to conduct synthetic reflectivity modelling.

About the Author(s)

Jason M. McCrank earned a B.Sc. in physics from the University of Alberta and a M.Sc. in geophysics from the University of Calgary. Jason’s M.Sc. thesis, completed in 2008, involved geophysical reservoir characterization of a coal bed in Alberta that had been flooded with CO2 for the purpose of enhanced CBM production and carbon sequestration. Jason now works for Shell Canada, Calgary, Alberta, Canada, where he has worked on seismic interpretation and the characterization of heavy oil and tight gas reservoirs. Jason is a member of CSEG, SEG and an APEGGA member in training. In addition to geophysics, Jason’s passions include mountaineering, ski touring and travel.

Don C. Lawton is a Professor of Geophysics and Chair in Exploration Geophysics in the Department of Geoscience at the University of Calgary. He previously served as Head of the Department from 1997 to 2002. His research interests include acquisition, processing and interpretation of multicomponent seismic data, seismic anisotropy, integrated geophysical and geological studies in complex geological settings, and geological storage of CO2. He is an Associate Director of the Consortium for Research in Elastic Wave Exploration Seismology (CREWES) and was recently appointed Theme Lead in Secure Carbon Storage for Carbon Management Canada, a $25M national centre of excellence. He is a past Editor of the Canadian Journal of Exploration Geophysics, and was a recipient of a Meritorious Service Award from the Canadian Society of Exploration Geophysicists (CSEG) in 1996 and the CSEG Medal in 2000. He is a member of SEG, AAPG, EAGE, CSEG, CSPG, ASEG, and APEGGA.

Cheran Mangat was born in Uganda, East Africa and has been a bit of a globetrotter since Idi Amin’s time. She is a fully qualified Geophysicist and has worked in various facets of the oil and gas industry from reservoir engineering, interpretation, seismic processing and reservoir characterization. Cheran joined Fugro-Jason in the UK in October 2001 as a Senior Geophysicist. In April 2004 she transferred to the Calgary business unit. Having been in a technical role for most of her career, she decided to make a 360-degree turn and join the business aspect of our industry. Cheran took the role of Business Development Manager for Fugro-Jason in January 2007. In 2011, Cheran moved to the technical team at GeoGlobal Resources, who is focused in the exploration for and development of petroleum and natural gas reserves in India, Israel and Colombia. She is a member of APEGA, CSEG and SEG.


McCrank, J.M., 2008, MSc Thesis, University of Calgary.

McCrank, J. and Lawton, D,C., 2009, Seismic characterization of a CO2 flood in Ardley Coals, Alberta, Canada: The Leading Edge, No. 7, Vol. 28

Schoenberger, M., Levin, F.K., 1974, Apparent attenuation due to interbed multiples: Geophysics, 39, 278-291.

Schoenberger, M., Levin, F.K., 1978, Apparent attenuation due to interbed multiples, II: Geophysics, 43, 730-737.

Torres-Verdin, C., Victoria, M., Pendrel, J., 1999, Trace-based and geostatistical inversion of 3-D seismic data for thin-sand delineation: An application in San Jorge Basin, Argentina: The Leading Edge, No. 9.

Van der Laan, J., Pendrel, J., 2001, Geostatistical simulation of porosity and risk in a Swan Hills reef, SEG International Exposition and Annual Meeting Expanded Abstracts, September 2001.


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