In unconventional reservoirs characterized by low matrix permeability, natural fractures (abbreviated as “fractures” for the rest of this chapter), even when sealed with calcite cement, can aid in the propagation of induced fractures during hydraulic fracturing treatment by acting as planes of weakness (Babcock, 1978; Fisher et al., 2002; Gale et al., 2007, 2014), but may also pose a hindrance to production because of resulting fluid loss and impeding the growth of new hydraulic fractures during hydraulic fracturing treatment (e.g., Gale et al., 2014). Placement of deviated wells perpendicular to the dominant orientation of the naturally fractured system is a standard methodology, which helps to intersect the maximum possible number of fractures (Gale et al., 2007; Sonnenberg et al., 2011). As such, a detailed understanding of natural fracture distribution and the relevant controlling factors may enhance the understanding and prediction of production anomalies that could be associated with the reactivation of natural fractures during or following the hydraulic fracturing treatment (e.g., Gale et al., 2007, 2014). Efforts directed toward the characterization and prediction of natural fracture systems on a regional scale have been pervasive and are commonly focused on seismic and structural attributes (e.g., Padgett and Nester, 1991; Narr, 1996; Ericsson et al., 1998; Bafia and Spencer, 1999; Pérez et al., 1999; Bai and Pollard, 2000; Gale et al., 2004; Kelley and Jones, 2013; Holman, 2014; Grossi, 2015) and petrophysical data (Bafia and Spencer, 1999; Dagistanova et al., 2011). Although stratigraphic and mineralogical aspects (e.g., lithology) are commonly recognized as a key factor in controlling fracture distribution (e.g., Ladeira and Price, 1981; Corbett et al., 1987; Helgeson and Aydin, 1991; Gross et al., 1995; Hanks et al., 1997; Dershowitz et al., 1998; Ericsson et al., 1998; Friedman et al. 1994; Underwood et al., 2003; Gale et al., 2004, 2007; Laubach et al., 2009, 2010; Sonnenberg et al., 2011) and can be the dominant control even over local structural deformation (Nelson and Serra, 1995; Hanks et al., 1997; Lorenz et al., 1997, 2002; Underwood et al., 2003), the lack of integration of fracture and mechanical stratigraphy into a high-resolution sequence stratigraphic framework may result in a limited understanding of the controlling factors of fracture distribution at the production scale, and a more limited application of relevant datasets in predicting subsurface fractured zones on an exploration or production scale. In addition, the impact of structural diagenesis on rock mechanical properties, from both a temporal and spatial perspective (e.g., Laubach et al., 2009), adds additional uncertainty when interpreting the mechanical stratigraphy during the time of formation of the natural fractures through the use of the present-day distribution of natural fractures.

In the unconventional “Mississippian limestone” play in the U.S. southern midcontinent, detailed core-based characterization of natural fractures has rarely been tied to a sequence stratigraphic framework. Correlating fracture distribution and intensity to a sequence stratigraphic framework may lead to the increased predictability of natural fractures in the subsurface, similar to which has been done in recent outcrop studies (e.g., Underwood et al., 2003; Cooke et al., 2006; Zahm and Hennings, 2009; Frost and Kerans, 2010; Zahm et al., 2010). The primary goal of this study is to characterize the type and distribution of natural fractures in the unconventional “Mississippian limestone,” and to tie key fracture parameters (e.g., type, intensity, spacing) with variations in mineralogy, facies types, and the established sequence stratigraphic framework to examine the controlling factors of natural fracture distribution, and to test the potential for prediction of naturally fractured intervals using a sequence stratigraphic approach.

The study area is located in north-central Oklahoma, part of the U.S. southern midcontinent (Figure 1A). During the Mississippian Epoch, this area was periodically covered by subtropical epeiric seas where a mixed carbonate–siliciclastic depositional system was developed. The general depositional system is interpreted to have been a distally steepened ramp (Childress and Grammer, 2019; Price and Grammer, 2019), with the study area spanning both proximal and distal depositional positions on the Cherokee platform (Figure 1B). Structurally, the study area spans the region east and south–southwest of the Nemaha fault zone (i.e., Nemaha uplift and Nemaha Ridge; Figure 1B, C), which is a key structural feature in north-central Oklahoma possibly active during the Late Mississippian to Early Pennsylvanian (Gay, 2003a, b).

An integrated approach was utilized to derive a comprehensive dataset on natural fracture attributes and distribution at various scales using five subsurface cores in north-central Oklahoma (Figure 1). The dataset was then interpreted within the context of a sequence stratigraphic framework (Figure 2). The surface of the cores serves as the primary data source that is accompanied by data from petrographic thin sections and micro-CT imaging. Fracture attributes such as height, kinematic aperture, spacing, termination style, and intensity were measured from the core surface. Because of the variable nature of kinematic aperture, which is defined as the accumulative opening of fracture, including opening space and cementation (Marrett et al., 1999), maximum values are estimated for subsequent analysis. To evaluate the fracture attributes beneath the core surface from a three-dimensional (3-D) perspective and at a finer scale, two 1.5 in. (3.8 cm) diameter core plugs were analyzed using micro-CT imaging technique at a resolution of 40 microns. Abundance of fractures is illustrated using fracture intensity, which refers to fracture count per meter of core (e.g., Ortega et al., 2006, 2010; Gale et al., 2014), and average fracture intensity is derived by dividing the fracture count by the corresponding footage of core of a selected interval (e.g., facies, depositional sequence). Mineralogical data from X-ray diffraction (XRD), such as calcite, quartz, and bulk clay content, add an additional facet to understanding the role of mineralogy in controlling the fracture distribution. These data were then compared with an established sequence stratigraphic framework to evaluate patterns that might assist with the prediction of fractures away from the well bore. The methodology utilized for fracture description and the definition of fracture attributes, as well as a means for distinguishing between induced fractures and natural fractures, are based on Kulander et al. (1990), Lorenz and Hill (1992), Nelson (2001), Lorenz et al. (2002), and Gale et al. (2007).

Seven lithofacies (five calcareous siltstones and two skeletal limestones) were defined based upon mineralogy, sedimentary structures, bioturbation intensity, and grain types (Figure 3), and interpreted to have been deposited on a distally steepened ramp. The lithofacies show repetitive vertical changes and apparent cyclicity that occur at various scales, which directly affect fracture distribution and reservoir compartmentalization. To present the rock data by highlighting potentially mineralogy controlled mechanical properties, five types of “petrophysically significant facies” (shortened as “P-Facies” for the rest of this chapter) are grouped based primarily upon the average XRD mineralogy and the extent of bioturbation (Figure 3). From P-Facies 1 to 5, general trends of increasing calcite content and decreasing clay content are observed, reflecting a generally shoaling-upward depositional trend (Figure 3). P-Facies 2, 3, and 4 exhibit similar mineralogical composition but were defined as separate facies petrophysically because of variations in the extent of bioturbation, which has been shown to possibly affect the resulting porosity and permeability values (e.g., Pemberton and Gingras, 2005; Gingras et al., 2012). The combination of the idealized vertical facies succession of depositional facies and the related succession of petrophysically significant facies (Figure 3) serves as the basis for constructing a hierarchical sequence stratigraphic framework, which is fundamental for integrating various datasets and testing the potential for increasing fracture predictability.

The glauconitic siltstone–fine sandstone facies is concentrated near the base of the Mississippian section and is volumetrically the most insignificant facies type. It varies from massive-bedded (Figure 4A, B) to locally laminated (Figure 4C), and is characterized by abundant fine sand-size glauconite grains (Figure 4D), which contributes to the distinct greenish color of this facies. Being subangular to subrounded and moderately sorted (Figure 4D), the glauconite grains possibly contain both transported and in situ components, the latter of which points to sediment starvation (Odin and Matter, 1981; Amorosi, 2012). Except for scarce Teichichnus (Figure 4B) and Zoophycos (Figure 4A) trace fossils, biogenic structures are rare. Similarly, the black-colored, massive-bedded siltstone facies rarely contains sedimentary or biogenic structures (Figure 4E) and is dominated by subangular to angular silt-size quartz grains (Figure 4F). With the highest average clay content (Figure 3), P-Facies 1 is interpreted as the facies with the lowest energy, being deposited in a distal outer ramp to basinal environment with generally restricted bottom waters.

The laminated siltstone facies is characterized by millimeter-thick, alternating clay-rich (darker colored) and calcite-rich (lighter colored) laminae (Figure 5A–C). The presence of millimeter-scale hummocky cross-stratification (HCS) in the calcite-rich laminae (Figure 5A) points to intermittent storm deposition as a potential mechanism of forming the calcite-rich layers (Cheel and Leckie, 1993). In this sense, variable proportions of the clay-rich and calcite-rich laminae (Figure 5A vs. 5C) suggest varying frequency and intensity of the storm events. The “rhythmic” appearance (Figure 5A, C) invokes the potential effects of tidal currents (Sarg, 2012), although the presence of HCS would likely negate that possibility. The laminated siltstone was likely deposited near storm wave base along the more proximal portion of the outer ramp relative to P-Facies 1.

The primary difference between P-Facies 3 and P-Facies 4 lies in the intensity of bioturbation: ranging from relatively localized and isolated burrow clusters in P-Facies 3 (Figure 6A–C) to those that are commonly connected and form a more homogenized burrow network in P-Facies 4 (Figure 7A–C). Trace fossil assemblages are characterized by a low diversity but high abundance of traces. The dominance of Phycosiphon (Figures 6A, 7A–C) suggests an environment with low energy and a fine-grained substrate (Goldring et al., 1991; Wetzel and Bromley, 1994; MacEachern and Burton, 2000; Bednarz and Mcilroy, 2009). Localized Zoophycos (Figure 6C) may suggest the occasional presence of coarser-grained substrates, possibly related to storm deposits (Pemberton and Frey, 1984; Pemberton et al., 1992; MacEachern and Burton, 2000). Teichichnus occasionally penetrates the storm-related, planar lamination (Figure 7D), illustrating their opportunistic nature in higher-energy environments. Additional evidence of fluctuating water conditions (e.g., oxygen levels and energy) includes variations in bioturbation intensity (Figure 7A) and the presence of carbonate-rich intervals (Figure 7B). Collectively, P-Facies 3 and P-Facies 4 are interpreted to have been deposited in a proximal outer ramp to distal middle ramp environment close to storm wave base.

P-Facies 5 contains the highest average calcite content (Figure 3) and is represented by two closely related subfacies. The massive-bedded facies (Figure 8A), which contains abundant peloids and crinoidal debris (Figure 8B), may be produced by storm- or earthquake-induced sediment liquefaction and transport by rapid sedimentation because of storms (Smith and Bustin, 1996; Boggs, 2006). For the hummocky cross-stratified, planar laminated packstone-grainstone facies, abundant HCS (Figure 8D) and planar laminations (Figure 8F) and scarce climbing ripples (Figure 8F) support an interpretation of rapid sedimentation during storms (Pemberton et al., 1992; Cheel and Leckie, 1993; Pemberton and MacEachern, 1997). Energy conditions during deposition fluctuated, as indicated by relatively low-energy mud drapes (Figure 8D, E) and burrowed zones (Figure 8C), both of which may represent quiescence following storms. High energy HCS beds in the massive-bedded facies (Figure 8A) points to the episodic nature of sedimentation (sensu Dott, 1983). P-Facies 5 is interpreted as a series of sand bodies that were deposited in the proximal outer ramp to distal middle ramp around storm wave base.

Based on morphological characteristics, four types of natural fractures were identified, including ptygmatic, vertical extension, shear, and zones with a mixed suite of fracture types (Figure 9). These fracture types are present in all of the cores examined in this study, indicating their likely presence throughout the “Mississippian limestone” play in the southern midcontinent. In terms of abundance, the ptygmatic fractures are the most abundant fracture type observed, and shear fractures are the least abundant (Figure 9). The fractures are primarily vertical to subvertical in orientation, and are generally sealed with calcite cement, although partially open vertical extension fractures occur locally.

Excluding the ones with a missing or broken fracture wall, the estimated kinematic aperture ranges up to 11 mm, with an average of 0.4 mm (Figure 10). In particular, kinematic aperture is variable from a 3-D perspective, as revealed in the micro-CT imaging (Figure 11). Including the fractures terminated at the core edge and because of the absent of core pieces, the measured height ranges from 1 mm to 710 mm, averaging 38.4 mm (Figure 12). Observed termination style of the fractures are categorized into tapering (at top and base) or abruptly terminating (at top or base) in a seemingly homogeneous portion of the rock, termination related to variations in mineralogy at top or base, and terminations at the core edge or because of missing core pieces at top or base. In general, tapering at both the top and base is the more common termination pattern, accounting for 53.8% of the fracture population (Figure 13). For fractures occurring in sets containing at least two individual fractures, which account for 19.5% of the fracture population, measured spacing ranges from 1 mm to 73 mm with an average of 18.2 mm (Figure 14), although the sampling bias related to the narrow core width (around 85 mm) should be considered.

The ptygmatic fractures are characterized by a folded morphology (Figure 15A). They occur as solitary fractures (Figure 15A), sets of parallel fractures (Figure 15B), and fracture bundles at the petrographic scale (Figure 15C), and are often discontinuous at a submillimeter scale (Figure 15D). When cutting through mud-rich laminae, the ptygmatic fractures commonly become increasingly contorted (Figure 16A). Brittle failure of the ptygmatic fractures (Figures 16B, 15D), variations in the direction of propagation (Figure 16C), and the abrupt occurrence of highly fractured intervals (Figure 16D) are locally present. In relation to bedding, the ptygmatic fractures terminate within both relatively thin (Figure 17A) and thick (Figure 16A) relatively rigid intervals and within relatively ductile siltstone facies where mineralogical variations are not distinct (Figure 15A, B). The fractures locally deform finer-grained laminae in both upward and downward directions, which can be observed at both core (Figure 17B) and petrographic scales (Figure 17C). However, most ptygmatic fractures commonly extend across calcite-rich layers (Figure 18A). In rare cases, densely arrayed short ptygmatic fractures span across silica-rich bands in P-Facies 5, some of which are seemingly related to the presence of silicified Zoophycos trace fossils (Figure 18B, C).

The vertical extension fractures are characterized by the absence of lateral offsets and straight fracture walls (Figure 19A). Locally, there are tall (e.g., over 80 mm in measured height) vertical extension fractures that, despite some calcite infill, are still partially open (Figure 19A). Without distinct evidence of dissolution at the contact between fracture-filling cement and fracture wall, the void space appears to be the result of a cessation of cementation. A variable relationship between the rates of calcite cement precipitation that occurred during and after fracture growth and variable rates of fracturing opening may be present (e.g., Gale et al., 2004, 2014; Olson et al., 2009).

Shear fractures (sensu Cooper and Lorenz, 2012) are characterized by a diagnostic “pinch-and-swell” structure (Figure 19B) that exhibits decreasing kinematic aperture at the inferred point where the shearing deformation occurred. In some cases, components of at least two of the three aforementioned fracture types coexist in a single fracture, which is herein defined as a “mixed” fracture type (Figure 19C). It should be noted that by “mixed type of fractures,” we do not intend to define a fracture as “mixed-mode” (opening and shear; Gale et al., 2014). Both the shear and mixed fracture types are relatively uncommon in the studied units (Figure 9).

The abundance of natural fractures is heterogeneously distributed and is controlled by the petrophysically significant facies types as well as the position within the sequence stratigraphic framework. In terms of facies control, P-Facies 5 (massive-bedded and HCS-planar laminated packstone–grainstone) exhibits the highest fracture count (Figure 20A) and average fracture intensity (Figure 20B) of all facies, which are correlated to the fact that P-Facies 5 has the highest average calcite content (Figure 3). This correlation between higher percentages of carbonate grains and more concentrated fracture distribution has been commonly documented and is thought to be associated with the higher strength because of the higher carbonate content or lower clay content in the more grain-supported textures in carbonate rocks (e.g., Corbett et al., 1987; Ericsson et al., 1988; Lorenz et al., 2002; Gale et al., 2007; Zahm et al., 2010). In addition, measured fracture spacing does not show a distinct separation among the P-Facies types (Figure 21).

The regressive phases of the interpreted third-order composite sequences generally exhibit a higher fracture count and average fracture intensity than the associated transgressive phases, particularly as P-Facies 5 commonly defines the regressive phases (Figure 22A, B). Fracture count and average intensity varies between sequences, invoking possible differences in the evolution of rock mechanical properties and intensity of structural deformation as external factors in addition to variable facies proportions, which are associated with fluctuations in the depositional system (e.g., water depth, oxygen level, and energy). Because the third-order sequences commonly show a clear upward-cleaning gamma-ray (GR) signature in this area, the correlation between fracture intensity and the sequence stratigraphy may enhance the predictability of fractures in the subsurface in nearby uncored wells, although clustering of the fractures should be considered when extrapolating core-based fracture datasets into well logs and inter-well spaces. Similar trends between the sequence stratigraphic framework and fracture distribution have also been observed in pure carbonate systems such as the Devonian reef outcrops in the Canning Basin in western Australia (Frost and Kerans, 2010). Similar to the observation in facies, measured fracture spacing shows similar values in the transgressive and regressive phases of the third-order sequences (Figure 23).

There are exceptions to the general pattern of fracture intensity being correlated to the third-order sequences, mostly seen when the third-order sequence may not be capped by packstone–grainstone facies, or when relatively abundant packstone–grainstone facies occur in the transgressive phases, most likely because of significant storm deposition. This illustrates the importance of rock data when calibrating well log data in predicting fracture distribution. At a finer scale, highly fractured intervals within the burrowed and bioturbated facies (Figure 15B, 16D) can further result in concentration of fractures in the transgressive phases. Similar patterns of higher fracture intensity in the “muddier” transgressive systems tract (TST) than the “grainier” highstand systems tract (HST) have also been documented in other units, for example, the Lower Cretaceous shallow-water limestones in south Texas (Zahm et al., 2010). Such fracture distribution is attributed to the lower unconfined strength and smaller bed or cycle thickness in the TST relative to those of the HST (Zahm et al., 2010). Variable evolution of rock mechanical properties may also play a role.

Mineralogical content, in particular, the calcite and clay content, directly controls yield strength and ductility and has been documented in several studies to affect fracture distribution (e.g., Corbett et al., 1987; Friedman et al., 1994; Gross, 1995; Rijken and Cooke 2001; Lorenz et al., 2002; Underwood et al., 2003; Lézin et al., 2009; Zahm et al., 2010). In this study, a general “first-order” mineralogical control of fracture distribution is observed. Average fracture intensity increases from P-Facies 1 to P-Facies 5 along with overall increasing calcite and decreasing clay content (Figure 24). Such a correlative relationship can be regarded as an average mineralogical representation of the facies as a control on fracture distribution. However, a detailed comparison between fracture count and whole-core XRD mineralogy (calcite and bulk clay content) is not successful. A key reason is the poorly defined sampling protocols (i.e., not aimed at specific facies or interval specific) that was used by the operator to sample once per unit interval or footage of core with variable sampling frequencies (from less than 1 m to more than 2 m [3–6.5 ft]). This nondirected sampling resulted in insufficient data points in thinly bedded intervals where facies variations frequently occur. Additional sampling bias is created when highly fractured intervals are skipped to filter out the potential impact of fracture-filling calcite cement on the mineralogy of the rock matrix. To test the relationship described earlier, high-frequency XRD data were collected from one of the cores in a 5 m- (16 ft-) thick “fourth-order” sequence to test the “high-resolution” (around 0.3 m [1 ft] per sample) mineralogical control on fracture distribution. This sequence contains all facies types except P-Facies 1 and exhibits well-constrained trend in facies-controlled fracture intensity. The results show a poorly constrained positive correlation between fracture intensity and calcite percentage, supporting the general premise that increasing calcite content leads to increased strength, even in the higher-frequency sequences or thinner units.

The highly folded ptygmatic fractures commonly occur in many of the unconventional reservoirs currently being worked in the continental U.S. (Gale et al., 2014), including the Barnett Shale of the Fort Worth Basin (e.g., Gale et al., 2007), the “Mississippian limestone” play in the southern midcontinent (this study), and the Bakken play in the Williston Basin (e.g., Sonnenberg et al., 2011). However, the mechanism of formation for ptygmatic fractures remains poorly understood. One possible scenario involves a critical condition when the rock is behaving as a ductile medium but still has the strength to break at a relatively early stage post-deposition (Figure 25). Ductile compaction of fractures is evidenced by the intense distortion along the fracture length (Figure 15A, B) and is possibly a product of the viscosity and resistant strength of both the mineralized fractures and the less competent host rock (Ramberg, 1959; Shelley, 1968). The increasing extent of distortion as the ptygmatic fractures cut through mud-rich laminae (Figure 16A), which can also be deformed (Figure 17B, C), supports this interpretation, pointing to the adaptive response of relatively ductile laminae to localized stress. In this sense, tortuosity of the ptygmatic fractures serves as a measurement of compressive strain (Ramberg, 1959; Shelley, 1968). Compaction at a relatively late stage is also likely present, as suggested by the brittle failure of the fracture (Figures 15D, 16B).

In contrast, the vertical extension fractures, which are characterized by the relatively straight fracture walls (Figure 19A), are inferred to be formed at both a relatively early (postdeposition) and late (postburial) stage as the rock obtains sufficient strength to break via tensile failure (Olson et al., 2007), reflecting a sense of displacement perpendicular to fracture wall and a pure opening mode (mode I; e.g., Olson et al., 2009). In the shear fractures, the “pinch-and-swell” structure (Figure 19B) may reflect mode II sliding with lateral shear stress being oblique or orthogonal to fracture wall (e.g., Olson et al., 2009). Although the dominant stress regime can be difficult to determine in the intervals of mixed fracture types, difference in rock mechanical properties (e.g., strength) at the time of fracture propagation can be evident where transformation of the fracture type occurs (Figure 19C).

Temporal variations of rock mechanical properties (i.e., structural diagenesis) result in the fact that the present-day fracture distribution may not reflect the rock mechanical properties when the fractures were initially formed (e.g., Gale et al., 2004; Shackleton et al., 2005; Olson et al., 2007; Laubach et al., 2009, 2010). Interpreted to have been formed respectively at a relatively early and an early (postdeposition) to late (postburial) stage, the ptygmatic and vertical extension fractures seem to favor different rock strength, which is likely further affected by the susceptibility of different facies to diagenesis. The primary porosity of the coarser-grained, higher-energy P-Facies 5 was preferably filled with extensive calcite cements following deposition (Figure 8B), which adds strength to the rock (e.g., Shackleton et al., 2005) and facilitates the formation of the vertical extension fractures at both a relatively early and late stage. Therefore, the vertical, cyclic occurrence of facies that exhibit various susceptibilities to petrophysically significant calcite cementation contributes to the compartmentalization of rock mechanical properties and overall fracture distribution (e.g., Gale et al., 2004). In this sense, the coexistence of the ptygmatic and vertical extension fractures in the same facies suggests the temporal evolution of rock mechanical properties related to the type and extent of diagenetic modification that results in the present-day fracture and mechanical stratigraphy observed in the play.

Absence of natural fractures in certain intervals adjacent to highly fractured zones, such as in P-Facies 5 (e.g., Figure 8A) which contains the highest average fracture intensity (Figure 24), points to the potential role of different burial conditions in controlling the spatial distribution of fractures (e.g., Laubach et al., 2009) and the uncertainties when tying fracture stratigraphy with sequence stratigraphy. Although such a temporal offset between the formation of ptygmatic fractures and present-day rock mechanical properties results in potential difficulties in revealing the rock and structural conditions at the time of fracturing and tying fracture stratigraphy with mechanical stratigraphy, pervasive diagenetic alteration of the mineralogical content in the rock matrix, such as dolomitization that can dominate over original mineralogy in controlling fracture intensity (e.g., Gale et al., 2004; Ortega et al., 2010), is not observed. This suggests that the distinct correlation between the fracture abundance and facies can be utilized to predict the fracture distribution using sequence stratigraphic approach in the subsurface Mississippian cores in this area. Although such a correlation could be overshadowed by postdepositional variations of rock mechanical properties (e.g., Olson et al., 2009), the grain texture, initial porosity, and mineralogical composition, all of which are affected by depositional environments and can be predicted by the relative positions in the sequence stratigraphic framework, are important to consider when addressing the issue of structural diagenesis when predicting and modeling fracture network.

In addition, structural diagenesis likely plays a key role in the clustering of fractures, which is difficult to evaluate when addressing the fracture dataset obtained from subsurface cores (e.g., Gale et al., 2014), and raises a challenge for upscaling the core-based data at a regional scale (e.g., Dershowitz et al., 1998). Among individual cores, the abundance of natural fractures is highly variable and can be partially explained by the various proportions of the petrophysically significant facies. Different potential for structural diagenesis (e.g., different spatial and temporal variations in rock mechanical properties; Laubach et al., 2009) affected by variations in facies distribution in different parts of the depositional system, as well as various activities relative to the structural elements in the area (e.g., fault distribution), is the other factor to consider. Although such sampling bias results in a challenge in predicting an exact fracture count in the subsurface, the correlative trend between fracture abundance and facies (Figure 20) and sequence stratigraphy (Figure 22) can provide insight for predicting the relative abundance of natural fractures in both cores and uncored wells of the “Mississippian limestone” play in this area through the utilization of sequence stratigraphy.

In addition to mineralogically controlled variations in the strength of the rock units, attributes of mechanical interfaces (e.g., thickness) can also be critical in controlling fracture termination, height, spacing, and intensity, especially in relatively massive-bedded stratigraphic units with subtle stiffness variations (Underwood et al., 2003; Corbett et al., 1987). A positive correlation between fracture height and bed thickness has been commonly recognized (e.g., Schultz and Fossen, 2002) and is also observed in this study, particularly in cases where short ptygmatic fractures are confined within relatively thin brittle layers (Figure 17A). Such a pattern of fracture confinement in layers with higher calcite content has been commonly documented in both cores and outcrops of fractured sedimentary rocks (Pitman and Sprunt, 1986; Helgeson and Aydin, 1991; Gross et al., 1995; Frost and Kerans, 2010; Zahm et al., 2010) and in unconventional reservoirs (e.g., Corbett et al., 1987; Friedman et al., 1994; Cooke et al., 2006), and is thought to be related to the internal deformation of ductile layers and the local opening or sliding of weak interfaces (e.g., bedding planes), the latter of which reduces the stress singularity at the fracture tip and results in subcritical crack growth (e.g., Olson et al., 2009). Consequently, the fracture tip loses propagation impetus because of dissipation of stress concentration, and the fractures are then restricted to the brittle layers, which possess lower tensile strength than the ductile layers (Corbett et al., 1987; Helgeson and Aydin, 1991; Cooke and Underwood, 2001; Renshaw et al., 2003; Cooke et al., 2006). More commonly, the step-over of fractures into adjacent mud-rich, relatively ductile intervals (e.g., Figure 18A) may be associated with the strongly bonded nature of the contact with a relatively high cohesion and friction coefficient, which promotes fracture propagation (Cooke and Underwood, 2001), and with the nonelastic deformation within the ductile interval (Rijken and Cooke, 2001). Based on laboratory tests, Friedman et al. (1994) proposed that the decreasing rate of fracture propagation and decreasing effective confining pressure, the latter of which may even transform the brittleness to ductility in mudrocks (Nygård et al., 2006), may also play a role. In cases where the ductile layers are relatively thin, stress at the fracture tip may not be dissipated by the internal deformation in the ductile layer so that fractures may continue to propagate into the next brittle layer (Figure 16A; Rijken and Cooke, 2001; Cooke et al., 2006), and waning propagation impetus at the fracture tip may be indicated by the decreasing kinematic aperture along the fracture length (Figure 16A). Although possibly affecting the average mineralogy of facies, the alternation of brittle and ductile layers can play a key role in determining the mechanical behavior of mudrocks (Gross and Engelder, 1995).

It has also been documented that thinner beds tend to contain more closely spaced short fractures than thicker beds with similar mechanical properties (e.g., McQuillan, 1973; Ladeira and Price, 1981; Corbett et al., 1987; Gross et al., 1995; Cooke et al., 2006; Sirat et al., 2007). In this study, closely spaced, relatively tall fractures (e.g., 15–30 cm [6-12 in.] in measured height) are present in thick (e.g., around 1 meter [3 ft]), massive-bedded intervals (Figure 16A), and singular, short fractures are common in thin calcite-rich layers (e.g., 1–5 cm in thickness; Figure 17A). However, an effort to delve into the relationship between bedding thickness and fracture attributes (e.g., height and spacing) was not successful, because of several reasons. The first and foremost is the difficulty to define beds in an unambiguous and definitive way, because of the common lacking of distinct variations in mineralogy (i.e., gradational contact) and common presence of thin laminae in a “bed.” Consequently, the thickness of “a bed” commonly equals an interval of one facies type that can be several meters (several feet) thick with gradational contacts and contains fractures with highly variable height and complex stacking patterns (e.g., the top tip of one fracture is commonly adjacent to the basal tip of another fracture; Figure 16A). Such loosely constrained control of bedding on fracture distribution is a major reason that fracture attributes are not reported along with fracture abundance in this study. As a result, the variable constraining relationship between fractures and bedding, when present, points to an overall poorly constrained relationship of bedding thickness to fracture height, further contributing to the uncertainties when relating bedding structures to mechanical layering. In addition, temporal evolution of rock mechanical properties and localized stress in intervals with similar mechanical properties may have further contributed to such a constraining relationship between bedding and fractures by producing a complex contrast of relative rigidity among units (e.g., Shackleton et al., 2005), which is illustrated by the deformed mud-rich laminae intersected by the ptygmatic fractures (Figure 17B, C). All of these factors may also be responsible for the tapering termination of fractures in zones seemingly without distinct mineralogical changes (Figure 15A, B) and the coexistence of various types of fractures in the same facies, illustrating the complexity in the attributes, distribution, and controlling factors of the fracture system, and the uncertainties when integrating fracture stratigraphy, mechanical stratigraphy, and sequence stratigraphy at a whole core scale. Timing of fracture development may also be a controlling factor, as “early” fractures tend to be independent of bedding thickness, and “late” fractures may be related to thickness of the brittle layers (e.g., Renshaw et al., 2003; Frost and Kerans, 2010).

Open fractures serve as key conduits for fluid migration, and therefore are crucial for reservoir permeability, especially when matrix porosity is extremely low (e.g., McQuillan, 1973). The localized presence of void space in the vertical extension fractures (Figure 19A) suggests that these type of fractures may be present as clusters in the subsurface similar to what has been inferred in the Barnett Shale and observed in the Austin Chalk (Gale et al., 2004; 2007). Therefore, these partially open vertical extension fractures may provide reservoir permeability at some scale, the quantitative extent of which is unknown. On the other hand, virtually, all of the ptygmatic and shear fractures are sealed with calcite cement and, therefore, likely to provide little, if any, primary contribution to reservoir performance. These mineralized fractures, however, may promote the propagation of induced fractures during hydraulic fracturing treatments by serving as planes of weakness and reactivation (e.g., Gale et al., 2007; Hu et al., 2014).

There are several potential limitations that must be considered in the application of the results reported in this study. The first and foremost relates to the nature of the core. The distance between the cores available for this study (ranging from several kilometers to tens of kilometers [several miles to tens of miles]) directly results in a sampling bias, which only reflects a clustered fracture distribution in the area. The narrow width of the core, which omits fractures with spacing wider than the core width, creates another sampling bias and potentially masks the “true” occurrence style and abundance of fractures (e.g., a singular fracture may be part of a fracture set not captured by the core; Figure 15A vs. B). All of these scenarios can lead to an incomplete picture of lateral fracture distribution, suggesting the scale-dependent nature of the fracture dataset. Gale et al. (2004) and Ortega et al. (2006, 2010) developed a scale-independent method to characterize fracture abundance with intensity and spacing data. Different from this study, their method was applied on outcrops and requires beds with a lateral extent that cannot be observed in core. In addition, average fracture spacing is commonly considered as the inverse of average fracture intensity in beds for closely spaced fractures (e.g., Ortega et al., 2006; Gale et al., 2007). Because of the highly variable and difficult nature of defining beds in much of these cores as discussed earlier, this trend was not observable in this study either. As such, the approach of counting fracture abundance from the perspective of unit length of core may be the best applicable method for the narrow subsurface cores, at least in this study.

The actual measurements of the fracture attributes may also be equivocal. For instance, true height of the fractures is underestimated when the fractures are observed terminating because of missing core pieces or at the core edge, which accounts for 23% of the fracture population (Figure 13). This effect can be worsened by an inclined instead of a vertical orientation of the core. Because of this potential sampling bias, measured height of the fractures is not compared with facies and sequence stratigraphic framework. The kinematic aperture data of this study is estimated using a hand held millimeter-scale ruler, possibly resulting in some level of data imprecision. Ortega et al. (2006) used a microscope-calibrated comparator for accurate measurement of kinematic aperture. However, considering the highly variable kinematic aperture exhibited by the sinusoidal ptygmatic fractures (Figure 15A, B) that commonly taper (Figure 13), this comparator has its inherent limitation for this study. Examination of the “true” kinematic aperture can also be affected by the occasional breakage of core along fractures and the angle of the core surface intersecting the fracture plane. In addition, fractures terminating beneath the core surface, which are common based on micro-CT imaging results (Figure 11A) and petrography (Figure 15D), cannot be captured in a core surface-based investigation. Along with the variable 3-D kinematic aperture (Figure 11B), these features, which have also been observed in micro-CT imaging results in the Cretaceous Eagle Ford Formation (Hu et al., 2014), illustrate the necessity of incorporating 3-D imaging techniques as well as petrography, both of which may help upscale the fracture data (e.g., Gale et al., 2004) and reinforce the comprehensiveness of a core surface-based fracture dataset. Therefore, integrating fracture data at various scales from well log, outcrop, core surface, petrographic thin-sections, and high-frequency micro-CT imaging would be most applicable to developing a comprehensive natural fracture dataset. Further integration of fracture data with high-frequency porosity–permeability and rock mechanical data can be used to develop a composite geomechanical–reservoir model (Gale et al., 2004; Sirat et al., 2007; Zahm et al., 2010).

The final major limitation involves the observation of the mineralogical control of the fracture distribution and intensity, which suggests that nonsystematic sampling for XRD analysis may fail to capture the smaller-scale variability in fracture intensity related to bedding. More selective and higher-frequency sampling would lead to a more reliable relationship between mineralogy and fracture distribution at the whole core scale.

This study illustrates the value of integrating a sequence stratigraphic approach into characterizing and predicting the distribution of natural fractures in the “Mississippian limestone” play in north-central Oklahoma. Four types of fractures are identified: ptygmatic fractures (most abundant), vertical extension fractures, shear fractures (least abundant), and zones of mixed types of fractures. The highly folded ptygmatic fractures are interpreted as forming relatively early via ductile buckling prior to lithification. As shown by micro-CT imaging, fractures are commonly discontinuous at the millimeter scale and terminate beneath the core surface, with a variable kinematic aperture from a 3-D perspective. The vast majority of the fractures are sealed with calcite cement with only local remnant void space observed in vertical extension fractures. These partially open vertical fractures likely contribute to reservoir permeability at some scale.

Results of this study indicate that the distribution of naturally fractured units in the “Mississippian limestone” play is heterogeneous and is controlled by both the facies types and the sequence stratigraphic framework. Among the five types of petrophysically significant facies, the massive-bedded and hummocky cross-stratified to planar laminated packstone–grainstone (P-Facies 5), which exhibits the highest calcite content, contains the most abundant fractures and the highest average fracture intensity, supporting the correlation between mineralogy and fracture distribution. From a sequence stratigraphic standpoint, the regressive phases of the third-order composite sequences exhibit higher fracture intensities than the associated transgressive phases when P-Facies 5 is the dominant facies. Because these third-order sequences can be identified by cleaning-upward GR patterns, the tie between fracture intensity and sequences that can be correlated sub-regionally provides a valuable tool to assist in the prediction of fractures in the subsurface away from cored wells.

This research was supported by the Oklahoma State University Industry Consortium on the Reservoir Distribution and Characterization of the Mid-Continent Mississippian Carbonates—A Major Unconventional Resource Play (American Energy Partners, Chaparral Energy, Chesapeake Energy, Devon Energy, Longfellow Energy, Marathon Oil, Maverick Brothers, Newfield Exploration, SM Energy, Samson Energy, Sinopec/Tiptop, Red Fork Energy, Trey Resources, and Unit Petroleum). We would like to thank the Oklahoma Geological Foundation for the Herbert G. and Shirley A. Davis Geology Fellowship and the American Association of Petroleum Geologists Foundation for the Grants-in-Aid award. Appreciation is also extended to the Powder XRD Core Facility at Department of Physics, Oklahoma State University, for assisting with the mineralogical analysis. The authors also gratefully acknowledge the valuable assistance from John C. Lorenz during the initial stages of this study and insightful comments from Julia F. W. Gale and one AAPG reviewer who helped significantly improve this manuscript.