Deepwater turbidite reservoirs hold some of the largest petroleum reservoirs and thus are important exploration targets. These reservoirs occur in the more distal reaches of the system where the suspended sediments are often very well sorted by their grain size, resulting in high porosity. In the absence of well control, prediction of lithology is based on our understanding of the environment of deposition. By identifying and mapping the diverse architectural elements of the turbidite system and placing them within the correct geologic framework, a skilled interpreter can predict which components of the system are more likely sand or shale prone. Seismic data and seismic attributes also provide insight into the connectivity or compartmentalization of different parts of the reservoir which can be used to estimate the number of wells needed to drain the reservoir.
The deepwater Taranaki Basin located off the northwest coast of New Zealand (figure 1) is covered by several high-quality 3-D seismic surveys that have been made available to the public by the New Zealand Petroleum Ministry. The Romney prospect in this basin has a closure area of 200 square kilometers in 1,600 meters of water depth and it has been estimated to contain 1,100-1,650 million barrels of oil-in-place or 1.7 to 2.7 trillion cubic feet of gas. Although drilling in New Zealand is currently on hold as part of their effort to mitigate climate change, the release of the data to public has resulted in the Taranaki Basin being one of the most widely studied basins in the world.
When 3-D seismic data image the complex internal geometries of channel systems as they do in this survey, they provide an important resource for studying channel fills and their associated lateral distribution patterns. We showcase our efforts at characterization of this deepwater channel complex in terms of seismic attribute applications, specifically using the unsupervised machine learning seismic facies classification techniques. Features such as levees, channel-fill and overbank deposits are more clearly defined on composite seismic attribute displays than on any single attribute by itself. In this study, we describe the turbidite facies and their geometry in detail in order to facilitate a better understanding of the channel architecture.
Only one well is drilled in the area and falls over the seismic volume but does not intersect the Miocene channel. However, this well provides critical lithologic information. The lack of additional well control precludes a quantitative interpretation. For this reason, our interpretation is based on the seismic response to the subsurface and our understanding of geologic processes.
Geologic Background
The deepwater Taranaki Basin developed as a result of multiple tectonic cycles that included the Mesozoic rifting, Late Cretaceous and Paleogene subsidence associated with seafloor spreading, and large-scale Neogene channel systems which transported excessive amounts of sediments into deeper water areas.
The existing evidence and basin models calibrated to the adjacent shallower water Taranaki Basin wells suggest that the expulsion of petroleum from the deepest parts of the basin began in Late Cretaceous and continues even today. Consequently, rocks older than 100 million years are considered to be the geologic basement. The Gondwana breakup from Early to Late Cretaceous initiated the development in the basement started with the opening of the Tasman Sea associated with the subduction of the Pacific plate along the eastern Gondwana margin. Coal-rich formations were deposited from the Cretaceous to the Paleogene in shallow and non-marine environments, which represent the source rocks in Taranaki Basin. The Late Cretaceous to Early Pliocene also saw the deposition of sandstone formations because of widespread marine transgression, which was followed by marine silts and muds. The marine shales deposited during the Early Cretaceous and Paleocene provide additional source rocks. The mudstone deposition took place from Cretaceous to Eocene as well as during the Miocene. These mudstones and the Oligocene to Miocene carbonates (when not fractured) are believed to provide a seal to the source rocks.
Significant uplift of the New Zealand North Island landmass occurred throughout Miocene, and coupled with climate-driven erosion, led to the development of the Neogene major submarine system that transported large volumes of sediment from onshore areas through the deepwater Taranaki Basin into the southern margins of the New Caledonia Basin (figure 1).
Available Seismic Data
The available seismic data were the Romney 3-D survey acquired by Anadarko New Zealand in 2011 made publicly available by New Zealand Petroleum and Minerals. We analyze a 700 square-kilometer subset that falls in water depths exceeding 1,000 meters. The bin size for the data is 12.5-by-25 meters, with a sample interval of 4 milliseconds. Although the Romney-1 exploration well drilled to a depth of 4,575 meters below the sea surface, it does not intersect the Miocene deepwater channel, but it provides useful lithology information. Figure 2 shows that the Miocene strata in Romney-1 includes claystone, siltstone and limestone. SEG polarity convention has been used for seismic data display, i.e., a positive reflection coefficient represents a peak. These lithoelements, when tied to the seismic data, exhibit strong amplitude, high continuity reflections for limestone to mudstone interfaces, and weaker amplitudes internal to the mudstone deposits. Bright seismic amplitudes seen within the channel are associated with shingled reflections indicating sandy fill during channel migration.
Conditioning Seismic Data for Attribute Analysis
Traditional seismic data are seen to preserve information with a frequency content going up to 50 or 60 hertz at the high end of the bandwidth. While such bandwidths might be acceptable for thicker conventional reservoirs, they could lack the needed resolution and be unacceptable for deepwater turbidite reservoirs. The advancements in seismic data acquisition and processing, coupled with computer capacities and speeds, provide cost-effective solutions for such objectives. For a decade-old vintage seismic dataset we need to apply amplitude-friendly poststack spectral – a balancing procedure to enhance the vertical and lateral resolution. The spectral balancing procedure of choice in this exercise is the method based on the peak frequency approach (discussed in the May 2014 and March 2015 installments of Geophysical Corner).
In this method the seismic data are first decomposed into time-frequency spectral components followed by the averaging of the power of the spectral magnitude over all seismic traces in the data volume spatially in the given time window. This step yields a smoothed average power spectrum. Next, the peak of the average power spectrum at each time sample is computed. Both the survey average time-variant power spectrum and its peak are used are used to design a single time-varying spectral-balancing operator that is applied to each trace in the data. Such spectral balancing is amplitude-friendly because a single time-varying wavelet is applied to the entire data volume.
Figure 3a and b shows a comparison of segments of seismic sections traversing the deepwater channel complex before and after frequency balancing, along with their frequency spectra. Notice the well-defined appearance of the reflections both outside and within the channel complex in figure 3b as well as the flattened appearance of the frequency spectra as compared with the input seismic volume.
This data volume was then put through structure-oriented filtering and attribute computation. To bring out the advantage of frequency balancing and structure-oriented filtering, we compute relevant attributes on the data before and after the two data conditioning processes and compare the results in the next section.
Coherence Attribute
Encouraged with the higher-frequency content of the seismic data, we first generated the multispectral energy ratio coherence attribute.
The energy ratio-based coherence algorithm is an eigenstructure-based coherence computation, where the original trace and its Hilbert transform are used to construct the covariance matrix, thereby minimizing artifacts at zero crossings, and providing the ability to use shorter analysis windows resulting in reduced stratigraphic mixing. Coherence run on band-pass filtered or spectral voice components often delineate edges at or near the tuning frequency of a given formation. In general, shorter, more vertically-limited faults and channel edges are often better delineated at higher frequencies, while through-going faults are often better delineated at lower frequencies. In addition to computing a covariance matrix from each bandpassed filtered version of the data and then computing coherence for each one, we can also add the covariance matrices and compute coherence from its sum, giving rise to multispectral coherence computation (discussed in the July installment of Geophysical Corner).
Figure 4 shows a time slice at t=3.18 s from this coherence volume. The channel features seen on the multispectral coherence display are well defined overall and therefore very clearly seen and exhibit higher signal-to-noise ratio. The best practice for stratigraphic objectives is to analyze stratal or even geochronostratigraphic slices through the data. In our example kthe structural dip is low such that the stratal slices are almost parallel to the time slices at the level of the channel. For this reason, we use time-slice displays to minimize any bias due to errors in the horizon picking.
In addition to multispectral coherence, we generated the following attributes for use in multiattribute analysis:
- Relative acoustic impedance is computed by continuous integration of the original seismic trace with the subsequent application of low-cut filter. The impedance transformation of seismic amplitudes enables the transition from reflection interface to interval properties of the data, without the requirement of a low-frequency model.
- Sweetness is a “meta-attribute” or one computed from others, which in this case is the ratio of the envelope to the square root of the instantaneous frequency. A clean sand embedded in a shale will exhibit high envelope and lower instantaneous frequency, and thus higher sweetness, than the surrounding shale-on-shale reflections.
- Spectral magnitude is the magnitude of each spectral component ranging within the seismic bandwidth of the data at a specific frequency increment obtained through spectral decomposition of input seismic data.
- Total energy is computed by measuring the energy associated with the sum of eigenvalues computed from the covariance matrix of the windowed seismic traces and helps isolate low energy chaotic reflectors from higher energy seismic responses.
- Reflector convergence is a relatively less common volumetric attribute, which is useful in the interpretation of angular unconformities. It is a measure of the vector change in reflector normal within a window, thereby mapping bed thinning and thickening, including the onlaps and offlaps used in seismic stratigraphic workflows applied to large workflows to large 3-D seismic volumes (discussed in the January 2012 installment of Geophysical Corner).
In figure 5 we show a segment of a vertical section from the seismic data volume depicting the main channel signature, and we try and correlate the individual signatures with the channel definition seen on the multispectral coherence time slice at t=3.18 seconds. This level is shown with a dotted line on the seismic display. Notice how the original channel position is marked with a black dashed outline and has been shown correlated with the channel definition on the coherence display. The gradual lateral migration of the channel can also be interpreted as shown with the red dashed lines on both the seismic and coherence displays. Similarly, the present-day position of the channel levees is shown with cyan dashed lines on both displays. Thus, the coherence attribute helps understand the original position of the channel as well as its lateral migration up to the present-day position. In addition, the equivalent attribute vertical displays as labeled, corendered with the seismic amplitudes in gray scale, have been included in the display. The individual attributes pick up different features of the channel architecture, and we expect all such component features to come together when we churn through them in machine learning applications discussed next.
ML Applications for Unsupervised Seismic Facies Classification
Machine learning uses mathematical operations to learn from the similarities and differences in the provided data and make predictions. Besides the supervised and deep learning machine learning techniques, there are two broad families of unsupervised machine learning algorithms. The first algorithm family includes dimensionality reduction algorithms such as PCA and ICA (discussed in detail in the August 2018 installment of Geophysical Corner) where N≥3 attribute volumes are represented by two or three component volumes. When plotted against a 2-D or 3-D color bar, the interpreter may be able to identify clusters, but the algorithm output is a continuum of data projected onto the lower dimensional hyperplane.
The second, unsupervised classification algorithm family attempts to explicitly cluster the data into a finite number of groups that in some metric “best represent” the data provided. Two methods that fall under this category are the self-organizing mapping and generative topographic mapping. In both algorithms the data are represented by neurons that lie on a deformed 2-D or 3-D manifold that fits the N-dimensional data where each attribute represents a dimension. In this manner, clusters that are near each other in N-dimensions will also be near each other when projected onto the deformed manifold. We show the application of GTM method to the Romney 3-D seismic volume, on both the input and frequency balanced versions.
Generative Topographic Mapping
The GTM method begins with an initial array of K N-dimensional Gaussian distributions. The Gaussians are centered on grid points mk with a fixed variance 1/β arranged on a 2D plane that best fits the data. These Gaussians represent the likelihood of any data point xr. The goal of GTM is to maximize this likelihood where at each iteration, the variance 1/β is decreased and the Gaussian centroids mk moved to better represent the data. Moving the centroids deforms our initial plane into a 2D surface (called a manifold in mathematics) that better represents the data. After convergence, GTM measures the probability that any data vector is represented by each of the cluster centers, thereby providing a measure of how confident we are in assigning a voxel to one or more clusters.
Essentially, GTM technique projects data from a higher dimensional space (7-D when seven attributes are used) to a lower 2-D dimensional deformed surface onto which an interpreter can draw polygons to form geobodies or project against a 2D color bar for better visualization. After the application of GTM to the data at hand, figure 6 shows the equivalent stratal display to the ones shown earlier, extracted from the crossplot volume generated with GTM-1 and GTM-2 volumes using the 2-D color bar shown alongside. Figure 6a shows the display for the GTM computation carried out when the input attributes were computed from the original seismic data volume. Spectral balancing improves vertical resolution and provides tuning effects that are more closely related to geology, removing much of the overprint of source wavelet. Figure 6b displays the equivalent display when the attributes were computed on the frequency balanced seismic data. Clusters seen on the display in figure 6b are better defined in terms of more colour separation and distinct definition, than the ones seen in figure 6a.
Conclusions
The seismic attribute computation carried out on the input seismic data as well as the frequency-balanced seismic showed more detailed definition of the channel complex features on the latter. Some of these attributes have revealed channel details not seen earlier. The workflow adopted for carrying out seismic facies classification through unsupervised machine learning by way of GTM application has clearly indicated the lateral variability of the channel complex. Though the analysis is qualitative at present, it paves the way for more detailed work as more well and other data become available.