Analysis of two 3D surveys, available well data, published outcrop data and subsurface information, as well as production data available from the state of Pennsylvania, demonstrates that wide-azimuth seismic is sensitive to variations in fracturing at the scale of individual pads or even individual wells. These variations in fracturing begin to explain why production varies significantly, even locally, within the Marcellus Shale gas play. Rose diagrams from quantitative fracture analysis using azimuthal seismic velocity volumes were compared with published data from Appalachian black shale outcrops and subsurface fracture models proposed in various papers to validate the results from subsurface data. These analyses provided insight into the rock fabric and the presence of systematic joints that likely affect production. There was a strong correlation between the low anisotropy and low heterogeneity of anisotropy and high estimated ultimate recovery (EUR). Additionally, interpreted fracture trend azimuths differed between areas of larger gas EUR and areas of smaller gas EUR as defined by decline curve analysis. Some perforations were likely to perform much better than others along the borehole, based on observed heterogeneity in the seismic profiles and map view.


Although it has long been understood that natural fracture systems are essential for achieving the best production in Marcellus Shale gas wells, methodologies for verifying the heterogeneities in these fracture systems in the subsurface are not well understood. Analysis of wide-azimuth, P-wave seismic velocity attributes at the reservoir level, and for specific laterals or proposed laterals, can provide this insight. The azimuthal variation in seismic velocities is an ellipse in the horizontal plane and can be characterized by three parameters: the fast velocity (Vfast), a perpendicular slow velocity (Vslow), and the fast velocity azimuth. From these parameters, we can also compute the proxy for percent anisotropy as the difference between the fast and slow velocities normalized by the fast velocity. Rose diagrams from quantitative fracture analysis using azimuthal seismic velocity volumes are compared with the published data from Appalachian black shale outcrops and subsurface fracture models proposed in various papers to validate the analysis derived from the subsurface data. Although this velocity anisotropy, measured as azimuthal variations in velocity, can reflect rock fabric or stress, we show evidence that the likely source of these anisotropies is the presence of systematic joints.

Published data and azimuthal seismic attributes show two primary joint sets, the J1 and a J2 sets, as well as neotectonic J3 joints that affect the Marcellus and other Devonian shales in the Appalachian Basin (Engelder et al., 2009). Evidence suggests that in organic-rich Devonian black shale intervals including the Marcellus, J1 and J2 joint sets formed in sediments at or near peak burial depth as a result of anomalous pressures generated during thermal maturation of organic matter (Lash et al., 2004). Although authors indicate that the east–northeast to west–southwest J1 joint set is generally restricted to the black shales, the younger north–northwest to south–southeast J2 joint set is described as being more likely to extend out of the black shales into overlying rock. The late cracking of oil to gas released much larger volumes of gas during the period when J2 joints were propagated. Although neither contemporaneous with nor genetically related to folding, J1 and J2 joints in black shales within the Finger Lakes area of New York strike approximately parallel and perpendicular, respectively, to the Alleghanian orogenic fold axes exposed there. Similar patterns can be seen in seismic anisotropy for the Analog 3D survey and for the Clearfield 3D survey, to the west, in Clearfield County, Pennsylvania.

Rose diagrams from azimuthal attributes are similar to those measured for joints in outcrop in southern New York (Fillmore Glen State Park) and in Central Pennsylvania. The J2 azimuths dominate in areas with higher gas estimated ultimate recovery (EUR) wells, which is attributed to more gas generation as well as longer joint length. High-EUR areas also have generally lower “interval Vfast velocities” and show evidence of J2 joints well above the top of the Marcellus in Hamilton Group gray shales. Areas with low-EUR wells have a more dominant J1 trend as well as other scattered azimuthal trends, along with higher anisotropy, and a higher standard deviation of anisotropy, which suggests greater heterogeneity. Areas showing rapid spatial variation in the predominant azimuth of anisotropy imply possible reservoir compartmentalization and heterogeneities that correspond with areas of lower production. These attributes offer a tool to high-grade drilling opportunities and improve production estimates for Marcellus wells.

Geologic setting

The Middle Devonian Marcellus black shale was deposited on continental crust in an interior seaway of relatively shallow water (<200m=656ft), under relatively high sea level conditions (Kohl et al., 2014). During the Middle Devonian, a microcontinent called Avalonia depressed the southeast edge of the Laurentia continental margin (now the Appalachian Basin) as a consequence of thrust loading in a highland at the edge of the continent. This created the accommodation space, in which the Marcellus and subsequent black shales were deposited. The seabed was depressed below a pycnocline, which is a boundary below which oxygenated seawater does not circulate. Organic material, mainly from marine algae, was preserved in this oxygen-starved environment to be buried to a depth favorable for oil and gas generation. The Marcellus, which is the oldest of the units in the Hamilton Group, consists of two major cycles of organic-rich shale accumulation, with the basal portion being particularly rich. The Marcellus and other middle Devonian black shales in the basin are limited in extent due to the tectonically controlled variations in relative sea level change. One model for the development of the richest black shale is a very rapid transgression possibly accompanying thrust loading with disrupted river systems, reducing sediment flux into the basin, and favoring the accumulation of rock with a high total organic carbon (TOC) content. Eventually, river channels organized to deliver clastic sediments at a higher rate, thus diluting the organic content of the shale. Each of the two cycles of black shale, the Union Springs and the Oatka Creek members of the Marcellus Formation, is bounded above and below by a carbonate, indicating improved oxygenation/circulation of a basin still receiving minimal clastic influx.

The axis of regional folding in the study area is dominantly east–northeast to west–southwest (parallel to J1), primarily as a result of slip on a regional décollement beneath the Appalachian Plateau within the Silurian salt below the Marcellus (Davis and Engelder, 1985). Additionally, younger faults, especially north–northeast to south–southwest faults that appear to cut section above the Tully, a carbonate that caps the Hamilton Group, may allow gas to leak from the Marcellus reservoir. Both fault trends may affect production from the Marcellus. They may be significant when they are connected to fractures near the borehole and may act to rob the stimulation by carrying fluids away from the borehole. For this reason, most horizontal Marcellus wells are drilled to avoid the larger faults. In the Analog 3D survey area, wells have been drilled with the laterals north–northwest or south–southeast to be perpendicular to the contemporary stress field, and have terminated before reaching regional east–northeast to west–southwest-trending regional fault cuts; however, some wells appear to have been cut by the younger faults.

A seismic cross section oriented parallel to dip is shown in Figure 1, which displays the structural style of the project area within the Hamilton Group and adjacent formations. This line, also nearly parallel to most horizontal well trajectories, shows the variable thickness of the Silurian Salt below the Marcellus, which controls much of the larger scale structuring in this part of the Appalachian Basin. Thrust faults with a décollement in the salt, similar to faults shown in green, may be limited to the Marcellus and lower Hamilton group. Large faults, similar to the red faults, parallel to large folds related to salt deformation, are avoided when Marcellus wells are drilled. Young faults (not shown) are often parallel to dip. These faults have likely reactivated pre-existing zones of weakness and have variable offset at the Onondaga and Marcellus. Some have significant offset at and above the Tully. These larger structural features are well understood and well documented in the literature. The focus of this paper, however, is a study of reservoir level joints and fracturing using seismic attributes.

Even though this play is a “resource play,” this study shows that the reservoir fracturing is heterogeneous and that not all well locations are likely to be economic. The Marcellus play in northern Pennsylvania and southern New York is a “shale gas” unconventional resource play, commonly described as having no obvious trap or seal and no water contact. Unlike the younger units in the Hamilton Group, the Upper and Lower Marcellus Shale contain high gamma and high TOC facies, which generated gas just prior to and during the Alleghanian Orogeny. These shale plays depend on fracturing to allow mobilization of the gas during stimulation, as they have low, poorly connected matrix porosity.

Azimuthal anisotropy

Seismic anisotropy refers to variation in elastic wave propagation velocity that is directionally dependent. Conventional P-wave processing algorithms ignore azimuthally dependent NMO by assuming heterogeneous media, resulting in a single, isotropic velocity applied to all traces in a common midpoint (CMP) gather. As described by Tsvankin and Grechka (2011), applying a single value of NMO velocity to the whole gather in a wide-azimuth 3D survey causes underestimation of velocity in the Vfast direction and overestimation of the velocity in the Vslow direction for horizontally anisotropic media. In addition to improving the overall stack, azimuthally dependent velocity attributes can provide insight into the anisotropy of these data, including analysis of potential fracture trends, and reservoir heterogeneity. For the Analog 3D project, interval Vfast, interval VfastVslow percent, and interval Vfast azimuth volumes revealed significant variations in subsurface trends and heterogeneities that have been correlated to production.

Wide-azimuth anisotropic processing corrects the time shifts remaining in gathers that are initially migrated with the best isotropic velocity. The velocities are obtained by a surface fitting of the observed traveltimes in offset-azimuth space using the azimuthal traveltime-velocity equation given by Grechka et al. (1999). The azimuthal variation of NMO velocity is obtained by first picking traveltimes and then inverting those traveltimes for the azimuthal velocity. The traveltimes are obtained by locally windowed crosscorrelations in the CMP-time-offset domain as summarized by Jenner et al. (2001), Jenner (2011), and described in detail by Jenner (2001). The results are traveltime picks every 20 ms for every trace in the CMP gather, which are then inverted in offset-azimuth space using the azimuthal traveltime-velocity equation given by Grechka et al. (1999). Those residual time shifts that are not corrected by the isotropic velocity are inverted to create root mean square (rms) velocity volumes, rms Vfast, and rms Vslow. More delayed traces are those that traveled in the “slow” velocity direction, and less delayed traces traveled in the “fast” direction. The azimuthal position of the maximum and minimum (fast and slow) shows the direction of aligned geologic fabric (in the Marcellus, fracturing) as fast, if it is primarily aligned in a single direction. And rms Vfast and rms Vslow are defined as the maximum velocities of an ellipse, which defines the velocity variation with the source-receiver azimuth. The rms calculations are similar to an average value because they represent the sum of anisotropies from the surface. Azimuthal variations are commonly corrected in gathers where anisotropy is present to improve the stack; however, these data can also be used to identify anisotropies in the geology.

Figure 2 represents an example subsurface point, with high rms anisotropy. Consider aligned fracture systems or other aligned changes in sedimentary fabric to be “speed bumps” for a wave traveling perpendicular to the aligned system in an anisotropic medium. If fractures or stresses are aligned unidirectionally, the traveltime is slower for a direction that crosses the speed bumps, and would give a Vslow direction perpendicular to the fractures or alignments. Conversely, Vfast is parallel to these speed bumps. The difference between the calculated Vfast and Vslow as a percentage is a proxy for percent anisotropy. If there is low anisotropy relative to the scatter of velocity, the difference between Vfast and Vslow is small, such that the azimuth of Vfast cannot be computed accurately. As with sonic scanner data, these azimuth values are only useful if the percent anisotropy is calculated to be more than approximately 3%. If anisotropy is less than 3%, azimuthal values from these data are within the margin of error, and should not be used to develop a proxy for anisotropy. If the difference between Vfast and Vslow is large, the “sinusoid” representing velocity variation with azimuth can be calculated with confidence. Then, a proxy for anisotropy, defined by a VfastVslow percent greater than 3% can delineate the azimuth of Vfast with minimal error.

Interval velocities are calculated from the rms velocities using a generalized form of the Dix equation (Dix, 1955; Grechka et al., 1999). This process effectively strips off shallow layers, removing anisotropy from layers above the zone of interest and computing the anisotropy between the top and the base of an interval. Because a low dip is assumed in the calculation of the interval velocities, areas with higher time dip gradient were eliminated from any analyses of the data. Interval velocity calculations are based on a window crosscorrelation calculating the time shift required for a least-squares fit. This is calculated at every sample in a user-defined sliding window. The result may be noisy, so smoothing is often applied. Smoothed data may cause artifacts in the result, especially for azimuth calculations, so interpreters may consider using unsmoothed volumes for calculations, keeping in mind that some data values will be anomalous due to noise. The anomalous anisotropy values mostly occur on the edges of the survey and in low fold areas where the azimuthal contributions to the bins are not ideal. They also occur in fault zones where raypaths cross large discontinuities and where imaging may be poor. Because of this organization and predictability, we determined that unsmoothed volumes could still be used for analysis with production, where these anomalies can be avoided.

Natural hydraulic fracturing and the Marcellus play

Natural fracture systems are essential for achieving the best production in Marcellus Shale gas wells. Analysis of wide-azimuth, P-wave seismic velocity attributes at the reservoir level, and for specific laterals or proposed laterals, provides insight into these natural fracture systems in the subsurface. Published outcrop and core data, combined with azimuthal seismic attributes all show two primary joint sets, the J1 and J2 sets, as well as a neotectonic J3 set of fractures that affect the Marcellus and other Devonian shales in the Appalachian Basin. For the Devonian black shale intervals, J1 and J2 joint sets formed in sediments when they were at or near peak burial depth as a result of anomalous pressures during thermal maturation of organic matter (Lacazette and Engelder, 1992). The east–northeast-trending J1 joints are more closely spaced, and best developed in the more organic-rich black shale units (Lash et al., 2004). Although east–northeast joints (J1 joints) are not well developed outside the black shale intervals, joints (northwest-trending J2 joints) that formed during the Alleghanian orogeny are found throughout the Upper Devonian shale sequence. The two fracture sets are crosscutting within the shales, which is important for the optimization of well placement and fracking of horizontal wells. The earlier J1 fracture set (east–northeast trending), that resulted from initial gas generation, is nearly parallel to the maximum compressive normal stress of the contemporary tectonic stress field, a coincidence.

The J1 joint set appears to be unique to gas shales. The J2 set appears to break out of the gas shales and populate the rock above those gas shales. The second joint set may appear 305 m (1000 ft) or even as much as 1219 m (4000 ft) above the gas shale. We interpret this to mean that a large enough volume of gas was generated, so the section above the gas shale became overpressured to the extent it also was hydraulically fractured. Therefore, the section above the gas shale became charged with high-pressure gas as well.

It is hypothesized that the variable length of J2 joints and fractures that extend into the gray shales above the Marcellus could be related to TOC. Because the J2 joints are related to the volume of gas generation, we should see changes in the length and presence of J2 joints across the survey area if TOC is changing, and the amount of gas produced was variable across the survey. Figure 3 shows, diagrammatically, the relationship between J1 and J2 fractures. The J2 joints, which are Alleghanian in age and perpendicular to Alleghanian folds, grew episodically during the time of maximum gas generation. Although the J1 joint set is described as being more closely spaced within the black shale interval, the J2 joint set, which may occur over a much larger vertical section, may be more “visible” to the seismic tool because it is longer or because it is more likely to be mineralized. This is described in more detail in the analysis of rose diagrams from the azimuthal velocity volumes that follow.

Virtually all horizontal wells within the Analog 3D survey area to-date have been drilled in the north–northwest to south–southeast direction. In the ideal case, this allows horizontal wells to cross and drain J1 joints where they are present. Subsequent staged hydraulic fracture stimulations run east–northeast, parallel to J1, thus cross cutting and draining J2 joints. In spite of being perpendicular to stress, horizontal maximum (SHMax), J2 fractures are partially mineralized, propping open the fractures, and allowing them to contribute significantly to drainage of the Marcellus gas (Wilkens et al., 2014).

Black shales in outcrops within the Finger Lakes area of New York carry east–northeast J1 joints, and J2 joints approximately perpendicular to the Alleghanian fold axes exposed there (Lash et al., 2004) (Figure 4). Similarly, azimuthal analyses for the Clearfield seismic survey, west of the Analog 3D, in Clearfield County, Pennsylvania, show that J2 azimuths from subsurface seismic change orientation from north–northwest to west–northwest to remain perpendicular to the oroclinal fold. J1 azimuths in the Clearfield 3D survey are similar to those in the Analog 3D study area with an east–northeast trend similar, but perhaps slightly more east–west than the present-day stress direction SHMax. The systematic changing of J2 orientation from north–northwest to west–northwest shows that the azimuthal variations are not related to stress. Present-day stress has been shown to remain consistent across northern Pennsylvania (Engelder and Gold, 2008). Thus, unless the rock fabric shows linear spatial variation, the joints and fractures are the likely source of anisotropy measured by seismic azimuthal anisotropy.

Anisotropy and productivity from Marcellus wells

Three-dimensional azimuthal analysis gives us the subsurface information that we need to understand the potential fracture systems within the Marcellus, as well as within zones adjacent to the Marcellus. The results of these analyses show how azimuthal seismic velocity attributes compare with the published outcrop data describing joint and fracture orientations in and above the Devonian black shale. These data may also be used to better understand the reservoir heterogeneities that occur in the reservoir due to natural hydraulic fracturing and tectonic stresses. With limited well data and production data, we cannot completely evaluate all the variables that affect production from the Marcellus. We can, however, find excellent correlation between published data, in particular, azimuthal analysis of joints and fractures in outcrop, and azimuthal data from subsurface seismic. Additionally, there are numerous observations showing consistency with regard to azimuthal velocity attributes and EURs shown from decline curve analysis.

A map of northern Pennsylvania and southern New York (Figure 4) shows rose diagrams from joint analysis of shales in the Finger Lakes area of New York, along with rose diagrams calculated from the interval Vfast azimuth for the middle Marcellus, Cherry Valley, for the Analog 3D and the Clearfield 3D (all azimuths, where anisotropy is greater than 3%, as well as with low dip gradients). J1 joints are consistently east–northeast regardless of the fold axis. The J2 joints (Figure 4a) show an orientation perpendicular to the Alleghanian fold axes exposed in the Finger Lakes area (Lash et al., 2004). In these shales, the northwest–southeast and north–northwest to south–southeast J2 joints dominate. The fold axes continue south into Pennsylvania, and they show a similar trend for the Marcellus Shale. The rose diagram for the Cherry Valley interval from the Analog 3D survey (refer to rose diagram, Figure 4b, right) shows azimuthal trends similar to those on the eastern side of the Finger Lakes district. The Clearfield survey, a large survey that lies on the oroclinal fold belt for the Marcellus, shows a J2 azimuth (Figure 4b, left) that changes orientation from north–northwest for the northern part of the survey (Clearfield 3D north) to west–northwest in the central part of the survey (Clearfield 3D central), and then to nearly east–west in the southern part of the survey (Clearfield 3D south), staying perpendicular to structure.

Outcrops of the Appalachian Valley and Ridge contain some J3 joint sets that correlate with stress orientation diagrams from the World Stress Map (Hancock and Engelder, 1989; Figure 5). The Tully and Upper Hamilton Group, above the interval that is influenced by natural hydraulic fracturing, does not show J1 or J2 azimuths. The rose diagram calculated from the Top Hamilton azimuths in the Analog 3D, shown in blue in Figure 6, matches this contemporary stress field trend, which is typically between 58° and 69° east.

For further analysis of these trends at the scale of a single pad, we summarize what has been observed so far on the regional scale. First, the older J1 trend, attributed to natural hydraulic fracturing, is shown to be in the range between 70° and 77° (east–northeast) for outcrop and subsurface data in the Analog 3D survey and the Clearfield survey. The J2 trend, related to episodic natural hydraulic fracturing during maximum hydrocarbon generation, is not perpendicular to the J1 trend, but is perpendicular to Alleghanian folds, and generally trends north–northwest in the analog survey. The younger J3 trend is shown to be between 58° and 69° (northeast), slightly more north of east than the J1 trend in this area. The northeast azimuths captured in the shallower Tully interval on the Analog 3D survey match very well with World Stress Map (Heidbach et al., 2008) trends shown in Figure 5. The Tully is well above any likely influence from J1 and J2 natural hydraulic fractures or joints originating in the Marcellus Shale.

Initially, the azimuth of interval Vfast was extracted for three horizons (instantaneous amplitude of interval Vfast azimuth), the Top Hamilton, the Top Marcellus, and the Cherry Valley. Values for azimuthal analysis were limited to areas with greater than 3% anisotropy to insure a statistical calculation that is above the noise level for these data. The vertical variation in azimuth indicated by these velocity attributes match the expected results for the geologic model. The J2 azimuths, dominant in areas with higher EUR wells, are attributed to higher gas generation and longer joint length, as shown in our subsurface model (Figure 3). Areas with low EUR wells show a more dominant J1 trend, with limited influence from J2, perhaps because the joints are shorter where less gas is generated.

Gas estimated ultimate recovery from decline curve forecasts

Decline curve analysis is a reservoir engineering empirical technique that extrapolates trends in production data. The most commonly used trending equations are those first documented by Arps (1944), who was an American geologist that published mathematical relationships for the rate at which production from a single well declines over time. The gas EUR assumed an economic limit when the forecasted gas rate fell below 20Mscf/day. The economic limit is the point in time at which the production is assumed to cease because it is no longer economic. The EUR was determined by adding the cumulative historical production with the forecasted cumulative production to the economic limit. Publicly available data for historic well production for the Marcellus wells in this study were available from the state of Pennsylvania for six-month intervals through December 2013. To obtain a monthly production rate for decline curve forecasting, the six-month historic production data were divided equally among the preceding six months creating a stair-step profile. Decline curve forecasts were made assuming that hyperbolic decline converted to exponential decline when the instantaneous decline rate reached 10% per year. An example of a decline curve for a well in the Analog 3D survey is shown in Figure 6.

Estimated ultimate recovery comparisons with seismic attributes

The EUR values for wells in the project were used to create a bubble map with bubbles located at the midpoint of the horizontal well trajectory, and are corendered with Cherry Valley Vfast velocity, shown in Figure 7. No production tests or other production data were available to determine the variability of production along the length of the horizontal well; however, lateral lengths are similar throughout the study area, and all wells were drilled by the same operator. All of the horizontal wells were drilled after the seismic data were acquired; thus, none of the anisotropy is due to well stimulation. Clusters of data were analyzed for 16 small subset areas (white rectangles) of similar EUR, to correlate with seismic attributes. In map view, Figure 7 shows that the higher EUR pads are all within the lower velocity fairway (for amplitude extraction on the Top Cherry Valley horizon) outlined by the red oval. It is hypothesized that the lower interval Vfast velocities may be due to the presence of gas and/or pervasive fracturing. Vertical variations in anisotropy in these subset areas were also analyzed. The seismic cross sections A-A′ and B-B′ shown in Figure 7 for subset areas 1, 2, and 5 (Figures 8 and 9) show these variations.

For area 1, seismic cross section A-A′ (Figure 8), rose diagrams for the Onondaga, Cherry Valley, Marcellus and 25 ms (61m=200ft) above the Marcellus show lower anisotropy overall and an azimuth of Vfast in the J2 direction. Figure 7 shows that the interval Vfast velocity of this area is generally low. Area 1 has the highest decline gas EURs and may be showing a strong influence from J2 joints above the Marcellus, which have been attributed to more prolific gas generation, and higher TOC. From 50 ms above the Marcellus and shallower intervals, the azimuth of Vfast shows a northeast trend, which is in the range of azimuths attributed to younger J3 stresses or fractures.

Areas 2 and 5, shown in seismic cross section B-B′ (Figure 9), have low and high EURs, respectively. Area 2 has mixed high and low velocities (instantaneous amplitude extraction calculated from interval Vfast for the Cherry Valley), and it is quite close to two significant fault systems that may have negatively affected production. The rose diagrams calculated for the interval Vfast azimuth at the Cherry Valley in area 2 show dramatically different results than in area 1. The low-EUR area on the left side of the section shows a considerable amount of scatter in the azimuths with J1 or J3 dominant for the Marcellus and other intervals. There is a J2 component, but it is smaller than the J1 or J3 component. The higher EUR area on the right (area 5) shows a strong, nearly north–south, J2 trend for the Marcellus with minimal indication of the J1 trend.

Figure 10 shows rose diagrams from several high-EUR areas. These wells show a dominant J2 azimuth. This dominance of J2 suggests that these joints or fractures contribute significantly to production. In this analog area, where J2 is nearly perpendicular to the contemporary stress field SHMax, it may seem counterintuitive that J2 fractures would be open and contribute to production; however, closer examination of J2 fractures in core shows that these fractures are mineralized. Wilkens et al. (2014) discuss evidence including increases in background gas concentrations while drilling, indicating that these J2 fractures are more hydraulically conducive due to partial mineralization. Mineralization acts to prop open the joints or fractures and preserve porosity as well as permeability within the fracture. Additionally, they confirm that dipole sonic logs demonstrate that where these partially open veins are seen in core, azimuthal shear anisotropy spikes to values up to 10%. This observation supports our statement that the J2 azimuths, perhaps due to higher anisotropy from partial mineralization, dominate in areas with higher EUR. Figure 11 shows rose diagrams from two low-EUR areas. These rose diagrams show numerous azimuthal trends and a high degree of heterogeneity that is not seen in higher EUR areas. Although some of the areas show a smaller (lower count) J2 trend, it is always subordinate to other trends including the J1 and J3 trend azimuths. A greater degree of heterogeneity also appears to be a consistent characteristic of lower EUR areas, which have smaller areas with a consistent azimuth, and may also be adjacent to larger fault systems.

Other azimuths at N30E appear to be related to a young north–northeast to south–southwest-trending fault system, which dominates the eastern part of the survey. Areas with high azimuthal gradient (a high rate of change in azimuth) and a high calculated standard deviation of anisotropy might also be considered areas that are more heterogeneous. There is a strong relationship between higher EUR and lower anisotropy, as well as low heterogeneity, which is represented by the standard deviation of anisotropy across a subset area. In general, the highest EUR areas have an average anisotropy of less than 5%, with a standard deviation of less than 1.5% in the Analog 3D area.


In Devonian shales, a strong relationship between azimuths derived from seismic anisotropy attributes and those azimuths calculated from joints in outcrop suggest that these seismic attributes are a valuable tool for subsurface analysis of joints and fractures. Furthermore, variability in seismic anisotropy, heterogeneity, azimuth, and velocity along individual well trajectories indicates that the relationship between these attributes and production are predictable despite their complexity. For the analog area, we have shown that high EUR wells occur when a predominance of north–northwest (i.e., J2) azimuths are indicated by seismic anisotropy, along with lower velocity, lower heterogeneity, and lower anisotropy. These attributes can be used as a tool for high-grade Marcellus drilling and for designing better completions and stimulations. Potentially lower EUR areas, with more seismic heterogeneity, higher anisotropy, and mixed azimuths, may be scheduled for drilling after better locations are completed. Areas with dominant J2 azimuths that are oriented more east–west, for example, in the southern part of the Clearfield 3D survey (Figure 4), where J1 and J2 joints are more oblique, may point to necessary adjustments in well trajectories for optimal drainage.

The area of our Analog 3D seismic survey is of such a limited size (approximately 49km2 is equal to approximately 30mi2) that we presume that thermal maturation was uniform within the Marcellus and did not cause pockets of higher gas content. Certainly, the variation of fracture complexity is spaced at a much smaller interval and in a much more complex pattern than variations in either depositional patterns or thermal trends permit. These observations lead to the conclusion that complex patterns in productivity have another cause, maybe local variation in TOC. However, the correlation between more complex fracturing as indicated by a heterogeneity in azimuthal trends for anisotropy and low EUR suggest that over the period because maturation during the Permian there has been slow leakage of gas on a local scale and that the presence of a complex joint pattern was responsible. A single systematic joint set (i.e., J2) is not interconnected, so the bulk permeability of the immediate Marcellus approaches matrix permeability (Hubbert, 1957). If these joints are partially mineralized and have a very high permeability, they still cannot leak unless interconnected (Pommer et al., 2013).

The most direct evidence that the Marcellus has leaked over geologic time is the presence of mineralization of some, but not all J2 joints (Evans, 1995; Evans et al., 2012). Unfilled and mineralized J2 joints are found side by side in outcrops and core (Evans, 1994; Engelder et al., 2009). By the time of propagation of J2 joints, unfilled and mineralized, dewatering by compaction, and maturation had reduced the water saturation of the Marcellus to an irreducible state (Lash et al., 2004). The coexistence of unfilled and mineralized joints indicates that capillary forces in a gas-charged section prevented pervasive penetration of formation water when it invaded from deeper, more porous, and permeable beds such as the Oriskany Sandstone. Invasion of water from below is indicated by the presence of some mineralized J2 joints. It is likely that the combination of unfilled and mineralized J2 joints provides for the prolific gas production seen in the Marcellus.

Slow leakage through a complex system of multiple crosscutting fractures does not necessarily signal a rock volume that will rapidly leak stimulation fluid as implied by some in the literature (Warner et al., 2012). For example, the J2 joints are more likely to have grown out of the Marcellus and into overlying rocks (Engelder et al., 2009). Yet, these are the joints that seemed to have been “sealed” most effectively since the early generation. Likewise, there is no indication that the less economic wells in the study area have completely leaked as indicated by a robust flow-back drive by high-pressure gas. This means that a relatively effective capillary seal was maintained within the vicinity of the Marcellus since the Permian despite a complex fracture pattern (Cathles, 2001). This is true even for those wells relatively close to the zone of cross-formation faulting (Figure 1). Perhaps the effective capillary seal might have formed as high in the section as the Tully Limestone; however, there is no indication of a complete pressure breach through the Tully, even where there is a major disruption of bedding on the seismic scale.

Clearly, P-wave velocity attributes provide a critical understanding of Marcellus fractures and veins in the northeast and north-central parts of Pennsylvania, where the Hamilton Group is relatively thick. These velocity attributes may require additional evaluation for Marcellus fracture identification in western and southwestern Pennsylvania or West Virginia as the Marcellus, and the Hamilton Group in general, thins. It is worth considering methodologies to improve the resolution of these analyses by sampling velocity at smaller intervals, using smaller crosscorrelation windows, and by using horizon-based azimuthal velocity analysis. In areas with significant structure, velocity attribute analysis errors could be minimized by using prestack depth-migrated gathers as input to the analyses.

The Analog 3D study shows seismic attributes that can be quantitatively analyzed to evaluate undrilled areas of interest in the Marcellus play. Table 1 shows a summary, sorted by decreasing EUR, for the Analog 3D area, comparing interval VfastVslow percent as a proxy for anisotropy, the standard deviation of anisotropy as a proxy for heterogeneity, and the dominant azimuthal component relative to EUR for each subset area within the Analog 3D survey. Areas with high azimuthal gradient (a high rate of change in azimuth) and a high calculated standard deviation of anisotropy are considered to be areas that are more heterogeneous. There is a strong relationship between higher EUR and lower anisotropy, as well as low heterogeneity, which is represented by the standard deviation of anisotropy across a subset area. In general, the highest EUR areas have an average anisotropy of less than 5%, with a standard deviation of less than 1.5% in the Analog 3D area. Using these criterion, risk is reduced in exploration plays in the northern Marcellus play. These attributes can be analyzed to prioritize a drilling schedule, and to avoid drilling potentially uneconomic wells.


The Analog 3D study shows seismic attributes that can be analyzed to evaluate undrilled areas of interest in the Marcellus play.

In this study, we show that

  1. Regional fracture trends inferred from seismic azimuths correlate with published joint/fracture trends measured in outcrop, in other subsurface data and with World Stress Map trends.

  2. There are vertical and lateral variations in the following:

    • azimuth,

    • interval Vfast velocity,

    • implied anisotropy,

    • heterogeneity.

  3. Interpreted fracture trends differ between areas with larger decline EUR and smaller decline EUR values, making subsurface anisotropy analysis a predictive tool.

    • Higher EUR is indicated in subset areas where:

      • There is low anisotropy and low heterogeneity of anisotropy.

      • Velocity anisotropy is less than an average of 5%.

      • The standard deviation of anisotropy is less than 1.5%.

      • The azimuth of Vfast is dominantly in the J2 direction.

  4. Some perforations are likely to perform much better than others along the borehole, based on observed heterogeneity in the seismic profiles and map view.

  5. Reservoir characterization described in the literature for the fracturing or joints induced by gas generation, specifically the J2 trend are supported by these analyses. J2 fractures that break into the gray shales above the Marcellus and other Devonian black shales may give clues to the volume of gas generated, and thus to the TOC. It has been shown that J2 azimuthal trends, which have been attributed to these joint/fracture trends persist above the Marcellus in areas that have higher EUR.

  6. Some fault and fracture trends appear to be related to recent fault movement, and may adversely affect production.


We thank B. W. Horn and S. Volteranni at ION Geoventures for their support of this paper, and ION GXT for permission to present and publish this paper. Also, thanks go to M. Wallace, ION GXT, and K. J. McDonough for their help on various technical issues. Finally, the authors thank IS Interpretation Services, Inc., the Department of Geosciences, Pennsylvania State University, and Solutions Engineering for permission to publish this paper.

Tanya Inks received B.S. and M.S. degrees in geophysical engineering from the Colorado School of Mines. Since 1993, she has been providing integrated geophysical and geologic consulting services in Denver, Colorado. She is currently involved with projects in several unconventional plays including the Marcellus Shale play and the Niobrara play. She worked as a processor for Geophysical Service Inc., and CGG prior to graduate school, and worked for Mobil Exploration and Producing US Inc., in Rocky Mountain and Gulf of Mexico Exploration for the first six years following graduate school. Since 1993, she has consulted for many clients, initially as the manager of Vector Interpretation Services and later (1998) as a partner in IS Interpretation Services Inc. In addition to her work in the Marcellus, she has contributed geoscientific expertise to exploration and field development projects in structurally and stratigraphically complex areas, such as the Bearpaw uplift and Disturbed Belt in Montana, The Greater Green River, Wind River and Big Horn Basins of Wyoming, Utah’s Uinta Basin, the North Slope of Alaska, California’s Sacramento and San Joaquin Basins, Oklahoma’s Arkoma Basin, as well as international projects in Columbia’s thrust belt, Chile’s Fell Block and Venezuela’s thrust belt. She is a longtime member of the SEG, AAPG, DGS, and RMAG.


Terry Engelder received a B.S. (1968) from Penn State, an M.S. (1972) from Yale, and a Ph.D. (1973) from Texas A&M. He is a leading authority on the recent Marcellus gas shale play. He is currently a professor of Geosciences at Penn State and has previously served on the staffs of the US Geological Survey, Texaco, and Columbia University. Short-term academic appointments include those of visiting professor at Graz University in Austria and at the University of Perugia in Italy. Other academic distinctions include a Fulbright Senior Fellowship in Australia, Penn State’s Wilson Distinguished Teaching Award, membership in a US earth science delegation to visit the Soviet Union immediately following Nixon-Brezhnev détente, and the singular honor of helping W. Alvarez to collect the samples that led to the famous theory for dinosaur extinction by large meteorite impact. He has written 160 research papers, many focused on Appalachia, and the research monograph Stress Regimes in the Lithosphere. In the international arena, he has worked on exploration and production problems with companies including Saudi Aramco, Royal Dutch Shell, Total, Agip, and Petrobras. In 2011, he was named to the Foreign Policy Magazine’s list of Top 100 Global Thinkers.


Edward Jenner received a B.S. in physics with astrophysics from the University of Birmingham, UK, and an M.S. in geophysics from the University of Leeds, UK. Since 2002, he has been the manager of ION’s Land Research and Development team in Denver, Colorado. After completing a Ph.D. at the Colorado School of Mines in 2001, he joined AXIS Geophysics and developed techniques for azimuthal velocity and amplitude variation with offset (AVO) analyses. Since then, he has continued to focus on anisotropy, AVO, and anisotropic imaging issues for land seismic data. In 2003, he was awarded the SEG Clarence Karcher Award for his work in the field of azimuthal anisotropy.


Bruce Golob received a B.S. in geophysics from Bowling Green State University and an M.S. in information technology (dual major in database design and object-oriented programming) from Regis University. He is a geophysicist with ION-GXT in Denver, Colorado, and he works as a seismic inversion specialist, well integration, and processing quality control resource. He also assists clients in interpreting azimuthal anisotropy volumes produced from wide-azimuth seismic. Before joining ION-GXT, he worked in Denver interpreting 3D seismic for an independent oil and gas company in the U.S. midcontinent. He was originally hired and trained by Amoco — working for 18 years within a team of geoscientists and engineers, generating economic drilling prospects worldwide, his last five years with the company working in Cairo, Egypt as senior consulting geophysicist. He is currently a member of SEG, AAPG, DGS, and RMAG.


Jacki Hocum received a B.S. in mathematics from Arkansas State University, and after a brief stint as a high school teacher, she began her geophysics career in 1981 at Western Geophysical, where she processed on- and offshore and transition zone data from around the world, including North America, South America, the Gulf Coast, Alaska, and North Africa. She is a senior processing geophysicist at ION GXT in Oklahoma City, and she has more than 30 years of experience working in the seismic industry. As a group leader with Western, she was instrumental in training new processors. Before joining GXT in 2009, she spent 10 years with ECHO Geophysical, where her role as a seismic processor has expanded to include customer support and technical marketing. Since becoming part of the GXT team, she has concentrated on North American processing, with an emphasis on Niobrara and Marcellus plays. She is a member of DGS and GSOKC.


Darien G. O’Brien received a B.S. in petroleum engineering from the Colorado School of Mines and an MBA from the University of Alaska. He is the director of engineering with Solutions Engineering in Lakewood, Colorado. His technical expertise is in domestic and international reservoir engineering, drilling and work-over assessments, reserves and economic evaluations, reservoir modeling, environmental issues, and technology development. He has evaluated numerous oil and gas opportunities throughout the country, ultimately making significant acquisition and divestiture recommendations. He has prepared independent proved, probable, possible reserve reports for banks and investment groups, and he worked with the U.S. Securities and Exchange Commission to obtain sanction for proved resources in emerging exploration and resource development plays. He has been an active SPE member throughout his career, serving as continuing education chairman for the Denver Section, SPE distinguished lecturer, chairman of SPE’s International Continuing Education Committee, and recipient of SPE’s Regional Service Award. He is a member of the Petroleum Technology Transfer Council’s National Board of Directors, Society of Petroleum Evaluation Engineers, and the National Society of Professional Engineers. He is a registered professional engineer in Colorado, Wyoming, and California.


Freely available online through the SEG open-access option.