Inversion and interpretation of seismic-derived rock properties in the Duvernay play

We have developed an interpretive seismic workflow that incorporates multicomponent seismic inversion, guided by structural mapping, for characterizing low-permeability unconventional reservoirs. The workflow includes the determination of a calibrated time-depth relationship, generation of seismic-derived structural maps, poststack inversion, amplitude-variation-with-offset analysis, and PP-PS joint inversion. The subsequent interpretation procedure combines structural and inversion results with seismic-derived lithologic parameters, such as the Young’s modulus, Poisson’s ratio, and brittleness index. We applied this workflow to a 3D multicomponent seismic data set from the Duvernay play in the Kaybob area in Alberta, Canada. Subtle faults are discernible using isochron maps, horizontal time slices, and seismic stratal slices. Fault-detection software is also used to aid in the delineation of structural discontinuities. We found that seismic-derived attributes, coupled with structural mapping, can be used to map reservoir facies and thus to highlight zones that are most favorable for hydraulic-fracture stimulation. By imaging structural discontinuities and preexisting zones of weakness, seismic mapping also contributes to an improved framework for understanding the induced-seismic-


Introduction
Unconventional plays represent resource fairways, including those characterized by low-permeability organic-rich rock formations, that are not economically producible by conventional drilling and completion methods. Low-permeability unconventional plays, such as the Duvernay play in the Western Canada Sedimentary Basin, are now routinely developed with multiple horizontal wells drilled from a single surface location and completed using hydraulic-fracturing technology to enhance permeability and enable economic rates of hydrocarbon production (Dusseault and McLennan, 2011). In planning an unconventional drilling program, well logs, cores, and seismic data can provide valuable information for horizontal well placement and the design of hydraulic-fracture stages. Perez and Marfurt (2015) describe the importance of the use of seismic to provide estimates of lithologic parameters such as the Young's modulus, Poisson's ratio (PR), and brittleness.
Prior to the widespread development of unconventional resource plays, Pennington (2001) introduces seismic imaging as a primary tool to deliver statistical constraints for reservoir development, within an emerg-ing framework of reservoir geophysics. In contrast to traditional approaches for seismic exploration, Pennington (2001) emphasizes how the calibration of 3D seismic models using well data, combined with rockphysics relationships, could be used to differentiate between competing reservoir models on the basis of lithologic character, fluid content, and in situ conditions such as pore (over)pressure. Dommisse (2013) applies similar concepts to unconventional reservoirs and demonstrates how stratal slicing can provide a powerful technique for seismic reservoir characterization. By focusing on lateral spatial variability in the seismic expression of various reservoir facies, Dommisse (2013) argues that stratal slicing can overcome bandwidthrelated limitations in vertical seismic resolution. Goodway et al. (2010) describe methods built on the use of amplitude-variation-with-offset (AVO) analysis to differentiate ductile shale reservoirs from brittle reservoirs in the Barnett Shale of Texas. Parameters known as lambda-rho and mu-rho were derived from well logs, crossplotted, and compared with data derived from seismic inversion (where LMR stands for the Lamé parameters, lambda, mu, and rho). Trends that emerged from this approach demonstrated that variations in rock properties, often attributed to the brittle behavior of reservoir rocks, exhibit coherent patterns that reflect reservoir quality. When seismic-derived LMR attributes are back projected onto a seismic section, trends in brittleness are observed and can be mapped spatially within a seismic volume. Perez and Marfurt (2015) show that this technique can be useful to extract the mineralogical content of the rocks. The importance of fault mapping is highlighted by Refayee et al. (2016) in a Utica Shale example. Parameters such as dip, similarity, and curvature were extracted from the seismic data and incorporated into an interpretive model. A "faultenhanced filtered" seismic volume was derived from the seismic cube (similar to semblance) and used to map major faults and image fracture networks. Refayee et al. (2016) propose that sweet spots within the reservoir are defined by areas of enhanced natural fractures and are mappable within the seismic volume. Meek et al. (2013) highlight the advantages of a multidisciplinary reservoir geophysics approach by combining results from microseismic monitoring, structural attribute analysis, and seismic petrophysics. Comparison of recorded microseismic events with curvature anomalies and other seismic-derived attributes such as Young's modulus revealed a strong correlation with the density of recorded microseismic events. Prospective areas for future development were then identified based on seismic attributes. Similarly, Rafiq et al. (2016) show how sets of attributes extracted from microseismic data can be correlated with curvature anomalies from 3D seismic data to partition a reservoir into depositional facies units. Chopra et al. (2017) conduct a seismic analysis within part of the Duvernay play to show that induced seismicity appears to follow preexisting basement faults. Schultz et al. (2017) investigate areas within the Duvernay play that are prone to induced seismicity from hydraulic fracturing and identify event hypocenters extending from the Duvernay into the underlying Precambrian basement. This study confirmed a correlation between the presence of basement fault systems and the location of induced seismic events. Igonin et al. (2017) use a passive recorded data set (colocated within our 3D/3C seismic survey) to analyze induced seismic and microseismic events in the study area.
Our study makes use of a 3D multicomponent seismic data set in the Kaybob portion of the Duvernay play, within a region where nearby induced earthquakes have occurred (Bao and Eaton, 2016). Analysis of this data set includes structural interpretation and mapping, fault detection, and reservoir characterization based on poststack and simultaneous P-P and P-S inversion. Various presentations of the data are used to illustrate how structure and inversion attributes can be used to identify areas that are most prospective for development, as well as areas where the existing fault architecture may pose a risk for induced seismicity.
The main objective of our study is to present a comprehensive workflow for interpretive processing and inversion of 3D multicomponent seismic data for unconventional reservoir geophysics. We start with the observation of the Duvernay/Leduc system in outcrop, and we use the observations of the outcrop to guide the interpretation. The faulting in the Duvernay Formation can be correlated and interpreted considering basement tectonics. The Gillwood Formation, affected by the same tectonic regime, provides a guide to the tectonic history, giving insight into deep-seated strike slip and vertical offset basement faulting. The study area has a history of induced seismicity (Igonin et al., 2017), which can independently be used to spatially identify fault activation. We apply this approach to a data set from the Duvernay play, demonstrating how seismic data can be used to map structural discontinuities and faults, thus providing insights into the tectonic history of the reservoir. We also show how seismic inversion applied to poststack and multicomponent data gathers (PP and PS) can be used as a tool to identify facies changes, fault boundaries, and potential geohazards such as basement faults.

Geologic Setting
The Duvernay Formation, a bituminous/argillaceous carbonate of late Devonian age in the Western Canada Sedimentary Basin, is emerging as a major resource play in North America (Creaney and Allan 1990;Hammermaster, 2012). The Duvernay is rich in organic matter and, depending upon thermal maturity and position within the basin, it produces gas, natural gas liquids, or oil (Switzer et al., 1994). It is also commonly believed to be the primary source rock for the Devonian Leduc reef, Nisku, and Wabamun carbonate plays (Dunn et al., 2012). With the advent of increasingly widespread horizontal drilling and multistage hydraulic fracturing, the Duvernay is being recognized as a world-class unconventional resource play (Davis, 2013).
The Duvernay Formation was deposited adjacent to several large Leduc reef complexes ( Figure 1). The Duvernay has an effective porosity of 6%-7% and an average total-organic-carbon (TOC) content of up to 4.5% (Chopra et al., 2017). It is correlative with the lower member of the Leduc Formation (Figure 2), and it is believed to be the source rock for most of the Leduc, Nisku, and Wabamun oil pools in Alberta (Dunn et al., 2012). Growth of Leduc reefs was terminated by sea level rise, during which the reef-building organisms were ultimately unable to keep up with the rising sea level and were drowned (M. MacKay, D. Eaton, and C. Clarkson, personal communication, 2018). The Leduc and Duvernay Formations are overlain by the quartz-rich Ireton Formation that forms a cap and seal to hydrocarbon reservoirs.
The depositional environment of the Duvernay Formation varied depending on where it was situated with respect to the Leduc reef margin. Factors that influenced the depositional environment included tides, storms, and sea-level changes. Interreef areas were protected from waves and tides, so that sediments in these regimes were deposited in a low-energy environment. In the subsurface, five lithofacies have been identified from cores (Dunn et al., 2012); these are argillaceous mudstones, bioturbated limestones, organic-rich siliceous mudstones, siliceous organic-rich mudstone, and mixed siliceous mudstones.
Variability in depositional setting is evident in Figure 3, which shows a Leduc reef (Presqu'ile Formation) outcrop characterized by multiple stages of reef growth, with adjacent Duvernay-equivalent (Perdrix Formation) strata. This outcrop exhibits lateral and vertical lithofacies variation within the section, including relationships relative to the reef margin. The reef mass visible in this exposure is 100 m thick and 250 m wide in section. The reef exhibits five distinct episodes of growth. To the north, the Perdrix Formation transitions into a quartz-rich shale. The increase in quartz content has an influence in rock properties; quartz-rich areas would have a higher Young's modulus and a lower PR, and they may be more brittle than the rocks near the reef (Cho and Perez, 2014). This outcrop shows, at a scale similar to features resolved by seismic images, how the Duvernay Formation lithology varies laterally and vertically. It follows that reservoir characteristics such as V P , V S , ρ, and the brittleness index (BI) exhibit similar lateral and vertical variability in the subsurface.
Brittleness has been identified as a key indicator of sweet spots in unconventional plays, i.e., areas that are developed with horizontal drilling and completed with hydraulic fracturing (Cho and Perez, 2014). In general, brittleness is a desirable property that defines which rocks will fracture and yet maintain sufficient strength for the fractures to remain open after placement of proppant. Cui et al. (2017) present various definitions of brittleness, including how they are derived and how the expressions are relevant to reservoir characterization. In this paper, we will use the definition of brittleness proposed by Rickman et al. (2008), which is based on V P , V S , and ρ.
A rock with a high BI will necessarily have a low PR and a high Young's modulus E. Within an established lithofacies framework, a strong correlation exists among quartz content, TOC, and brittleness of the Duvernay Formation (Dunn et al., 2012). The higher silica content of the Duvernay comes from deposition at distal, low-energy areas. Soltanzadeh et al. (2015) provides evidence that clay content also plays a significant role in the brittleness of the Duvernay and Ireton Formations. Moreover, hydrocarbon generation has changed rock properties such that the high-TOC Duvernay tends to be more brittle.

AVO-compliant processing
The 3D multicomponent seismic data used in this study were acquired and processed in 2015. The data were processed in an AVO-compliant manner (Lee et al., 1991;Allen and Peddy, 1993) by using processing steps designed to preserve relative amplitudes for PP and PS data sets. The AVO-compliant processing flow Figure 1. Hydrocarbon maturity windows for the Duvernay petroleum system (modified from Creaney and Allan, 1990), showing the location of the study area. The stratigraphic column on the right shows the Middle and Upper Devonian regional stratigraphic nomenclature, highlighting several units that are discussed here. Interpretation / May 2018 SE3 used to generate the PP data volume is generally similar to a conventional seismic-processing flow, with a few key differences. Spectral balancing was applied in a surface-consistent manner (Taner and Koehler, 1981) by scaling low-amplitude frequency bands using correlated stack equalization. Noise attenuation was then performed in the cross-spread domain (Calvert et al., 2008). Traces with residual anomalous amplitudes (i.e., outliers) were manually edited from the data. Prestack time migration (Fowler, 1997) was then applied, followed by a radon multiple attenuation (Kelamis et al., 1990). Five-dimensional interpolation and normalization (Chopra and Marfurt, 2013) was then applied to the data to fill in missing traces and to normalize the data to the bin center.
Essentially, the same AVO-compliant workflow was applied for PS data processing, with several adaptations that account for the specific raypath geometry of converted waves (Stewart et al., 2002) in addition to S-wave splitting. In particular, the bin centers were determined using common conversion point binning (Eaton and Lawton, 1992). S-wave splitting analysis was performed to determine the time shift and rotation angle between the fast and slow qS arrivals. Amplitude compensation was performed using the approach of Jin et al. (2000), by applying corrections that reverse the S-wave splitting effect along the upgoing raypath. After applying these processing steps, 5D interpolation and normalization was performed, with the PS output bins normalized (interpolated) to coincide with the location of the PP bins.
The version of the stacked data used in this study for interpretation and prestack inversion was limited to a maximum 40°P-wave incidence angle, based on ray tracing using a velocity model derived from a control well (well A) that is discussed below. The corresponding gathers then had residual normal moveout and trim statics (Hoeber et al., 2005) applied to the PP and PS prestack data volumes.

Structural mapping
Sonic and density logs from control well A were used to generate a synthetic seismogram to achieve a precise tie to the seismic data ( Figure 4). During this process, the source wavelet was also estimated. Well A was logged to the Middle Devonian Gilwood member of the Watt Mountain Formation to a depth of 3607 m. This synthetic seismogram tie establishes a time-depth relationship, which was subsequently used to create angle gathers, perform PP-PS registration, and provide a calibration point for the seismic inversion. Figure 5 shows a west-east seismic profile extracted from the PP stacked data volume. The density and sonic curves from well A are overlaid to illustrate the correlation of these log curves with the seismic profile. The Duvernay, Gilwood, and Precambrian seismic horizon are correlated based on the tie to well A. A horizon marking the top of the Precambrian basement was also picked, based on a regional correlation from a deep crustal seismic profile (Eaton et al., 1999) that passes close to this 3D seismic survey. Structural discontinuities are evident as breaks or sharp bends in the picked horizons ( Figure 5). Figure 6 shows a perspective view of time slices from the PP data volume that intersect the Duvernay and Precambrian horizons, respectively. The location of well A relative to the 3D seismic surface is also shown. Both time slices reveal amplitude trends that define a prominent set of north-south-oriented features. In the case of the shallower time slice, terminations of these prominent features occur along a northeastward trend, defining a secondary set of lineaments.
Time-structure maps were created throughout the data volume for all of the horizons listed above. Representative time-structure maps of the second  Figure 7. For all of the Devonian stratigraphic horizons as well as the top of Precambrian basement, the overall structural dip is toward the northwest. For the shallower horizons, this trend changes to a regional southwest dip, which is oriented toward the Rocky Mountain deformation front. The 2WS surface also exhibits a distinctive set of northeast-trending lineations. This pattern is also present on the Duvernay horizon, albeit to a lesser extent, whereas the aforementioned north-south structure trend is strongly expressed. This north-south structural fabric is not seen above the Ireton Formation; at greater depth, it extends through the Gilwood and Precambrian horizons.
A time-structure map on the Gilwood horizon is shown in Figure 8. The Gilwood member of the Watt Mountain Formation is a clastic unit that was deposited during uplift of the Peace River Arch (O'Connell, 1994). Seismic mapping of this unit reveals a conspicuous feature that is interpreted as a meandering river channel that may represent a drainage system that is linked to the uplift of the Arch. Figure 8 also shows images that were created using a coherence-and curvature-based algorithm for fault detection, followed by an adaptive principalcomponent analysis. This software delineates and displays lateral discontinuities in the seismic data and highlights several features that are further enhanced by 3D rendering and shading. These features help to elucidate the meandering channel-like feature, which is cross-cut by a prominent north-south linear feature. As discussed below, we interpret this feature as a possible strike slip fault with left-lateral displacement; based on the channel offset, the net displacement along the fault appears to be approximately 1 km.

Seismic Inversion
The seismic horizons, well ties, and structural discontinuities discussed in the previous section are used to guide the inversion process. The theoretical basis for the inversion procedures was presented by Hampson et al. (2005).

Poststack inversion
Our inversion workflow begins with poststack inversion, in which an initial model of acoustic impedance (i.e., the product of velocity and density) is constructed and then perturbed until a satis-factory maximum-likelihood fit is achieved between predicted and observed stacked seismic traces. Sonic and density logs from well A were used to build a 1D acoustic impedance model, which was then extrapolated and adjusted throughout the data volume by tracking the Precambrian, Gilwood, Swan Hills, Ireton, and Wabamun horizons. This produced the initial acoustic impedance model, which is illustrated in Figure 9a. This model provided low-frequency (0-12 Hz) components that are absent Figure 4. The PP synthetic-seismogram well tie for the interpretive workflow, showing the extracted source wavelet. The integrated two-way times from the V P log at well A have been shifted to match the observed times. Zero-offset synthetic seismograms are duplicated for clarity; those plotted in blue are from the zerooffset synthetic. Those plotted in red are the actual seismic traces at the wellbore. A good fit for the latter synthetic seismogram is evident by overlaying it onto the PP data volume at the location of well A, as shown on the right. This correlation provides confidence for the horizon picks and time-depth relationship. The location of well A relative to the seismic survey is indicated in Figure 6. in the stacked seismic data, due to unavoidable bandwidth limitations. Figure 9b illustrates the poststack inversion results obtained using poststack seismic inversion. A sample set of input traces from the PP data volume is shown in Figure 9c. The model-based inversion procedure uses a linear regression method (Russell and Hampson, 1991) with a linear-programming approach to obtain a maximum-likelihood waveform fit, using the estimated source wavelet from the well tie. The application of poststack inversion increases the resolution of smallscale features, vertically as well as laterally, and also extends the dynamic range of the model (i.e., the difference between the maximum and minimum acoustic impedance values).
Stratal slices can be used to facilitate geologic interpretation for a data volume. This approach is based on extraction of a horizon-parallel slice through the impedance volume, which can then be visualized in map or perspective view. Unlike an image-time slice (Figure 6), amplitude or property variability that is apparent on a stratal slice represents internal variations that occur across a constant geologic-time surface. An example is shown in Figure 10, in which a grayscale rendering is used to accentuate subtle features. This stratal slice was extracted 8 ms above the base of the Duvernay zone, generating an image that reveals features that are broadly representative of the lithologic variability seen in the outcrop image shown in Figure 3.

PP-PS joint inversion
The next step in our workflow procedure involves PP-PS joint inversion. Figure 11 shows a small-scale PP seismic section showing the seismic events used for registration. There were no dipole logs available on the seismic survey that could be used to calibrate the V P ∕V S relationship; there was only one deep well location having a suitable sonic/density log (well A, which was drilled to the Gilwood Formation). The seismic events correlated for the purpose of registration were Colorado, 2WS, Montney, Wabamun, Ireton, Swan Hills, Gilwood, and Precambrian. The same events were correlated on the PS data set and used to perform the PP-PS registration. A crucial step in this process is registration of the PP to PS data, which requires correlation of seismic events in the PP and PS data sets that correspond to the same geologic formation. In practice, if we assume (for the sake of simplicity) that a constant V P ∕V S ¼ 2 is applicable, the PP and PS sections can be easily compared by displaying the PP data at 1.5 times the scale of the PS data. This simple scaling relationship can be derived easily from the formula for average V P ∕V S for the interval between two horizons (Stewart et al., 2002): where Δt PP and Δt PS are PP and PS isochrons (time difference) between horizons that have been correlated between PP and PS sections. If the upper horizon is taken as the surface, then it follows from equation 1 that if V P ∕V S ¼ 2, then the surface-to-depth isochron ratio is Δt PS ∕Δt PP ¼ 1.5. Equation 1 has been used to determine the interval V P ∕V S for various pairs of horizons, as summarized in Table 1. Table 1 calculates V P ∕V S values varying in and approximately 2.0 for the correlated seismic events suggesting the character correlations for registration  were reasonable, and they fall within the theoretical limits for V P ∕V S ratios (Castagna, 1993). Figure 12 shows a large-scale plot showing the result of the PP PS registration, showing the data tie over the zone of interest. Figure 12a shows the PP seismic data, and Figure 12b shows the PS data displayed in PP time.
The approach for simultaneous prestack inversion uses angle gathers for the estimation of P-impedance, S-impedance, and density. The algorithm is based on the assumption that the linearized approximation for reflectivity is valid. In particular, the Aki-Richards equation (Hampson et al., 2005) is used, and it determines the reflectivity as a function of the incident angle. The PP and PS gathers are converted from offset gathers to angle gathers, based on the velocity and ray trace tie from well A. The coregistered PS and PP angle gather data were used as inputs into the joint inversion. The inclusion of the V S data means that the Sreflectivity and impedance is a direct calculation based on the PS section. The outputs of the joint inversion step are V P , V S , and ρ, from which Young's modulus E, PR, and BI were calculated using: where Z S is the S-wave impedance. The BI is a function of E and PR and is given by where w is a weighting factor that controls the relative importance of Young's modulus and PR. Here, w ¼ 0.5 is used. A high BI is considered to be a desired lithologic trait for hydraulic fracture completion (Rickman et al., 2008). The geographical area covered by the 3D/3C survey has very little well control at the Duvernay Formation, with four horizontal wells drilled to date. This area is pres-  Interpretation / May 2018 SE7 ently in the process of evaluation and development. As development proceeds, the weighting factor w may have to be adjusted to better explain reservoir performance. Figure 13 presents stratal slices that depict E and PR extracted from the same Duvernay stratal interval as shown in Figure 10. These stratal slices highlight several distinct features. For example, the northeastern part of the map is characterized by relatively low E and high PR; this combination is indicative of relatively lower brittleness, which is less desirable for unconventional reservoir development. Near the central part of the survey, a region with high E values is bounded to the east by an abrupt north-south-trending edge. This edge appears to be aligned with an interpreted deep-seated fault. To the west of this edge, the maps of E and PR exhibit a high degree of variability with pronounced quasilinear features at various orientations. Figure 14a shows a crossplot of E versus PR for the Duvernay interval, including all data values from the joint inversion volume from the top Duvernay horizon to the base Duvernay horizon. For comparison, the points are overlaid onto theoretical values of E versus PR (Cho and Perez, 2014) as shown in Figure 14b. The Duvernay points cluster within a relatively confined region near the middle of the theoretical values for quartz-clay and quartz-limestone mixtures. The Duvernay values also appear to fall approximately along a trend of constant bulk modulus, as indicated on the graph by a set of black curves.
The crossplot in Figure 14a has been used to group the data into four distinct seismic-lithology zones on the basis of PR. A seismic cross section in Figure 15 displays the PP-PS inversion results, in which each data point is assigned a color according to the parameter ranges of these seismic lithology zones. Although the zones are defined in Figure 14 based on the range of Duvernay values, these have been applied to data points that fall outside of the Duvernay interval. Along this west-east profile, the Duvernay interval is classified as either zone 3 (blue) or zone 4 (red), in which zone 4 may be considered to be more brittle because it is characterized by a higher E and lower PR. The fault appears to be marked by a pronounced lateral transition to higher brittleness on the west side. The potential significance of this transition is discussed below. Figure 16a shows a plot of the BI calculated using equation 3. The major north-south fault is shown, as is the eastwest seismic reference line, and the tie to well A. Warm colors are high in the BI, cold colors are low values. Figure 16b shows interpreted near-vertical basement-rooted faults seen on all horizons from the Precambrian to the Duvernay. This system of faults was picked and correlated based on matching seismic displacements from the Ireton to the top of the Precambrian basement. Combined, these two maps show key elements in evaluating the Duvernay structure controls and lithologic properties.

Discussion
Although standardized resource-play development schemes, such as horizontal drilling on a regular grid and use of uniform completion designs, may be well-suited to homogeneous reservoir conditions, they are not optimal for the type of geologic complexity that is likely to prevail in the case of the Duvernay play. For example, a uniform pattern of wells may fail due to reservoir heterogeneity and the presence of faults, whereas uniform perforation intervals may not be suitable, particularly when a wellbore crosses a facies boundary with varying brittleness. The stratigraphically equivalent Perdrix Formation in Figure 10. Stratal slice showing acoustic impedance variations within the Duvernay Formation, created using a time horizon 8 ms above the base Duvernay reflector and extracting the amplitude from the poststack inversion volume. Some interpreted trends are highlighted on the right. (a) An uninterpreted slice with arrows highlighting prominent features and (b) an interpretation of three of the stratal boundaries. There are several more stratal features seen on these displays. Figure 11. Horizon correlation used for registration of PP and PS data volumes. This display includes all seismic events used for the PP-PS registration, tied to the integrated sonic log.
SE8 Interpretation / May 2018 the southern Canadian Rockies (Figure 3) exemplifies, in a spectacular manner, how changes in the Duvernay lithofacies may be linked to increments of reef growth. Outcrop analysis also reveal how the off-reef depositional environment varies over spatial scales of tens to hundreds of meters, which is compatible with resolution capabilities of seismic data. Stratal slices from 3D seismic inversion volumes, as depicted in Figures 10, 13, and 16, reveal coherent patterns of parameter variability that likely reflect varying reservoir conditions due, in part, to distinctive lithofacies arising from episodes of Leduc reef building.
The Duvernay Formation is also highly structured, and development of this unconventional play must consider structural controls along with depositional heterogenity. Depending upon their present-day permeability structure, faults could represent either a barrier or a pathway for fluids in the subsurface. The presence of faults can thus exert a major influence on reservoir characteristics, even if the accumulated offset is relatively small. The expression of such a fault system in a 3D seismic volume depends, to a large degree, on whether past tectonic activity was syn or postdepositional. Syndepositional faulting could lead to some step-like facies change that is localized across a fault; for example, if the total organic content of the Duvernay Formation varies systematically across a fault, this would lead to changes in E, PR, and the BI (Soltanzadeh et al., 2015). On the other hand, postdepositional faulting could generate seismically detectable linear features that are aligned along a fault. Such linear features could be associated with a damage zone around seismically invisible fault core, or from fluid alteration in proximity to a permeable fault. Hence, three pertinent ways in which an active fault system could affect lateral facies distribution are (1) syndepositional faulting could create a lower energy depositional environment on the downthrown sides of the fault in a marine setting, (2) postdepositional faulting could be manifested by a damaged or alteration zone, and (3) strike-slip faulting could transport material laterally such that different facies are juxtaposed against each other.
This study highlights various imaging approaches that can be used for identi-  Interval and surface-to-depth V P ∕V S ratios determined using equation 1. These values were measured from the well A tie on the PP and PS reflectivity data. The row in blue is the measured PS times, the column in green is the measured PP times. fication and crossvalidation of putative fault systems, based on inversion of multicomponent 3D seismic data. Faults with small vertical offsets can be manifested in profiles as seismically detectable structural discontinuities, or more commonly as fold hingelines ( Figure 5). Such features identified on seismic profiles extracted from 3D volumes can be crossvalidated, on the basis of spatial coherency on time slices or horizon timestructure maps (Figures 6-8). This type of anomaly is also amenable to seismic mapping using curvature attributes (Rafiq et al., 2016;Chopra et al., 2017), as well as automatic discontinuity/fault detection with adaptive principal-component analysis. Lateral offset of paleogeographic features such as a river channel ( Figure 8) provides a tool to measure lateral fault offset, which is otherwise very difficult to quantify in a horizontally layered medium. In this study, seismic images imply that strike-slip faulting was active during or after formation of a Middle Devonian Gilwood river channel.
As argued above, seismic lithology parameters may provide important clues for distinguishing syn-and postdepositional faulting. For example, a step-like facies change across a fault, as evident from PP-PS inversion results for an inferred north-south feature (Figures 15  and 16), is interpreted here to indicate that syndepositional tectonic activity may have taken place during Leduc reef buildup and deposition of the Duvernay. Various linear anomalies that lack such a step-like change are also documented in this study, suggesting that the Duvernay Formation in this area may be transected by a fault system that includes younger (i.e., postdepositional) elements. On the other hand, care in the interpretation of these seismic anomalies is needed because some prominent linear features (Figure 13a) may be artifacts arising from the acquisition footprint.
Defining the fault architecture is thus crucial to well planning. The objective of a fracture stimulation is to inject propant into the Duvernay Formation, open up the fracture porosity, and thus produce hydrocarbons. In the presence of a critically stressed basement fault, the fault may become a thief zone; for example, in a neighboring part of the Duvernay play, Bao and Eaton (2016) find evidence for fluid loss and associated induced seismicity defining  Figure 10. The east-west lineations in Young's modulus to the south of well A are probably artifacts due to an acquisition footprint.  (Cho and Perez, 2014) showing that the Duvernay points fall approximately along a line of constant bulk modulus and cluster near the middle of the range for quartz-clay and quartz-limestone mixtures. The highest BI corresponds to high E and low PR and thus would plot toward the top left corner of the graph. Finally, as noted by Cho and Perez (2014), a desired combination for brittleness is a high E and a low PR. Hence, a desired reservoir is easy to fracture but sufficiently strong to hold the fractures open long enough for placement of proppant to occur. Attempting to induce a fracture in a rock with low E would produce poor results because the rock material would collapse back into the induced fracture, sealing it off. Similarly, attempting to induce a fracture in a high-PR material is more likely to produce ductile deformation, with a lowered tendency for a fracture to form. In the absence of faults, fracture stimulation should thus be targeted toward areas of high brittleness while avoiding more ductile facies. Although straightforward in concept, the use of BI has many complications and is not straightforward. Nevertheless, the PP-PS seismic inversion in this study shows promise for mapping seismic lithology parameters to aid in the development of Duvernay reservoirs.

Conclusion
The Duvernay play is prone to induced seismicity and is likely to exhibit considerable lithologic heterogeneity in the subsurface, at scales that are seismically resolvable. Consequently, we recommend that planning for drilling and hydraulic-fracturing well-completions programs should consider reservoir facies and presentday fault structure. In this study, we have developed a comprehensive workflow for interpretive inversion of multicomponent 3D seismic data, guided by structural analysis. Our workflow requires AVO-compliant data processing and makes use of poststack inversion to obtain an acoustic-impedance volume, followed by correlation of horizons between PP and PS sections and prestack PP-PS AVO inversion. The output parameters from the inversion, consisting of V P , V S , and ρ, are used to calculate data volumes containing Young's modulus E, PR, and BI. Various presentation formats for the inversion results, including perspective views and stratal slices, provide insights for geologic interpretation of the data.
The inversion results reveal spatially coherent patterns of seismic lithology parameter variability that support several key interpretations. Lithofacies variations within the Duvernay Formation likely reflect episodes of growth of proximal Leduc reef complexes. These processes lead to variations in E, PR, and BI, which indicate that some areas within the survey volume are better suited than other areas for hydraulic fracturing. Fault systems that are characterized by relatively small vertical offsets are expressed by structural discontinuities and, more commonly, curvature anomalies defined by fold hingelines. The largest fault is a steeply dipping structure with a north-south strike; this fault produced approximately 1 km of leftlateral strike-slip offset, based on inferred displacement of a Middle Devonian Gilwood river channel. A step-like discontinuity in BI, derived from the PP-PS inversion results, provides evidence that for syndepositional motion on this fault. Other seismic anomalies are quasilinear features that lack such a step-like discontinuity; these anomalies are interpreted as damage or alteration zones near postdepositional faults that belong to a complex fault network in this region.
From an economic perspective, basement faulting may produce an undesirable result when a fracture stimulation is attempted nearby. The objective of a completion program is for the hydraulic fracture system to remain within the formation. The presence of preexisting critically stressed faults can produce the undesireable result of induced seismicity. This comprensive analysis can provide the needed data to effectively develop the Duvernay Formation unconventional play in the Kaybob area of Alberta.

Acknowledgments
We are grateful to Arcis Seismic Solutions for providing the 3D multicomponent data used in this study and for PP and PS AVO-compliant processing that was used for this study. We thank the sponsors of the Microseismic Industry Consortium and Consortium for Research in Elastic Wave Exploration Seismology (CREWES) for their financial support of this study. This work was also funded by Natural Science and Engineering Research Council of Canada (NSERC) through grant nos. CRDPJ 461179-13, CRDPJ 474748-14, IRCPJ/485692-2014, and IRCSA 485691. We thank B. Russell for his technical advice and consultation as well as P. Pederson, staff and students at the University of Calgary for technical assistance. We thank S. Reid for assistance with graphics