Although chemostratigraphic methods have been used in petroleum exploration since the 1980s (Scholle and Arthur, 1980), recent studies have focused on evaluating the potential for characterizing unconventional hydrocarbon reservoirs at high-resolution scales (Rowe et al., 2012; Koch et al., 2014). Chemostratigraphy can be defined as the study of the temporal and spatial geochemical variability of the rock record with a goal to establish genetically related, mappable rock units, defined in terms of their unique geochemical composition (Ramkumar, 2015). The importance of establishing a genetic, and often an age-specific, relationship between rock units is a fundamental aspect of chemostratigraphy as it allows that unit to be correlated both vertically and laterally with enhanced confidence.

Although the midcontinent Mississippian system is known for its complex reservoir lithologies and vertical and lateral heterogeneities, recent studies (LeBlanc, 2014; Price, 2014; Childress and Grammer, 2015; Jaeckel, 2016; Vanden Berg, 2016; Childress and Grammer, 2019; Vanden Berg et al., 2019) have demonstrated the ability to predict the highly cyclic nature of the Mississippian units in Oklahoma at multiple scales (Fischer, 1964). A repeatable pattern comprised of shallowing-upward facies successions and a hierarchy of depositional cyclicity is recognized throughout the basin, enhancing the predictability of the depositional system overall from lateral facies distribution to vertical reservoir compartmentalization. Our findings suggest that chemostratigraphic correlations of carbon and oxygen stable isotope curves are related to the depositional hierarchy and can be applied at the basin scale. Overall, the high-resolution geochemistry has provided a revised interpretation of the facies-based hierarchy first defined by LeBlanc (2014) based upon a new set of cycle characteristics that incorporates diagenetic susceptibility, more frequent depositional cyclicity related to facies stacking patterns, and a potential global chronostratigraphic tie.

Carbon and oxygen stable isotopic data is a traditional method used for chronostratigraphic correlation and paleoceanographic studies on a global scale (Bruckschen et al., 1999; Mii et al., 1999; Saltzman, 2003; Batt et al., 2007; Buggisch et al., 2008; Saltzman and Thomas, 2012); however, recent studies have tested its potential for regional applications aimed at correlating stratigraphic sections in the midcontinent. Koch et al. (2014) reported three important shifts in δ¹³C in the Anadarko Basin of north-central Oklahoma, which they correlated with the lower Tournaisian (Kinderhookian), upper Tournaisian through middle Visean (upper Kinderhookian–Meramecian), and uppermost Visean (lower Chesterian) substages. Koch et al. (2014) reported similarities in trends between the Anadarko shelf profile and the global type section of the Mississippian–Pennsylvanian subperiod boundary at Arrow Canyon, Nevada, and concluded that δ¹³C values overall were not facies dependent and instead reflect global changes in seawater chemistry during the Mississippian.

On a regional scale, Koch et al. (2014) considered broad energy level differences in depositional environment and proposed widespread meteoric diagenesis as the source for excursions in their shallowest depositional facies (uppermost Mississippian–lowermost Chester Group), but did not relate the geochemical signatures within the context of a high-resolution sequence stratigraphic framework. Idealized facies successions and their resulting sequence stratigraphic hierarchy are fundamental in understanding sea level–driven cyclicity within a stratigraphic section, and ties to a sequence stratigraphic framework would lead to enhanced predictability away from the core being utilized in the study (Grammer et al., 2004). This predictability can then be applied to different orders of cyclicity and isotopic signatures can be correlated at multiple scales. This study’s comparison of isotopic signatures to a previously defined sequence stratigraphic framework (LeBlanc, 2014) is crucial for capturing both large and small-scale geochemical variability and evaluating the correlation potential of these signatures at multiple scales.

The vertical and lateral heterogeneity of Mississippian reservoirs has led to nomenclature inconsistencies that complicate correlation efforts and often cause the system to be defined lithostratigraphically. This lithostratigraphic approach often leads to the misidentification of units and a relatively simplistic view of reservoir distribution and internal communication within a complex depositional system. Recent midcontinent studies have shown that the Mississippian time interval represents a time-transgressive system that is characterized by repetitive facies belts that are deposited in time equivalent packages and migrate as a result of base level rise and fall (LeBlanc, 2014; Price, 2014; Childress and Grammer, 2015; Jaeckel, 2016; Vanden Berg, 2016). Multiple orders of eustatic sea level change during the Mississippian result in lateral discontinuity of reservoir units and in instances where biostratigraphic data and strong chronostratigraphic ties are lacking, stable-isotope chemostratigraphy may be a useful method for correlating depositional units in the subsurface and to facilitate accurate regional correlation.

The single core being evaluated in this study is interpreted to have been deposited on a distally steepened ramp of the Anadarko shelf that experienced multiple orders of sea level fluctuation resulting in facies that were deposited mostly in normal marine conditions with some facies deposited in more restricted, possibly oxygen deficient waters (LeBlanc, 2014). The Adkisson #1–33 (98 m [323 ft] long) is located in Logan County, Oklahoma, on the Cherokee Platform due east of the Nemaha uplift (Figure 1). The Adkisson #1–33 contains contacts with the underlying Devonian Woodford Shale and overlying Pennsylvanian shale and is interpreted to represent a full Mississippian section (Figure 2).

Overall, the Mississippian represents a transitional period from greenhouse conditions dominant during the Devonian to icehouse conditions, which were dominant during the Pennsylvanian and Permian (Read, 1995). In addition to climatic and oceanographic changes, a transition from greenhouse to icehouse times implies the enhanced effect of glacioeustatic sea-level change during which obliquity signal (40 k.y.) within the Milankovitch band is assumed dominant (Read, 1995). This signal, along with the combined effects of eccentricity (100 k.y. and 400 k.y.) and precessional signals (19–23 k.y.), is then superimposed on larger, third-order cycles (0.5–5 m.y. duration), which have been shown to form the smallest scale of depositional sequences that are most confidently correlated within recent midcontinent studies (LeBlanc, 2014; Price, 2014; Childress and Grammer, 2015; Jaeckel, 2016; Vanden Berg, 2016; Childress and Grammer, 2019; Vanden Berg et al., 2019).

Well defined, third-order composite sequences have been identified as global in scale and also reflect the greenhouse to icehouse transition (Figure 3). The Tournaisian–Visean (Kinderhookian, Osagean, and Meramecian) is characterized by relatively long third-order sequences (avg. 3.6 m.y., 2.3 m.y., and 2.0 m.y. respectively), which may reflect low-amplitude sea-level change (10 m [33 ft]) during the Early Mississippian (Read, 1995). Alternatively, late Visean–Serpukhovian strata (Chesterian) is characterized by relatively short, and more frequent third-order sequences (average 1.3 m.y.), which may reflect enhanced glacioeustatic effects inducing rapid and high-amplitude (75–100 m [246–328 ft]) sea-level change during the Late Mississippian (Read, 1995). This has important implications for the prediction of depositional cyclicity as slight changes in the pattern of shallowing-upward facies successions, the resulting depositional hierarchy, and the potential for subaerial exposure should be expected from Early to Late Mississippian.

Bulk sediment samples were drilled using a microdrill from the backside of core slabs approximately every 0.3 meters (1 ft) and analyzed for stable isotopes of carbon (δ¹³C) and oxygen (δ¹⁸O) A clean surface was prepared for each sample by drilling away surface residues to prevent contamination from drilling mud or other chemicals associated with the drilling and coring process.

High-resolution sampling allows for the identification of small-scale heterogeneities and is similar to Koch et al. (2014) wherein bulk sediment samples were drilled approximately every 0.48 meters (1.6 ft). Stable isotope analysis was conducted at the University of Miami (RSMAS) Stable Isotope Laboratory using a common acid bath interfaced to a Finnigan-MAT 251 mass spectrometer. Data are reported using the conventional per mil notation (‰) relative to Vienna Pee Dee Belemnite (VPDB). The precision for this method is ±0.1‰.

In addition to high-resolution analysis, high-density sampling across third-order sequence boundaries defined by LeBlanc (2014) was conducted. Five 5-sample transects (0.9–2.0 cm [0.35–0.8 in]) were measured across each third-order sequence boundary, including the basal contact with the Devonian Woodford Shale and the upper contact with the Pennsylvanian shale.

Recent work suggests that midcontinent Mississippian facies and their vertical successions best represent deposition on a distally steepened ramp (Price, 2014; LeBlanc, 2014; Childress and Grammer, 2015, 2019; Jaeckel, 2016; Vanden Berg, 2016; Vanden Berg et al., 2019). The predominant facies types present in the Adkisson #1–33, as defined by LeBlanc (2014), range from fine-grained, glauconitic sandstones to crinoidal-brachiopod packstones and grainstones (Figure 4) and have been interpreted to reflect a shallowing-upward succession in close proximity to fair weather wave base in a relatively distal position on the ramp. Changes in bioturbation indices and diversification of faunal assemblages support the hypothesis that facies present at the base of the succession were likely deposited in low-energy, potentially restricted conditions and facies present toward the top of the succession were deposited in normal marine conditions with improved circulation. Each facies is briefly summarized herein; for full facies characterization, see LeBlanc (2014).

Facies 1 is a fine-grained, glauconitic sandstone characterized by well-rounded to subrounded, well-sorted, silt-size quartz grains (30–60%) and rounded to subrounded, moderately sorted, sand-size glauconite grains (40–50%; Figure 4). Although glauconite can form in various marine environments, it is commonly interpreted to indicate deposition within an anoxic, submarine environment extending over a large area characterized by normal salinities, low energy conditions, and low sedimentation rates (Middleton et al., 2003). Facies 1 is only present in the lowermost sections of the core and is interpreted to represent deposition in a restricted, low-energy environment with limited circulation and reflect the initial stages of flooding on a regional scale.

Facies 2 is a burrowed, calcareous siltstone characterized by weak calcareous laminations and localized millimeter-scale burrows (Figure 4). Limited abundance and low diversity of skeletal grains (siliceous sponge spicules and thin-shelled brachiopods) indicates that conditions were not suitable for a diverse, normal marine fauna (Flügel, 2010). Alternations in the degree of preservation of laminations due to millimeter-scale burrowing suggest fluctuating environmental conditions between oxygen-poor and oxygen-rich conditions (Ekdale et al., 1984). Facies 2 is interpreted to represent outer to distal ramp deposition below fair weather wave base in a low-energy, restricted environment that represents intermediate conditions between the underlying glauconitic sandstone facies and the overlying bioturbated, calcareous siltstone facies.

Facies 3 is a bioturbated, calcareous siltstone characterized by abundant fine sand-size peloids (30–40%), silt-size quartz grains (average 40%), a moderate bioturbation index (2–4 BI), and local, centimeter-scale vertical burrows (Figure 4). The presence of more diverse skeletal material (brachiopod and crinoid fragments), overall moderate bioturbation index, and presence of vertical burrows suggest well-circulated, normal marine conditions during deposition. Facies 3 is interpreted to represent mid- to outer-ramp sediments deposited under low to moderate energy settings.

Facies 4 is a peloidal packstone to grainstone characterized by abundant fine sand-size peloids (50%), silt- to fine sand-size quartz grains (35%), and skeletal debris (15%) composed of crinoid, brachiopod, and lesser amounts of bryozoan fragments (Figure 4). The presence of diverse and more abundant skeletal material suggests well-circulated, normal marine conditions during deposition. Facies 4 is interpreted to represent deposition within the mid-ramp or distal portion of the ramp crest environment near fair weather wave base, proximal to skeletal shoals.

Facies 5 is a skeletal packstone to grainstone characterized by disarticulated, very fine- to medium sand-size brachiopod, crinoid, bryozoan, and siliceous sponge spicule skeletal debris (Figure 4). Cross-bedding, local hummocky and swaley cross stratification, and inclined bedding planes are commonly observed in core. Intermittent muddy intervals, burrowed firmgrounds, bioturbated intervals, and shell beds are also occasionally observed. The presence of cross bedding and truncation surfaces suggests that Facies 5 may have been deposited in an active portion of a skeletal shoal; however, intermittent mud-rich intervals suggest a more likely interpretation of a distal expression of a skeletal shoal complex. Hummocky cross-stratification indicates that deposition was occasionally influenced by storm processes (Middleton et al., 2003), which has likely caused autocyclic deposits in some vertical successions. Facies 5 is distinguished from Facies 4 by its relative larger grain size and the presence of high-energy bedforms and is interpreted to represent deposition within the mid-ramp or distal portion of the ramp crest environment.

By utilizing an idealized stacking pattern, LeBlanc constructed a sequence stratigraphic framework composed of three orders of cyclicity from third-order sequences to fifth-order high-frequency cycles (LeBlanc, 2014). An alternative interpretation composed of higher orders of cyclicity has been proposed for the same framework consisting of one second-order sequence, four third-order sequences, and 18 fourth-order high-frequency cycles (Figure 2). Due to a lack of biostratigraphic data for the Adkisson #1–33 core, the hierarchy only reflects relative cyclicity without a means to temporally constrain the periodicity. This alternative first-order interpretation served as the framework for isotope sampling and provides a fundamental tie back to sequence stratigraphy allowing for a true chemostratigraphic approach to be applied.

Three hundred and sixty-four bulk sample powders yielded δ¹³C and δ¹⁸O values from −4.54‰ to +3.53‰ and −8.28‰ to −0.41‰, respectively. These values are somewhat consistent with Mississippian seawater values from Mii et al. (1999); however, many δ¹⁸O values and some δ¹³C values plot outside of this range in relatively depleted fields (δ¹⁸O: −4.0‰ to −6.25‰ and δ¹³C: +1.0‰ to −4.5‰; Figure 5).

Incorporating high-resolution stable isotopic data with gamma ray signatures and facies relationships resulted in a new, chemostratigraphically driven interpretation of third-order cyclicity for the Adkisson #1–33 (Figure 6). The criteria for choosing new third-order sequence boundaries are high gamma ray values (> 50 API), shallowing-upward facies successions, and depletions <1.5 meters (5 ft) below and enrichments <1.5 meters (5 ft) above newly defined boundaries in δ¹³C and δ¹⁸O values. Thirteen third-order sequences have been identified using these criteria and suggest more frequent third-order cyclicity than previously interpreted (Figure 2). If all four North American stages are present, the reinterpreted section also suggests a loose correlation with globally recognized third-order sequences from Haq and Schutter (2008; Figure 3). Of the 13 newly identified third-order sequences, 2 are hypothesized as Kinderhookian, 3 as Osagean, 4 as Meramecian, and 4 as Chesterian. Deviations from this pattern (e.g., abbreviated Kinderhookian, Osagean, and Chesterian strata) may be due to nondeposition or erosion associated with the Mississippian–Pennsylvanian unconformity, which is discussed in a later section.

There are observable trends in the vertical profile of data from the Adkisson #1–33, despite the limitations due to (1) age of the samples (late Paleozoic), (2) bulk sample type, and (3) diagenetic modifications. Overall, despite samples being derived from bulk sediment, some of the δ¹³C and δ¹⁸O values agree well with Mississippian seawater values derived from nonluminescent brachiopod shell data reported by Mii et al. (1999; Figure 5). Mii et al. (1999) report δ¹³C and δ¹⁸O stable isotope profiles that likely reflect unaltered Mississippian seawater due to their strict preservation criteria (thin section petrography, SEM imaging, cathodoluminescence, and trace element analysis) in brachiopod shell sampling.

There is a significant portion of the data, however, that plots outside of the suggested Mississippian seawater field, typically as more depleted δ¹³C and δ¹⁸O values. This effect is more pronounced in δ¹⁸O values, which likely reflects the enhanced susceptibility of oxygen isotope alteration during water–rock interactions (Lohmann, 1988; Banner and Hanson, 1990) and often leads to the interpretation of an unreliable δ¹⁸O curve. The data in this study, however, suggests that a pattern of depleted δ¹⁸O values (~−5‰) in shallow-water facies at the tops of shallowing-upward successions indicates an important chemostratigraphic relationship that may reflect mixed meteoric input near third- and potentially fourth-order sequence boundaries. Moving average curves were constructed to minimize noise related to oxygen isotope susceptibility and also support the hypothesis of higher frequency cycles, especially within larger, potentially third-order packages (Figure 7). Relatively depleted δ¹⁸O values are similar to a meteoric water line (~−5.5‰) established from soil-formed microcomponents in subaerially exposed horizons in the Mississippian Holder Formation in New Mexico (Goldstein et al., 1991; Figure 5). Although there are no significant subaerial exposure horizons or pedogenic features present in the Adkisson #1–33, recent work in Kingfisher County, directly northwest of this study area, has revealed multiple, well-developed exposure horizons characterized by dissolution pipes and chert breccia in a single core (Flinton, 2015). This suggests that Mississippian strata in this study area are likely to have been affected by some level of meteoric input and the difference (~0.5‰) in the δ¹⁸O data from the Adkisson #1–33 likely reflects a mixed value rather than a pure meteoric signal.

These findings are consistent with the diagenetic history established for the Adkisson #1–33, in which two types of early calcite cement were identified using thin section petrography and cathodoluminescence (Mohammadi et al., 2019a). Bladed, partially isopachous cement was reported most prevalent in Facies 4 and Facies 5 and interpreted to reflect rapid precipitation in a marine environment. Syntaxial blocky cement, reported most prevalent in Facies 5, was interpreted to represent precipitation in reducing, marine phreatic environment (Mohammadi et al., 2019a). The higher position of facies within the stacking pattern and enhanced probability of being affected by mixed marine–meteoric waters likely caused a relatively depleted final δ¹⁸O value to be recorded in some of these samples. This depletion is due to the fractionation between ¹⁶O and ¹⁸O that occurs during evaporation and precipitation processes wherein meteoric waters contain a higher percentage of ¹⁶O relative to ¹⁸O and results in a lighter, depleted δ¹⁸O value (Lohmann, 1988). The enhanced diagenetic susceptibility of skeletal-rich facies may also play a role in recording relatively depleted ¹⁸O values as syntaxial blocky cement most commonly occurs as overgrowths on crinoid and brachiopod fragments (Mohammadi et al., 2019a). Alternatively, micritic-rich facies near the base of third-order sequences are less commonly affected by isopachous and fringing cements (Mohammadi et al., 2019a). Basal facies tend to plot within the seawater field established by Mii et al. (1999) and suggest the recording of a relatively unaltered seawater signal.

In addition to the patterns observed at the third-order scale, tentative ages have been assigned to the midcontinent profiles of this study using a compilation of literature data sets hypothesized to reflect global carbon cycling (Grossman et al., 1993; Bruckschen et al., 1999; Mii et al., 1999; Saltzman, 2003; Batt et al., 2007; Buggisch et al., 2008; Koch et al., 2014). Although it was not possible to direct the current samples to selective cements or shells that have not undergone diagenesis, the patterns and trends of isotopic composition within the bulk samples are comparable to published curves of the individual Mississippian stages and show similar variations in isotopic composition that are often attributed to local paleogeographic factors in other Mississippian sections (Grossman et al., 1993; Bruckschen et al., 1999; Mii et al., 2001; Grossman et al., 2008). Figures 6 and 8 depict each of the proposed North American stages that are characterized by specific δ¹³C trends that resemble assumed global patterns.

The proposed correlation of the Adkisson #1–33 data set (depths 5819.85–5760.6 ft [1.77–1.75 km]) to the Kinderhookian stage records a positive δ¹³C excursion from −4.54‰ near the base to +2.53‰ near the top (difference +7.07‰). This is the most significant excursion of the data set and resembles the well-known positive shift in the early to middle Tournaisian (Bruckschen et al., 1999; Mii et al., 1999; Saltzman, 2003; Buggisch et al., 2008; Koch et al., 2014; Figure 8). This excursion, often covariant with δ¹⁸O, is hypothesized to represent a significant alteration in carbon cycling, the causality of which is strongly debated among different authors. Although the midcontinent excursion in this study begins with the most depleted δ¹³C of published values (−4.54‰) and only achieves a slightly positive maximum value (+2.53‰), the magnitude of the shift (~+7‰) is consistent with previous studies. Koch et al. (2014) report a similar excursion from −2.3‰ to +3.3‰.

A marked δ¹³C depletion from ~+1.5 to −0.25‰ (difference −1.75‰) occurs near the middle of the proposed stage and is consistent with a middle to late Kinderhookian negative excursion present in global profiles (Mii et al., 1999; Saltzman, 2003; Koch et al., 2014). Recent outcrop-based carbon and oxygen stable isotope work in the Tri-State region of northeastern Oklahoma, southwestern Missouri, and northwestern Arkansas suggests this negative excursion may be related to the deposition of the late Kinderhookian Northview Formation (Sessions et al., 2019). The excursion has been dated through the use of high-resolution biostratigraphy (Boardman et al., 2013) and similarities in the isotope curves suggests the presence of middle to late Kinderhookian strata in the Adkisson #1–33. The subsequent enrichment seen in the midcontinent data of this study is somewhat consistent with the slow enrichment of the late Kinderhookian Pierson Formation near +2.0‰ δ¹³C of Sessions et al. (2019). The similarity in observed trends, including the major positive excursion early in the stage, the negative excursion in the middle to late stage, and the stabilization near +2–3‰ later in the stage, suggests a correlation with established global profiles and recent Mississippian midcontinent work in the Tri-State region (Mii et al., 1999; Saltzman, 2003; Koch et al., 2014; Sessions et al., 2019).

The thickness of the proposed Kinderhookian unit (18.3 m [60 ft]) is consistent with the depositional model of diachronous, prograding carbonate wedges proposed for the northeastern Oklahoma, southwestern Missouri, and northwestern Arkansas Tri-State region for Lower to Middle Mississippian strata (Boardman et al., 2013). Outcrop-based sequence stratigraphy in this proximal position suggests the Kinderhookian thickness ranges from 3 meters (9.8 ft) to 5.5 meters (18 ft; Childress and Grammer, 2015, 2019); however, the relative distal positioning of the Adkisson #1–33 likely recorded a thicker section of Kinderhookian strata as prograding clinoforms continued to the southwest into the basin. In addition to variable sedimentation rates and preservation related to the more distal positioning, syndepositional faulting associated with the Nemaha uplift may also have created additional accommodation space that allowed sediment to be funneled into local depositional lows (LeBlanc, 2014) resulting in the deposition of a relatively thick Kinderhookian section. In their work to the northwest in Garfield County, Koch et al. (2014) also report a thick Kinderhookian interval (36.6 m [~120 ft]) based on global patterns recognized in δ¹³C.

The proposed boundary between the Kinderhookian and Osagean stages is characterized by a high gamma ray signature (>100 API units) and the occurrence of a very thin (6.1 cm [2.4 in]) bed of Facies 1 (glauconitic sandstone). Isotopic data does not show significant enrichments or depletions associated with this boundary; however, lack of geochemical evidence is common across this horizon, particularly in other midcontinent studies (Mii et al., 1999; Koch et al., 2014; Figure 8). Therefore, the assignment of this boundary relies primarily on gamma ray and facies data, both of which support a significant alteration in deposition from relatively well oxygenated Facies 3 (characterized by very low gamma ray values, <5 API units) to relatively restricted, oxygen-poor Facies 1 (characterized by high gamma ray values, >100 API units).

Using the criteria established by Amorosi (2012), the glaucony present at this horizon is characteristic of a simple omission surface. Preliminary energy-dispersive x-ray fluorescence data suggests the glaucony is relatively immature, characterized by low percentage (2–3%) of K₂O and is reflected in the samples light green “nascent” hue. Although the glaucony is not representative of a mega condensed, globally correlative surface, this may represent a relatively proximal expression of the same surface. In this case, the relatively thin (~6 cm [2 in]) deposit would not have experienced the same magnitude or duration of restriction and would have relatively immature glaucony compared to exclusively basinal deposits. Induced coupled plasma mass spectrometry and oriented XRD mounts are needed to confirm K₂O weight percentages and clay mineralogy with confidence. Additionally, Facies 1 is not present after its occurrence at this horizon, which may further support a significant deepening episode.

The proposed correlation of the Adkisson #1–33 data set (depths 5760.6–5703.7 ft [1.75–1.74 km]) to the Osagean stage is characterized by a relatively stable δ¹³C profile between +2.0 and +3.0‰ with high-amplitude changes (up to 1‰ shifts in 3.05 m [10 ft]) occurring in the uppermost 7.6 meters (25 ft; Figure 6). This relatively stable profile is similar to that of Koch et al. (2014) wherein late Tournaisian and early Visean strata maintain an average 1‰ range in δ¹³C values for the entire stage (Figure 8). Conversely, Buggisch et al. (2008) and Saltzman (2003) report stability only during the late Osagean and more common excursions (up to 2‰ δ¹³C) near the basal Kinderhookian boundary. The variability in the uppermost portion of both δ¹³C and δ¹⁸O profiles may be related to more frequent depositional cyclicity driven by sea-level change and enhanced facies variability due to the diagenetic susceptibility of shallower water facies. Although δ¹³C values are less susceptible to diagenetic alteration compared to δ¹⁸O (Banner and Kaufman, 1994), repeated exposure to relatively depleted, mixed marine–meteoric water may have caused a slight (0.5–1‰) lowering of δ¹³C. This is consistent with the systematic depletion of both δ¹³C and δ¹⁸O in the 1.5 meters (5 ft) below the Osagean–Meramecian boundary, which may be evidence of repeated exposure to depleted waters. Carbon isotopes also tend to be covariant with δ¹⁸O over this interval, which may suggest the influence of higher frequency, potentially fourth-order cyclicity (discussed in δ¹⁸O High-Frequency Cycles and Fourth-Order Potential). The δ¹⁸O curve also shows a relative depletion from base to top, which may reflect more mixed marine-meteoric input toward the top of the section overall. These isotopic trends, in conjunction with a cleaning upward gamma ray signature, are consistent with the overall shallowing-upward facies succession toward the end of the stage and the diagenetic history for shallow-water facies.

The proposed boundary between the Osagean and Meramecian stages is characterized by the above-mentioned depletion directly below and marked enrichment directly above in both δ¹³C and δ¹⁸O values (Figure 9). Although there is a general lack of geochemical evidence for this boundary in the literature, the enrichment may signify an important flooding surface. Katz et al. (2007) report δ¹³C curves of Lower Mississippian strata across third-order sequences that exhibit the highest values at maximum flooding surfaces in relatively distal positions within the depositional system. The proposed boundary is characterized by the most enriched values of the data set (+3.53‰ δ¹³C and −0.41‰ δ¹⁸O), records +1.6‰ and +4.79‰ shifts in δ¹³C and δ¹⁸O respectively, and based on facies present, occurs in a relatively distal position on the distally steepened ramp. Based on the newly defined chemostratigraphic framework, this enrichment may represent a horizon on the composite sequence scale, larger than the third-order scale excursions reported in Katz et al. (2007), but smaller than the single, second-order sequence of the entire Mississippian. Smaller positive enrichments (+0.5‰ to +1‰), similar in scale to Katz et al. (2007), also occur in the transgressive portions of the newly defined third-order sequences in the Osagean strata leading up to this boundary. This may represent several deepening upward episodes evidenced by high gamma ray signatures (50–100 API units) and the repeated occurrence of Facies 2. Overall, this additional evidence suggests the expression of relatively restricted, potentially deeper water horizons with the most significant occurring at the proposed stage boundary. This is consistent with regional findings of Curtis and Champlin (1959) wherein an early Meramecian marine transgression was suggested to have caused the submergence of most of Oklahoma and resulted in deposition of relatively restricted facies.

The proposed correlation of the Adkisson #1–33 data set (depths 5703.7–5585.4 ft [1.74–1.70 km) to the Meramecian stage is characterized by a relatively stable δ¹³C profile between +2.0 and +3.0‰ in the lowermost 9.14 m (30 ft), a steady depletion from +3.04 to +0.27‰ (difference −2.77‰) in the uppermost 23 m (75 ft), and a marked depletion from +2.05‰ to +0.44‰ (difference −2.49‰) occurring in the last 1.5 m (5 ft; Figure 6). A similar depletion is recorded in upper Meramecian and lower Chesterian strata with the most depleted δ¹³C values occurring near the Meramecian–Chesterian boundary (Mii et al., 1999; Saltzman, 2003; Batt et al., 2007; Koch et al., 2014; Figure 8). It has been suggested that the Meramecian–Chesterian δ¹³C decline may be due partly to the growth of the organic carbon reservoir and partly to intensified eastern Panthalassan upwelling in response to the closure of the seaway between Laurussia and Gondwana (Popp et al., 1986; Mii et al., 1999). Alternatively, a minor regression proposed by Curtis and Champlin (1959) in the late Meramecian may be the result of potential facies driven effects related to diagenetic susceptibility of shallow-water facies causing an overall depletion in δ¹³C values. This pattern is somewhat consistent with the current data set as the proposed Meramecian stage is characterized by the highest percentage of shallow-water facies and has the most consistently depleted δ¹⁸O values of any other stage. However, the late Meramecian of this study is characterized by a relatively higher percentage of deeper water facies and relatively enriched δ¹⁸O values compared to the early and middle portions of the stage. This suggests influences other than repeated exposure to mixed meteoric waters, at least at the top of the section near the Chesterian boundary. Global influence may have had more of an impact on isotopic signature at this portion of the stage, which is consistent with seawater signal preservation in deeper water facies.

The proposed boundary between the Meramecian and Chesterian stages is characterized by a marked depletion in δ¹³C (−1.74‰ over 1.2 m [~4 ft]) just below and a marked enrichment in δ¹³C (+1.55‰ over 0.91 m [~3 ft]) just above (Figure 6). A lack of a distinct shallowing-upward facies package and relatively enriched δ¹⁸O values beneath the boundary suggest a deviation from the pattern of potential mixed marine–meteoric input at third-order sequence boundaries seen downsection. Clean gamma ray signature (~15 API units), occurrence of Facies 4, and significantly depleted δ¹³C at the top of the Meramecian suggest a carbonate-dominated deposit that may have been altered due to repeated mixed marine–meteoric input; however, enriched δ¹⁸O values do not agree with this hypothesis. It may be possible that the relatively enriched δ¹⁸O values, which are consistent with Mississippian seawater field established by Mii et al. (1999), and slightly depleted δ¹³C values collectively represent microbial sulfate reduction despite their occurrence in a relatively shallow facies. This may be associated with the suggested deepening in the early Chesterian (Curtis and Champlin, 1959) evidenced by the large deposit of relatively restricted Facies 2 and high gamma ray signature (>50 API units). Slightly depleted δ¹³C values may also reflect times of enhanced oxidation of organic matter (Swart, 2015), which may be more consistent with glacioeustatic changes of the Late Mississippian. The discrepancy in the two curves may also suggest the dominance of a global δ¹³C signal during the late Meramecian, which may not be recorded in δ¹⁸O values.

The marked enrichment in δ¹³C above the boundary is of similar magnitude (+1.5‰) reported in late Visean–early Serpukhovian strata, which is hypothesized to represent a climatic link between enhanced organic carbon burial and potential onset of glaciation during the Serpukhovian (Buggisch et al., 2008). Additionally, basinal subsidence is suggested to have been much greater during the Chesterian relative to earlier Mississippian stages (Curtis and Champlin, 1959), which may have enhanced the distribution of relatively restricted facies and caused increases in the burial rate of organic carbon. This is consistent with the large (3.65 m [~12 ft]) deposit of Facies 2 present at the base of the proposed Chesterian section, which is also characterized by the highest gamma ray values (100–250 API units) of the entire data set. The marked depletion in δ¹⁸O values just above the boundary may be related to mixed marine–meteoric signatures in deep-water facies associated with rapid and high-amplitude sea-level fall. Large drops in sea-level related to the onset of polar glaciation, which were likely more frequent during Late Mississippian, may have caused deeper water facies to be uncharacteristically exposed to mixed marine–meteoric waters. Alternatively, the rapid depletion may be evidence of hydrothermal input (Goldstein et al., 2019); however, the magnitude of the shift would likely be in the −10‰ range as opposed to only −6‰ to −7‰. Recent fluid inclusion work (Mohammadi et al., 2019a, b) suggests the likelihood of at least some alteration by basinal fluids of midcontinent Mississippian sections, which may also cause relatively depleted isotope values. Additional analysis of this proposed boundary is needed to delineate the high-frequency changes and discrepancy between the two isotope curves.

The proposed correlation of the Adkisson #1–33 data set (depths 5585.4–5496.35 ft [1.70–1.68 km]) to the Chesterian stage is characterized by an early-stage, gradual enrichment in δ¹³C from +0.27‰ to +2.15‰ (difference + 1.88‰), a mid-stage stabilization at approximately +2.0‰ (~0.4‰ range), and a late-stage, gradual depletion from +1.89‰ to +0.97‰ (difference −0.92‰; Figure 6). This high-resolution pattern is somewhat unique compared to global data sets; however, the general trend of relatively depleted (~+1‰) δ¹³C values is similar to patterns recognized in the midcontinent (Mii et al., 1999; Koch et al., 2014; Figure 8). The base of the Chester in this study was defined by facies stacking patterns in conjunction with stable isotopic data and is interpreted here to be related to the moderate depletion (~2‰) seen in the Arrow Canyon data set (Figure 8). The more marked shift (~3–4‰) above this boundary observed in both Arrow Canyon and in the Anadarko Basin is not clearly discernible in this data set, although it is recognized that some shifting of the data is feasible based upon issues related to core depths being offset. The use of an integrated approach including facies stacking patterns, gamma ray values, and moderate depletion aided in the assignment of this proposed boundary. The gradual depletion in δ¹³C in the uppermost 10 meters (33 ft; Figure 6) may be related to an increase in higher order cyclicity; however, relatively enriched δ¹⁸O values again do not support the hypothesis of increased exposure to depleted, mixed marine–meteoric waters. Alternatively, the more frequent preservation of Mississippian seawater values suggests less diagenetic alteration in shallow-water facies as both δ¹³C and δ¹⁸O values tend to plot within the Mississippian seawater field of Mii et al. (1999). The lack of mixed marine–meteoric signature in Facies 5 may support a Chesterian section frequently affected by storm-related events. In this case, the repeated deposition of Facies 5 in a relatively distal position with a lack of mixed meteoric influence would produce a relatively stable δ¹⁸O profile consistent with Mississippian seawater values that have not been diagenetically altered. This is consistent with increased basinal subsidence suggested for Chesterian strata, thus this portion of the section may actually record a relatively deep deposit of large autocyclic packages between fair weather and storm wave base.

The most depleted δ¹⁸O value of the entire data set occurs just below the proposed Chesterian–Pennsylvanian boundary (−8.28‰) and records a negative shift of −5.25‰ over only 0.15 meters (~0.5 ft). δ¹³C values are also significantly depleted (−1.65‰) and record a shift of −2.62‰ over the same interval. Although there is evidence for subaerial exposure of the section in different parts of the region (Flinton, 2015; Jaeckel, 2016), no physical indication of exposure is observed in the Adkisson #1–33, although such exposure could also affect the depletion in the isotopic values. Extensive fluid inclusion work in Mississippian midcontinent rocks in south–central Kansas suggests late hydrothermal alteration of the upper Mississippian associated with the retardation of hot basinal brine migration at the low porosity, low permeability Pennsylvanian shale contact (Goldstein et al., 2019). Although this is characteristic of only one sample, the two samples above the boundary are also relatively depleted (δ¹³C: −1.17, −0.82‰ and δ¹⁸O: −4.43, −3.43‰), especially in δ¹³C relative to the majority of the data set.

Traditionally, δ¹⁸O data from bulk sediment without comparison to a suite of targeted sample types (cement, unaltered shells, micrite) has not been considered for characterizing stable isotope values of Mississippian strata. Although there is considerable variability in the δ¹⁸O data of this study, with values that range from −0.41‰ to −8.28‰ and can change rapidly (up to 5.25‰) over short transects (0.3 m [<1 ft]), there is a discernable pattern of high-frequency cycles that are characterized by relatively enriched bases (~−2.0‰) and relatively depleted tops (~−5.0 to −6.0‰). Forty-four of these cycles, characterized by depleted δ¹⁸O values near the top, have been identified and are hypothesized to reflect high-frequency, potentially fourth-order, shallowing-upward cycles that have been affected by mixed marine–meteoric input at cycle tops (Figure 6). Shallowing-upward cycles tend to be capped by shallow-water facies and plot in relatively depleted fields relative to the Mississippian seawater box (Figure 10). The suggestion of repeated exposure to relatively depleted waters is consistent with subregional findings of Flinton (2015) wherein multiple subaerial exposure horizons were identified in a single core. Although there are no exposure horizons present in the Adkisson #1–33, the potential for exposure to mixed marine–meteoric waters is probable within this area of the depositional system and likely occurs at multiple scales related to the stratigraphic hierarchy.

High-frequency cycles may or may not be characterized by a cleaning upward gamma ray signature, which suggests a resolution limitation of gamma ray curves at this scale. Although some cycles are picked more confidently than others, patterns seen in δ¹⁸O values may be extremely useful in defining higher order cyclicity that may not be predicted using facies and gamma ray data alone, especially in distal areas of the depositional system where subaerial exposure surfaces are lacking. This has important consequences when considering the potential for cemented, flow baffling horizons that can compartmentalize reservoir units and may not be identified based on stacking patterns alone.

Overall, the relative distal position of the study area has important consequences for isotopic signatures reported herein. It is difficult for samples only experiencing marine diagenesis (and potentially some mixed meteoric) to become significantly altered from original δ¹³C values, thus facies deposited below average sea-level fall for the particular period are preferred for chronostratigraphy because they avoid most diagenetic effects (Swart, 2015). This suggests the δ¹³C curve in this study may indeed retain a global signal despite being derived from mixed input, bulk sediment. Although there are shallow-water facies present, there are no significant exposure horizons in the Adkisson #1–33, thus the study area may represent a relatively ideal location for developing a chronostratigraphic tie. Assuming a full Mississippian section, the Adkisson #1–33 also lacks significant exposure despite the switch from probable low-amplitude (10 m [33 ft]) sea-level change during Early Mississippian to high-amplitude (75–100 m [246–328 ft]) changes during Late Mississippian (Read, 1995), unlike cores in neighboring Kingfisher County where several horizons are developed in a single core. This may support syndepositional faulting associated with the Nemaha uplift wherein the downthrown block (Logan County) created more accommodation space and funneled sediment into a local depositional low whereas the upthrown block (Kingfisher County) experienced simultaneous episodes of subaerial exposure (LeBlanc, 2014). Evidence for this type of subregional variability in diagenetic history, indicated here by carbon and oxygen isotopic values, may be crucial for understanding reservoir compartmentalization and porosity evolution related to meteoric and mixed marine-meteoric input over relatively short distances.

The integration of high-resolution carbon and oxygen stable isotope profiles with gamma ray curves, high-resolution sequence stratigraphy, and diagenetic history has revealed key patterns that can be used to delineate potential global signatures as well as probable Milankovitch band cyclicity with a chemostratigraphic approach. A refined interpretation of the Adkisson #1–33 suggests more frequent third-order cyclicity than previously defined which is more consistent with global ties to the individual North American stages and within the Mississippian overall. Incorporation of biostratigraphic data (Hunt, in progress) will help further constrain the accuracy of these isotopic trends and further aid in correlation efforts.

Although there are a significant number of data points that plot outside of the Mississippian seawater field established by Mii et al. (1999), the incorporation of the diagenetic history suggests new fourth-order boundaries related to mixed marine–meteoric input at cycle tops. Overall, this high-frequency pattern recognized in δ¹⁸O values may provide important insight where facies data cannot, especially in distal areas of the depositional system where subaerial exposures are lacking. A fundamental component in understanding reservoir distribution is understanding the level at which subsurface units can be compartmentalized, especially in a mixed lithology system. The integrated chemostratigraphic approach applied in this study provides key insight in enhancing the predictability of subsurface units in terms of reservoir heterogeneity and how it relates to geochemical trends and high-resolution sequence stratigraphy.

This research was supported by the Oklahoma State University Industry Consortium on the Reservoir Distribution and Characterization of the midcontinent Mississippian Carbonates—A Major Unconventional Resource Play, the Oklahoma Geological Foundation, the Oklahoma City Geological Society, and AAPG Grants-in-Aid. Devon Energy provided access to core and associated wireline log, XRD, and TOC data sets. Carbon and oxygen stable isotopes were measured at the University of Miami Stable Isotope Laboratory—Rosenstiel School of Marine and Atmospheric Sciences (RSMAS). OSU Mississippian consortium members include: American Energy Partners, Chaparral Energy, Chesapeake Energy, Devon Energy, Longfellow Energy, Marathon Oil, Maverick Brothers, Newfield Exploration, SM Energy, Samson Energy, Sinopec (Tiptop), Redfork Energy, Trey Resources, and Unit Petroleum. This is Boone Pickens School of Geology Contribution 2017-54.