Integrated shale-gas reservoir characterization: A case study incorporating multicomponent seismic data

Introduction

The Horn River shales in northeast British Columbia contain approximately 70% of Canada’s total estimated shale-gas resource and the largest accumulation of gas in North America (as detailed in the chart in Figure 1). Nexen Energy ULC holds multiple leases in the area, but this case study focuses on a specific project where a multicomponent 3D seismic survey was acquired in 2012.

The resource objectives in the Horn River Basin (see the location map in Figure 2) are the Devonian organic shales of the Muskwa Formation, and the Otter Park and Evie members of the Horn River Formation (Figure 3). In the study area, the top of the Muskwa is at a depth of approximately 1850 m and the top of the Lower Keg River platform (base of the reservoir zone) ranges from 1980 to 2110 m. Commercial quantities of gas can be extracted from these tight shale units using multistage fracturing technology. For this technology to be effective and economic, the wells must be positioned to access high-quality reservoir rock, and well-completion techniques should be optimized based, in part, on the rock properties. The purpose of this 3D seismic attribute derivation and analysis study was to predict rock properties in the seismic data volume, which can then contribute to effective decision-making regarding drilling and completion.

PS data are not as commonly used by interpreters as PP data, mainly due to additional uncertainty in subsurface positioning (laterally and vertically), complications in processing related to accurate estimation of velocities and statics, anisotropic effects that are much more pronounced than in PP data, and most of all, the severe attenuation of high frequencies in the near surface, resulting in considerably lower resolution than in the PP data (Stewart, 2009). For these reasons, PS data on their own are typically neglected as an interpretive tool. Although the PS integration workflow used in this case study has shown promising results elsewhere (Weston Bellman, 2014), a secondary purpose of this project was to determine if any additional benefit to the quantitative analysis of PP data in this play was realized by combining the PS information with the P-wave data.

Quantitative interpretation (QI) is the process of relating quantitative seismic data characteristics, such as amplitude and phase measurements, directly to rock and fluid properties, and using these derived relationships to convert seismic data to geologic predictions (Avseth et al., 2005; Weston Bellman, 2006; Hunt et al., 2012). This is a deceptively simple statement; I say deceptive because the execution is complex, sensitive to subjective parameter choices, and dependent upon assumptions (that may, themselves, be inclined to interpretation biases). For these reasons, I informally refer to the process as qualitative QI. To minimize biases, the workflow used in this project was constrained as much as possible by known quantities, and subjectivity was mitigated by parameter scanning and quantitative measures of assessment wherever possible (this kind of constraint and assessment can be difficult in practice when there are few wells with uncertainty of their own).

Of course, at some point, we do want to impose bias and subjectivity because that is where human creativity and intuition are incorporated, enriching the final interpretation. The workflow described here acknowledges and incorporates the idea that there needs to be separation between process (data cleaning, manipulation, and computation) and interpretation (assigning meaning to inert components; bringing them to life).

The main elements in this QI workflow comprise the following general categories:

Quantitative analysis of the PP data in this study was completed first, resulting in comprehensive predictions of facies and properties, such as the Young’s modulus, Poisson’s ratio, and relative brittleness index (Weston Bellman and Leslie-Panek, 2014). Subsequently, the converted-wave component of the 3D seismic data was analyzed, in addition to a 3D vertical seismic profile (VSP), which were, in turn, integrated into the QI workflow (Weston Bellman et al., 2016). This paper will describe the combined workflow developed for multicomponent integration and compare the resulting predictions with and without incorporating the PS data.

Analysis

Well data

To understand geologic properties in the context of elastic properties (seismically measurable quantities), the wells require thorough analysis for quality and reasonable petrophysical interpretation. The first element in the QI workflow, well data analysis, starts with a visual inspection of recorded wireline logs and the computation of elastic-property logs. Figure 4 shows a log layout with the intervals identified in three levels of detail: The group (labeled on the right), the formation (labeled on the left and colored), and the very fine scale unit distinctions identified in the tops. All scales are important for analytical input to and validation of the QI workflow. The logs shown on this display comprise acquired logs, such as gamma-ray, density and sonic, plus elastic-property logs (computed from sonic and density logs), including P-impedance, Young’s modulus, and Poisson’s ratio.

The crossplots shown in Figure 5 represent elastic properties computed from log data and colored by geologic properties of interest. For these plots to be meaningful, the data must first be assessed for accuracy and completeness. The well-data provided by the client for this project consisted of 25 wells, of which 13 penetrated the Devonian zones and included sonic information. Five of these wells contained dipole sonic logs, which are important for the computation of most elastic properties. Detailed quality control, including assessing the repeatability and consistency of log measurements in similar lithology between wells, was carried out. In addition, data-distribution plots and crossplots of original and computed curves were inspected to identify any issues regarding the log-data integrity for specific wells. Slight variations in the relative positions of the log distributions between certain wells highlight data that fall outside the average (Figure 6). Without a reasonable geologic explanation for the discrepancy, the anomalous data require petrophysical adjustment or are excluded from crossplots and used with caution in the seismic calibration.

Quality-filtered log data from multiple wells, therefore, are plotted on the final crossplots, for the purpose of observing regional or well-to-well variations that can be attributed to geology, not data issues, and for establishing high confidence in geologic property trends.

The essence of QI is to use these distinguishing relationships to predict geologic attributes using elastic properties that are derived from seismic data, a process that is described by several authors in a variety of play types (Goodway et al., 1997; Andersen and Gray, 2011; Iverson et al., 2013). The challenge of QI is to derive the elastic properties from seismic data accurately enough and precisely enough to use the relationships determined from well logs literally.

Seismic data

In practice, of course, the literal translation between well data and seismic data is never perfect. Seismic data are very different from log data. Although representing the same properties of the same rocks, they are acquired in completely different ways — in different frequency bands, with different resolutions and different scales. Log data contain detailed vertical information and very poor horizontal information; what seismic data lack in vertical precision, they makes up for with lateral sampling and sheer volume of information, especially when multicomponent data are included. In addition, the workflow required to determine a property, such as Poisson’s ratio from seismic data is not nearly as straightforward as the simple computation from log data. Multiple steps in the QI workflow, which are each inherently uncertain, are necessary to achieve a reasonable estimate.

We have to accept that seismic representations of elastic properties are always averaged responses with large error bars; similarly, log data are also subject to uncertainty and error. A workflow, such as the one described in this paper, that allows a periodic “manual override” or interpreter judgment, therefore combines the best of all worlds (literal, approximate, and creative). Nevertheless, the aim of the QI workflow is to make the seismic attributes the best they can be with the best input data, the best assumptions, the best parameter selections, and the best tools; naturally, as fast as we can (to be useful).

Seismic data analysis for QI requires, for all but the simplest quantitative objectives, prestack data, which are designed for sampling redundancy in one aspect (midpoint). In other aspects, such as offset and azimuth, PP prestack data provide a rich source of varying perspectives on that midpoint, or spatial location in the earth. PS data provide similar perspectives, but there is more ambiguity in the precise location of the reflections.

With this kind of redundancy, and in the absence of absolute values; trends, least-squares solutions and reasonable constraints are necessary. Much of the constraining information comes from the wells in the form of $V_{\mathrm{P}}/V_{\mathrm{S}}$ relationships, low-frequency models of P-impedance, S-impedance or density, time-depth relationships, and attribute dependencies between the main elastic properties. Appropriate yet cautious use of well data in these steps is critical to the outcome. Additional discussion about the use and pitfalls of well-data analysis can be found in Weston Bellman (2015).

The multicomponent stacked seismic data available for this study are displayed in Figure 7. Comparing the two wave modes — PP and PS — visually reveals differences in character, frequency content, and amplitude variation. The prestack data are shown in Figure 8. The vertical time scale of the PS data in these figures has been adjusted to approximately match that of the PP data. Nevertheless, it is apparent that the PS data are, even at this condensed scale, of a lower dominant frequency. The amplitude spectra of each stacked data type (in the original time scales) are shown in Figure 9 and highlight the significant difference in original frequency content. An interesting feature to note (Figure 8) is that although there are obvious multiples present in the prestack PP data (these are complex residual interbed multiples remaining after several passes of multiple removal), the PS data are “clean” in that regard. This observation alone provides encouragement that there may be unobscured signal in the PS mode. But, the question remains if that information will be sufficiently valuable in the low-resolution context of this data type.

In addition to potential shortcomings related to the low PS resolution, another source of error in assimilating converted-wave data is the uncertainty in choosing joint inversion parameters. This uncertainty can be compounded by inaccurate time-registration between the two data types, which is a prerequisite for joint inversion. The workflow designed for this project (described in more detail by Weston Bellman, 2014) mitigates these issues by first inverting and deriving prestack and poststack attributes on the separate volumes in their native time domain, depth registering these attributes using the velocity information from the VSP, and, finally, combining the PP and PS attributes in multiattribute analyses in the depth domain.

To accomplish the first step in this joint-analysis procedure, AVO attributes, including P-reflectivity (Rp), S-reflectivity (Rs), and density reflectivity (Rd, representing the relative change in density across an impedance boundary), were computed from the prestack, AVO-compliant, PP PSTM gathers using and comparing several methods (Aki and Richards, 1979; Shuey, 1985; Fatti et al., 1994). Similar prestack PS attributes — as approximations to Rs — were computed from full and offset-limited stacks of the PS gathers.

The dipole sonic logs, which were used indirectly for prestack data calibration and low-frequency S-impedance models in the PP workflow, were also important at this stage for synthetic seismogram (synthetic) character comparison and matching with the PS wave data, and shear velocity modeling. For wells without the benefit of a VSP, the sonic and density logs were relied upon for the determination of appropriate synthetic character ties and depth-time relationships with the PS seismic data. Because the PS seismic traveltimes represent a combination of P-wave velocities ($V_{\mathrm{P}}$) and S-wave velocities ($V_{\mathrm{S}}$), the shear-velocity log cannot be used alone to represent the appropriate velocity for creating the PS synthetic. Consequently, a combined velocity log (Vps) was derived using the following function (Todorov and Stewart, 1998):

Synthetic seismograms were created for each well with dipole logs, from the density and derived Vps logs. This method represents zero-offset reflectivity, which strictly speaking does not exist for PS data. The more conventional method for deriving PS synthetics is to model the converted-wave traces for a specific range of angles and stack the synthetic gather to create a synthetic trace for correlation with the real PS data (Lawton and Howell, 1992). However, Figure 10 shows the tie at the VSP well using the above equation, which with minimal adjustment, matches the PS seismic character and the PS corridor stack (with a significantly higher frequency wavelet) very accurately in relative amplitude and phase. Furthermore, the graph shown in Figure 11 illustrates the empirical relationship between the reflectivity obtained using the two methods for multiple angle ranges of the modeled synthetic trace. The relationship is not quite linear, but close enough to justify the use of the equation as a reasonable approximation for a relative reflectivity PS synthetic. In fact, there may be more error introduced by uncertainty in the choice of maximum angle for the modeled approach than that of the linear scaling approximation.

These well ties were used to accurately identify geologic events on the PS data, which were necessary for horizon mapping, event correlation between the PP and PS data, and the creation of accurate velocity models for depth registration.

The various AVO-attribute stacks were then converted from reflectivity volumes to impedance properties using model-based inversion algorithms. A prestack inversion was also performed on the PP data during the PP-only analysis phase. These inversions made use of appropriate low-frequency models derived from log data (inversion examples are shown in Figure 12). There are now six S-impedance derivations to evaluate: prestack inversion of PP data, inversion of the AVO Rs from PP data, poststack inversion of the PS PSTM stack, and inversions of three ranges of offset limited stacks of the prestack PS PSTM gathers. How do we decide which version comes closest to the truth? We can compare with well data of course, but there are only a few dipole sonic logs from which to derive S-impedance logs. Fortunately, in this particular study, we can also compare with inversions of the 3D multicomponent VSP (which is described in more detail later) to get some direction and confidence in the correct choice.

The PP and PS workflows to this point had been conducted in parallel, but in their separate original time domains to avoid compounding uncertainty introduced through the time-registration step. Jointly calibrated velocity models were constructed using time-depth relationships derived from the multicomponent VSP and P- and S-wave velocity measurements from dipole sonic logs. These models were used to convert all computed attribute volumes from their separate time domains to a common depth domain. When registered in depth in this manner, the calculated elastic properties can be correlated to each other, and to well logs, with confidence. This is beneficial for crossplots and mathematical multiattribute analyses.

Additional attributes were computed from various combinations of the AVO/inversion attributes, such as $V_{\mathrm{P}}/V_{\mathrm{S}}$, Young’s modulus $E$, Poisson’s ratio $\nu$, and a relative brittleness index (Rickman Mullen et al., 2008) calibrated to the field-specific minimum and maximum of $E$ and $\nu$. These attributes were derived separately, compared, and used in the classification described later, with and without using the PS data inputs.

To carry out the combined PP and PS analysis, all seismic attributes in the depth domain were included in a mathematical multiattribute procedure to predict density (Russell et al., 1997). Density is difficult to predict directly from seismic data, yet it is a useful geologic property on its own and a fundamental component of other derivative elastic properties such as Young’s modulus. Figure 13 shows the density prediction resulting from the multiattribute analysis using only PP attributes compared with the updated result that includes PS attributes. Although the PP-only result was considered accurate (compared with the match with density logs at wells), it is apparent that by including PS attributes, it is indeed possible to increase the resolution and accuracy of the prediction. At first glance, it looks like this may be a noisier result, but the additional detail is corroborated by similar detail in the S-impedance derivations (particularly when the midoffset inversion is used), and the comparison with the well logs is consistently better. Furthermore, the next section, describing the use of the 3D VSP, illustrates that the geology is, as usual, more complex than originally thought.

3D VSP inversion and rock-property derivation

Analysis of the multicomponent 3D VSP, including inversion, was done at the same time as the seismic workflow to create high-resolution, arguably more accurate (due to the general belief that VSP data are more reliable than surface seismic data), impedance estimates for comparison with and validation of the surface-seismic equivalents.

The offset data for the 3D multicomponent VSP, acquired at a key well with dipole sonic logs, were interpreted and inverted for high-resolution imaging and comparison with the surface seismic. Trace spacing on the migrated 3D VSP was 7.5 m or one-fourth the spacing of the surface seismic (30 m) trace interval. Figure 14 shows a profile through the 3D PP VSP compared with the VSP corridor stack and equivalent surface seismic data at the same scale. There is an obvious increase in detail and resolution visible on the VSP, revealing structural and stratigraphic features that are almost invisible on the surface seismic. The wavelet extracted from the 3D VSP data confirms frequencies up to 100 Hz, which helps account for the temporal resolution, but the increased spatial sampling is also a significant factor in improving the VSP image.

The 3D PS VSP (reverse polarity) in PP-time is shown in Figure 15 compared with the PS surface seismic data. In spite of apparent migration artifacts, a similar increase in resolution is apparent in this comparison because the PS VSP has equivalent frequency content of the PP VSP of up to 100 Hz, but the converted-wave seismic data contain maximum frequencies of only 55–60 Hz (when adjusted approximately to PP time).

A synthetic was created for the VSP well, as described earlier, and correlated to the PP-VSP and key horizons were interpreted based on the synthetic tie. A comparison of one of those horizons (the Mid-Devonian carbonate), between the VSP and surface seismic interpretations is shown in map view in Figure 16, illustrating the structural and stratigraphic detail contained in this high-resolution data. An offset in the structural features is apparent between the VSP and the surface seismic; the reasons for this offset are unclear and, at the time of analysis, were still being investigated by the VSP processors.

The main horizons were smoothed and used to guide a low-frequency P-impedance model derived from the well logs for inversion of the PP-VSP. Similarly, using the S-impedance log computed at the well location and key horizons interpreted on the PS-VSP data, a low-frequency S-impedance model was derived and used to constrain the inversion of the PS-VSP in PP-time. Additional attributes were computed from the derived P-impedance and S-impedance, including $V_{\mathrm{P}}/V_{\mathrm{S}}$, lambda*rho and mu*rho volumes.

Figure 17 shows the comparison between the P- and S-impedance attributes for VSP and surface-seismic data. Interestingly, the S-impedance derived from the VSP has comparable character and resolution with the inversion of the surface seismic midoffset PS data. Interesting, because the PS-AVO equations (Aki and Richards, 1979) describe the influence, relative to angle, of S-impedance and density properties of the rocks (the P-impedance influence is negligible). Models (Figure 18) show that the mid-range angles (or midoffsets) are most affected by S-impedance and that the density influence increases in the far-offsets (Reine and Tilson, 2016). Moreover, the far offsets are the most likely to be affected by stretch and residual NMO (a little trickier to predict and resolve in PS processing) and there is very little S-wave conversion on the near offsets. Hence, the midoffsets do seem to contain the optimum shear impedance information, as we discovered empirically.

Because of this evaluation, the midoffset S-impedance was chosen for the subsequent computation of elastic properties for the final step of the process: correlation to geologic properties using crossplot classification, described below.

Reservoir classification

The quantitative process is now complete and the qualitative interpretation can begin (although certain parts of the process did require some interpretive judgement). In this stage of the workflow, everyone (geophysicists, geologists, engineers, and managers) can provide knowledge and experience to investigate geologic concepts and consequences, which are constrained by known conditions at existing wells and by the various quantitative products derived to this point. Using a combination of well-data crossplots and interpreter intuition, conceptual insight and observation, meaning can be assigned to the various elastic properties derived from the seismic data.

This task starts by circling back to the relationships derived on the crossplots of well data, which will guide the interpretation of the seismic attributes. We can now overlay the seismically derived elastic properties on the equivalent well-data plots and define categories, using cutoffs or polygons, of multiple seismic attributes that directly represent geologic properties. Figure 19 shows the well data and seismic data on a series of layered crossplots. The clusters and trends displayed by the seismic attributes are similar to the well data; however, as described earlier, the literal application of log-derived relationships to the seismic data is generally not ideal due to inaccuracies and uncertainties inherent in seismic data and log data, despite our best efforts. Therefore, an interactive method is introduced that allows appropriate adjustments and calibration of the predicted geologic property categories to match known conditions at well locations.

The formations of interest are first identified and isolated using the pair of attributes most effective at separating large-scale formations (Figure 19a). These formations of interest can then be further subdivided by unit (Figure 19b) on a different pair of attributes. Finally, each unit was assigned a property range (brittleness) using a third pair of appropriate attributes, in this case, Young’s modulus and Poisson’s ratio (Figure 19c). The classified output volume is interactively updated with each class and subclass, so that the interpretation can be refined by inspection and the match at well locations can be optimized.

Figure 20 shows the final classified volumes (with and without using the PS data) after all data have been integrated and all steps in the process completed and optimized. In this particular version of the “qualitative QI,” geologic intervals have been classified and the Muskwa and Upper Otter Park have been selected for further subdivision by relative brittleness index. The actual geologic intervals and computed relative brittleness index are shown at the well location, confirming that the additional detail apparent in the integrated multicomponent result is reasonable.

Conclusions

The objective of this QI project was to characterize the Horn River and Muskwa reservoir intervals and surrounding rocks as accurately as possible with the highest resolution possible, while simultaneously balancing and minimizing uncertainty in the well and seismic data. A secondary objective was the investigation designed to assess the benefits of incorporating the PS data into the QI workflow, and, if possible, improve the prediction of seismic attributes and bulk geomechanical properties. Multiple sources of seismic data were available for this purpose, including 3D-3C surface seismic data, dipole sonic logs, and a 3D-3C VSP.

The results show significant variation in rock properties within the reservoir zones of interest. There are also distinct differences between certain zones, such as the Evie and the Muskwa/Otter Park, and the Upper Otter Park and Lower Otter Park. Although the 3D seismic detects these variations, the 3D-3C VSP resolves significantly more detail, revealing complex structural and stratigraphic features. Comparisons between the surface seismic and VSP horizon structural interpretation and rock property prediction show potential areas of discrepancy that are being investigated further. Nevertheless, the complexity revealed by the 3D VSP should be considered when interpreting surface seismic data in general.

Incorporating the PS data into the QI workflow shows improvement in accuracy and resolution, in spite of initial skepticism that this lower frequency, more uncertain wave mode could add value. Consequently, valuable additional detail is indeed apparent in the combined facies and property predictions, in some of the layers of interest, which contributes to greater resolution and confidence of rock-property predictions (as assessed by the multicomponent VSP and the small number of well ties). However, this apparent improvement can only be definitively confirmed with additional accurate log data (including measured dipole logs) and results from future drilling.

The case study described in this paper has shown that with thoughtful and comprehensive data collection in advance of a shale-gas development project, detailed integrated analysis can add considerable value to the characterization of the geologic properties necessary for efficient development. With attention to detail in well analysis and accurate seismic attribute derivation, calibration and classification, PS data have been merged with PP data to provide improved results in the QI workflow — this in spite of the apparent shortcomings and array of uncertainties in converted shear data.