Interpretation of Paleozoic paleokarst features in the Arkoma Basin of Oklahoma using 3D seismic and well logs Open Access

Paleokarst regions worldwide are repositories for hydrocarbons, mineral deposits, and groundwater. Time structure maps were generated for the Ordovician Viola Limestone, Mississippian Caney Shale, and Pennsylvanian Jefferson Sandstone and Wapanucka Limestone. Isochron maps indicate pronounced visible sinkhole time thickening in the Viola-Caney and Caney-Jefferson intervals relative to the Jefferson-Wapanucka. Sinkhole features in the Viola exhibit mappable structural depression, characterized by lower positive amplitude, higher seismic variance, and most-negative curvature. Curiously, spatially coincident sinkhole features in the shallower Wapanucka display the opposite characteristics relative to adjacent areas that have not been modified, namely, higher positive amplitude and lower seismic variance with no mappable time structure relief. Seismic amplitude analysis based on well logs and Gassmann modeling indicate that the Viola has a reduction in limestone acoustic impedance inside sinkholes that allows estimation of increased porosity near 10%. Identical analysis for the Wapanucka suggests that no reasonable alteration of the limestone acoustic impedance alone can account for the observed amplitude behavior, implying that the limestone and overlying shale must be altered in sinkhole areas. Some of these interpreted sinkhole features coincide with vertical pipe structures with up to 490 m (1610 ft) vertical extent, diameter up to 520 m (1700 ft), and separation of at least 460 m (1510 ft). We interpret the Viola sinkhole features and associated vertical pipes to be part of a mature epigene karst system. Conversely, the shallower and more subtle Wapanucka sinkholes we interpret as related to an immature mixed karst system with epigene and hypogene elements. Our study indicates for the first time the seismic evidence of pipe features that extend both below and above the Viola and the presence of Wapanucka sinkhole features in the Arkoma Basin of Oklahoma, which provides a better understanding of paleokarst occurrence and its possible impact on resource exploration.

Paleokarst is karst that is not hydrologically connected to the current earth’s surface and buried by younger sediments (Ford and Williams, 2007). Hydrocarbons, minerals, and groundwater are found in paleokarst reservoirs. Paleokarst furnishes information about past geologic and hydrologic conditions, sea level and climatic changes (Palmer and Palmer, 2011), can cause damage to property and civil engineering works (Waltham and Fookes, 2003), as well as lost circulation and complete loss of mud in hydrocarbon drilling (Andre and Doulcet, 1991; Lomando et al., 1993; Zhao et al., 2014).

James and Choquette (1988) note that the development of karst landforms occurs by external and inherent factors. The external factors include climate (precipitation, evaporation, and temperature), base level (relief and elevation, sea level, or local water bodies), plant life, and duration of time; and inherent factors such as structure and stratigraphy (strata attitude, unconfined or confined aquifers, and structural conduits) and lithology (fabric and texture, bedding thickness, fractures, stratal permeability, mineralogy, and bulk purity).

Two broad categories of karst development are recognized, called epigene and hypogene (Palmer, 1991; Klimchouk, 2015). Epigene, or meteoric, karst is associated with an unconformity surface involving the interaction of meteoric water with carbon dioxide from soil organic matter to form carbonic acid originating at or close to the earth’s surface (Palmer, 2007; Klimchouk, 2015; Milad and Slatt, 2017). Soil biogenic activity increases with temperature in terrains at low altitudes and latitudes, such as humid, temperate, and tropical regions. Hypogene karst is associated with carbonate dissolution by confined, deep-seated hydrothermal fluids (Palmer, 1991; Loucks, 1999; Klimchouk, 2007, 2009a, 2009b), or the acceleration of epigene processes can also produce similar fluids (Palmer, 1991). Those fluids can include sulfuric acid breached from hydrocarbon oilfields (Hill, 1995), as well as high-temperature and -pressure igneous basement hydrothermal fluids migrated along faults (Palmer, 1991; Burberry et al., 2016), or thermal convection of hydrothermal fluids (Wright and Harris, 2013). Unlike epigene karst, hypogene processes are independent of climate (Klimchouk, 2009a, 2009b).

Sinkholes are closed depressions of subsurface drainage diagnostic of epigene karst topography (Waltham et al., 2005). Collapse breccias with infill sediments are often present in sinkholes (Loucks, 1999), as well as open shafts into cave networks (Waltham et al., 2005). Figure ¹ illustrates a generalized karst model. Vertical karst pipe structures become connected by hydrothermal flow, tectonic activity, mineralization, and collapse (Waltham et al., 2005; Sullivan et al., 2006; Sun et al., 2013; Burberry et al., 2016) with hypogene-formed sinkholes enhanced during subaerial exposure (Sullivan et al., 2006; Burberry et al., 2016). Wright (2016) notes that the circular collapse features that occur in hypogene networks may be interpreted as surface sinkholes originally associated with meteoric karst. Sinkholes and associated pipe features have been identified from 3D seismic data in the Fort Worth Basin (Hardage et al., 1996; McDonnell et al., 2007), the Persian Gulf (Burberry et al., 2016), the Pearl River Mouth Basin, China (Story et al., 2000; Heubeck et al., 2004; Sun et al., 2013), and Florida (Cunningham and Walker, 2009; Cunningham, 2015; Cunningham et al., 2018). In 3D seismic data, karst pipes are often seen to narrow upward, develop a cylindrical to vertical conical geometry, and exhibit a spectrum of disruption of stratal seismic reflections from localized sag features to completely chaotic (Cartwright et al., 2007; Sun et al., 2013). Pipes are believed to have formed over an extensive time period (Waltham et al., 2005).

There have been limited seismic studies of paleokarst in the Arkoma Basin, Arkansas-Oklahoma. Brinkerhoff (2007) uses a waveform classifier to distinguish the various stages of karsting, specifically paleocave development, incipient karst collapse, and noncollapse regions in the Ordovician-Devonian Hunton Limestone in the Arkoma Basin of Oklahoma. Moser (2016) uses curvature to map sinkholes with an average diameter of 630 m (2070 ft) in the Mississippian Boone Limestone in the Arkoma Basin of Arkansas. Milad and Slatt (2017) map sinkholes in the Hunton and Viola Formations in Pottawatomie County, Oklahoma, on the Cherokee platform 92 km (57 mi) northwest of the current study area, finding sinkhole diameters of 350–700 m (1150–2300 ft). Using the same 3D seismic survey as the current study, Kumbalek (2015) maps and identifies Viola paleokarst expressed as sinkholes with an average diameter of approximately 280 m (780 ft) that occurred in only 4.1% of the $460km2$ (approximately $180mi2$⁠) survey area.

This paper reports the first seismic mapping in the study area of the Mississippian Caney Shale, the Pennsylvanian Jefferson Sandstone, and the Pennsylvanian Wapanucka Limestone, as well as the Ordovician Viola Limestone. We identify and measure sinkhole and vertical pipe features in the Arkoma Basin of Oklahoma using horizon time structure maps and isochron maps, optimized seismic attribute volumes of variance, curvature, and amplitude maps, sinkhole feature amplitude analysis is calibrated to the Gassmann equation to form a predictive rock-physics model, and we extend Kumbalek’s (2015) Viola sinkhole analysis. Evidence is presented for paleokarst collapse that extends below the Viola and into shallower horizons and the first published description of sinkhole features in the Wapanucka Limestone. This study has broad applications in paleokarst science and hydrocarbon exploration.

The study area is located in the western Arkoma Basin, a peripheral foreland basin formed by collision of the North American and Gondwanan plates during early Mississippian through middle Pennsylvanian time (Suneson, 2012). It is a structural-sedimentary basin that covers much of eastern Oklahoma and western Arkansas and stretches south to the Choctaw Fault (Amsden, 1984). Figure ² shows the study area, Arkoma Basin and adjacent basins along the Ouachita fold belt. Figure ² shows a cross section across the Arkoma Basin and the Ouachita fold belt. Surface rocks of the western Arkoma Basin trend east–northeast with regional northwest dip (Berry and Trumbly, 1968). The youngest beds are visible on the northwest edge of the basin, whereas the oldest beds occur near the Choctaw Fault. The Wapanucka Limestone and older rocks dip regionally to the southeast (Berry and Trumbly, 1968). Depositional thinning in the Atoka and younger formations are evident in growth structures (Berry and Trumbly, 1968), whereas south-dipping faults cut through early Pennsylvanian and older rocks to define the basin (Perry 1994). Compressional folds show substantial structural changes in the southern region of the basin adjacent to the Ouachita front (Berry and Trumbly, 1968; Suneson, 2012), and drape anticlines are present in the northern Arkoma Basin over normal faults (Suneson, 2012).

The deposition of Cambro-Ordovician Arbuckle Dolomite and basal sandstone occurred in a gradually subsiding platform near a geosyncline located to the southeast receiving some input of coarse clastics. During Simpson time, the region was subjected to an influx of clastic sediments that formed the mid-Ordovician Joins and Oil Creek sandstones and shales, and their northern equivalents, the Burgen-Tyner sequence. In the south and southeastern shelf area, carbonate production was high during McLish and Bromide times with marginal amounts of shale and sandstone. The absence of coarse clastic rocks and a stable platform initiated the deposition of the upper-Ordovician Viola Limestone and Sylvan Shale, and the Silurian-Devonian Hunton Limestones (Arkoma Basin Study Group, 1961). In the study area, the Viola has an average thickness of 50 m (163 ft).

The deposition of the Viola Limestone occurred in an extensive shallow epicontinental sea with no apparent orogenic activity in south-central Oklahoma (Wengerd, 1948; Mairs, 1966). After the deposition of the Viola, the epeirogenic tectonic movement exposed the Viola Limestone to subaerial erosion (Wengerd, 1948). Sykes et al. (1997a) and Sykes (1997b) note that the timing of karst development in the Viola is pre-Pennsylvanian in age, with vugs, solution-enlarged fractures, and channels in the upper Viola (Welling/Fernvale) suggesting dissolution before deposition of the Sylvan. The presence of sphalerite, copper sulfides, and pyrite with asphalt has been reported in the Viola west of the study area in Pontotoc County (Sykes et al., 1997a), indicating some hydrothermal activity and associated hypogene karst likely due to movement of brine and petroleum below the organic Sylvan Shale, a confining and effective aquitard unit over the Viola.

A marine transgression led to deposition of the Sylvan Shale (Amsden, 1984), with an average thickness in the study area of approximately 29 m (96 ft). The overlying Hunton Limestone does not exceed 2 m (8 ft) thickness in the study area and thins from the southwest to the northeast due to local erosion. Shelf subsidence of the Hunton surface resulted in the buildup of the transgressive Misener Sandstone and Woodford Shale, with the Woodford thickness of 52 m (171 ft). The Mississippian Mayes-Caney Shale sequence records a clastic advancement from the south, with a rapid thinning of the shale to the north across the platform due to slower deposition and rapid subsidence of the basin to a southward thickening of the Pennsylvanian Caney Shale, also called the Goddard or Springer Shale (Elias, 1956; Arkoma Basin Study Group, 1961). The average thickness of the Mayes-Caney Shale is approximately 146 m (482 ft) and the Goddard Shale is approximately 54 m (176 ft).

Jefferson Sandstone lenses are found at the edge of the platform in the upper section of the Pennsylvanian Springer/Caney Shale (Arkoma Basin Study Group, 1961). The Jefferson is composed of more than one sandstone that divides and amalgamates suggesting bar facies and variable depositional surroundings in a marine environment (Andrews, 2007a, 2007b). The Jefferson Sandstone has an average thickness of 34 m (110 ft) in the study area. Cromwell Sandstone deposition occurred with amplified movement to the south in a stable environment. Thin shale streaks in the sandstone indicate variability in subsidence rates. The average thickness of the Cromwell Sandstone is approximately 50 m (164 ft).

The Wapanucka Limestone formed in shallow waters before initiation of basin subsidence characterized by superficial and localized movements of the seafloor with a slow rate of deposition (Arkoma Basin Study Group, 1961; Suneson, 2012). The average thickness of the Wapanucka is approximately 47 m (153 ft). Before the advancement of the Atoka Sea, Morrowan rocks were subject to erosion that increased northward across the basin. During Atokan time, the deposition of coarse clastic rocks occurred throughout the basin with increased subsidence. A northward transgression occurred depositing younger shallow marine sands and shale over older Atoka units in the subsiding trough (Arkoma Basin Study Group, 1961). Figure ³ shows the stratigraphy of well C in the study area.

The 3D seismic and wireline data from three wells used in this study were made available by Devon Energy. Figure ⁴ shows the seismic data coverage and key well locations. Table 1 shows the formation tops and thicknesses encountered in the three wells. The data straddle the Hughes-Coal County line in southern Oklahoma. The seismic data have a 1 ms sample rate, 2.7 s record length, and bin size of $33.5×33.5m$ (⁠$110×110ft$⁠), and they consist of prestack time-migrated data with 798 east–west crosslines and 698 north–south inlines. The processing datum is 274.3 m (900 ft) with a replacement velocity of $3048m/s$ (⁠$10,000ft/s$⁠) and areal coverage of $470km2$ (⁠$180mi2$⁠). Fourier analysis indicates the minimum and maximum frequencies of 13 and 102.5 Hz at negative 20 dB, with a dominant frequency of 57.5 Hz (Figure ⁴). Vertical seismic resolutions for the Viola and Wapanucka are 27 m (89 ft) and 25 m (82 ft), respectively. Well A is located in Hughes County with total depth (TD) of 1847 m (6059 ft) in the Woodford Shale. Well B in Coal County had TD of 2417 m (7931 ft) to the base of the Viola Limestone. Well C, also in Coal County, had TD in the Simpson Group (McLish) at 2469 m (8102 ft).

Three wells, herein called A, B, and C, were used to correlate seismic events to geologic formation tops. Table 1 shows the formation tops and thicknesses in the three wells. A synthetic seismogram generated in well B is shown in Figure ⁵. This well was used because it had a long run of sonic and density that reached the Viola. Check-shot data were not available. A zero-phase 200 ms wavelet (taper 25 ms) was extracted in an 1100 ms time window based on field traces in a $10×10$ bin area centered on the well location. A time shift was applied to the synthetic to match the field seismic data, but no stretch/squeeze was required.

Wireline log plots were generated for the Viola and the Wapanucka Formations over a 91 m (300 ft) interval. This interval started 100 ft (30m) above the carbonate formation tops in wells C and B. Figure ⁶ shows the gamma ray (GR) and mineralogy rock fractions in well C. Figure ⁷ shows GR, acoustic (DTCO), and shear (DTSM) velocities in well B.

The mapped horizons of interest are shown in the yellow circles labeled V, C, J, and W representing the Ordovician Viola Limestone, Mississippian Caney Shale, Pennsylvanian Jefferson Sandstone, and Pennsylvanian Wapanucka Limestone, respectively (Figures ⁵, ⁸, ⁸, ⁹, and ¹⁰). The horizons are positive amplitude reflections. Time structure maps were generated for the interpreted intervals (Figure ¹¹–¹¹), along with isochron (time-thickness) maps for intervals among the Viola-Caney (VC) (Figure ¹²), Caney-Jefferson (CJ) (Figure ¹²), and Jefferson-Wapanucka (JW) (Figure ¹²).

To optimize imaging of karst-related features, seismic attribute parameter tests for variance were performed on a 1400 ms cropped seismic amplitude data volume (⁠$54km2$ [$21km2$]) covering wells B and C. All variance calculations used a $3×3$ bin operator. Two triangular weighted time filters were tested (five and 15 samples), as well as with/without dip correction of two types (horizontal variance and variance computed along a dipping plane). The dipping plane method uses principal component analysis (PCA) with a directional parameter (inline, crossline, and vertical scale) of 1.5 and a 0.6 plane confidence threshold. PCA dip correction was for confidence >0.6, whereas other regions were processed with horizontal variance dip correction (Gersztenkorn and Marfurt, 1999). Dip guided smoothing was the final variance parameter tested in conjunction with the operator size length and dip corrections. In total, six variance volumes were computed and examined for optimum detail at the Viola (Figure ¹³) and Wapanucka (Figure ¹⁴) horizons. A visual inspection determined that the optimum variance parameters for the Viola were those of Figure ¹³, which were then applied to the entire survey to extract horizon slices for the Viola, Caney, and Jefferson (Figure ¹⁵–¹⁵). The variance parameters of Figure ¹⁴ were deemed optimum for the Wapanucka and were applied to the entire survey to generate the horizon slice of Figure ¹⁵.

The most-positive and most-negative curvature volumes (Chopra and Marfurt, 2007) were computed using an operator size of (nt, nx, ny) = (12, 1, 1). Horizon slices were extracted from each curvature volume on four interpreted horizons (Figure ¹⁶–¹⁶). Amplitude maps were generated for the four horizons (Figure ¹⁷–¹⁷). Rose diagrams of faults/lineaments for the Viola and Wapanucka from the curvature volumes are shown in Figure ¹⁸.

Wireline logs in well C show that the carbonate rock fraction is higher in the Viola (Figure ⁶) than the Wapanucka (Figure ⁶): Specifically, mineralogy fractions indicate that the Viola has an average of 88% carbonate, 9% quartz, and 3% clay, whereas the Wapanucka average composition is 77% carbonate, 14% quartz, and 8% clay. Other mineralogy fractions are negligible. Bogli (1980) notes that the presence of impurities such as clay and quartz in limestone lowered the capacity for karstification, implying that the Viola Limestone has greater karst potential than the Wapanucka. However, we note that other factors such as climate, hydrology, and the structural setting may be preponderant. The mineralogy logs also indicate the presence of coal in the shale section overlying the Wapanucka.

The Viola has an average acoustic velocity (DTCO) of $6224m/s$ (⁠$20,420ft/s$ or $48.97μs/ft$⁠) and an average shear velocity (DTSM) of $3216m/s$ (⁠$10,551ft/s$ or $94.78μs/ft$⁠) (Figure ⁷). The Wapanucka average DTCO is $5872m/s$ (⁠$19,265ft/s$ or $51.91μs/ft$⁠) and an average DTSM is $3072m/s$ (⁠$10,079ft/s$ or $99.22μs/ft$⁠) (Figure ⁷). The higher velocities for acoustic and shear in the Viola compared to the Wapanucka are consistent with a higher carbonate fraction. The results from the mineralogy rock fractions and velocities indicate that the Viola has a higher potential for karst development than the Wapanucka.

In this paper, the term “pipe” refers to a disrupted, semichaotic volume of seismic data, “sinkhole” means a concave-upward depression across a seismic reflection event occurring in carbonate, and “sag” means a quasicircular concave-upward depression in siliciclastic rocks.

Figure ⁸ and ⁸ shows dip and strike geoseismic sections, respectively, through well B, which were used to generate the synthetic seismogram of Figure ⁵. The mapped horizons are shown in the yellow circles labeled as V, C, J, and W representing the Viola Limestone, Caney Shale, Jefferson Sandstone, and Wapanucka Limestone, respectively. All of these horizons are positive polarity events representing a soft to hard response at the formation boundary. Vertical pipe features are indicated by bounding dashed white lines and are visible on the horst block (Figures ⁸, ⁹, and ¹⁰) but not adjacent graben blocks (Figure ⁸). We observe that sinkholes and sags are often associated with these pipe features.

The faulting architecture consists of normal faults with drags and folds. The faults compartmentalize the section into horst and half-grabens (Figures ⁸–¹⁰). In map view, the predominant faults strike northeast–southwest, and other strike orientations include west–east, northwest–southeast, and north–south. We observe that these faults compartmentalize the study area into five separate fault blocks (Figures ¹¹, ¹², ¹⁵, ¹⁶, and ¹⁷). For convenience, the fault blocks are named beginning from the north to the south of the study area as follows: HG1, G1, H1, HG2, and HG3, where HG is a half-graben, G is a graben, and H is a horst.

From rose diagrams in Figure ¹⁸, two principal fault orientations are seen in the Viola along N40°E–N50°E and N50°E–N60°E, whereas only one is evident for the Wapanucka along N50°E–N60°E. This follows the regional trend of the Ouachita fault. This predominant northeast axis along the regional trend of the Ouachita fault suggests that the orientation of tectonic activity is consistent from Ordovician Viola time to Pennsylvanian (Morrowan) Wapanucka time.

Sinkholes on the Viola and sags on Caney reflections are sometimes associated with pipe features. The pipes are subvertical with a probable narrowing upward (Figures ⁸–¹⁰). Sinkhole and sag features are visible above these pipes. Low amplitudes and disrupted reflections characterize the internal configuration of the pipe. Similar features are known in the Pearl River Mouth Basin (Sun et al., 2013) and Fort Worth Basin (McDonnell et al., 2007). In our data, some pipes extend above 1.2 s, below the Viola into the Simpson Group (Figures ⁸, ⁹, and ¹⁰), and possibly extend downward to the acoustic basement. However, this is not clearly visible on the amplitude section due to deep image and resolution limits. Outside of pipe features, reflections show greater continuity. As expected, the volumetric variance is greater in the pipes than adjacent undisturbed data volumes (Figures ⁹ and ¹⁰). Taken together, these observations suggest rock fracturing and/or dissolution. These low-amplitude, high-variance pipe features indicate collapse and infill, which we interpret as probable breccia pipes (Waltham et al., 2005).

Away from the pipes, the surrounding host rock has low variance and consistent amplitude, which we interpret to be unkarsted rock that has not undergone significant dissolution or collapse. Vertical faults are likely bounding the sag/pipe features and may have served as conduits for deep-seated hydrothermal fluids migrating during the Ouachita Orogeny in the Pennsylvanian (Kupecz and Land, 1991), or they may be due to meteoric water that percolated along fracture networks enhancing carbonate dissolution. Fracture networks may be linked with vertical faults bounding pipe and sinkhole features.

The structural high in the study area rises toward the west from the H1 horst block. The Viola and Caney maps (Figure ¹¹ and ¹¹) show structural relief that highlights circular to elliptical sinkhole features. The Jefferson Sandstone and Wapanucka Limestone time maps (Figure ¹¹ and ¹¹) do not exhibit any mappable sag or sinkhole structural relief.

Isochron maps for the VC, CJ, and JW are shown in Figure ¹²–¹² with sinkholes features indicated by the red arrows. VC and CJ isochores show thinning over sinkhole features. Subtle visible lineaments trending N70°E are seen in the CJ isochron in juxtaposition with sinkholes and sags. The JW isochron shows very subtle sag/sinkhole features indicating that collapse and dissolution may have been active during this interval.

On the Viola horizon, an optimum variance was achieved with a 15 ms vertical window and dip-guided smoothing, bringing out fine detail on sinkholes in the red oval area of Figure ¹³ relative to the other parameter choices. For the Wapanucka, optimum variance parameters were 5 ms vertical window and dip guided smoothing (Figure ¹⁴). For both horizons, variance shows faults with a higher definition than the associated horizon time or amplitude map. The Viola (Figure ¹⁵) and Caney (Figure ¹⁵) show high variance inside, and low variance outside, the sinkhole and sag features. A possible subtle circular feature is observed on the H1 block on the Jefferson Sandstone (Figure ¹⁵). Sinkholes in the Wapanucka (Figure ¹⁵) do not show a well-defined variance compared to the Viola sinkholes or Caney sags. The Wapanucka variance is high around the edges of the sinkhole and low within sinkholes.

Curvature accentuates faults in the study area, showing up-thrown fault blocks with positive curvature, and downthrown fault blocks with negative curvature. Viola and Caney horizon corendered most-positive and most-negative curvature maps (Figure ¹⁶ and ¹⁶) reveal positive curvature on the rim of sinkhole/sag features and negative curvature inside them. For the Jefferson Sandstone, the curvature shows some subtle evidence of sags (Figure ¹⁶) The curvature maps of the Caney, Jefferson, and Wapanucka reveal northeast–southwest lineaments expressed on the H1 horst that are also visible on the CJ isochore (Figure ¹², the yellow arrows).

Figure ¹⁷–¹⁷ shows horizon amplitude for the Viola Limestone, Caney Shale, Jefferson Sandstone, and Wapanucka Limestone, respectively, which reveals clear sinkhole/sag features on all horizons except the Jefferson. We have noted elsewhere (Aboaba and Liner, 2018, 2019) that the Viola amplitude (Figure ¹⁷) shows strong positive outside sinkholes and very low to negative within sinkholes. Conversely, the Wapanucka amplitude is seen to be weak positive away from sinkholes and stronger positive inside sinkholes. In the vicinity of well B, we were able to combine amplitude, log data, and the Gassmann (1951) theory to investigate these relationships as explained below.

Table 2 shows the analysis results for acoustic impedance, reflection coefficient, and amplitude. The key results are: The Viola calculates to a 27% impedance decrease from sinkhole exterior to interior (out-to-in), whereas the Wapanucka calculates to an 82% impedance increase from out to in.

Further investigation using the calibrated Gassmann curve for the Viola (Figure ¹⁹) shows that an increase in porosity of approximately 10% can account for the computed impedance drop inside sinkholes. Thus, the inferred acoustic impedance drop for the Viola is consistent with a reasonable porosity increase related to karst activity leading to sinkholes.

The calibrated Wapanucka Gassmann plot (Figure ¹⁹) shows that the maximum limestone acoustic impedance does not exceed 18 SI, but our estimated impedance from well B and amplitude ratio is 28 SI for the sinkhole interior. We conclude that no reasonable alteration of the Wapanucka Limestone by itself can explain the observed amplitude behavior. It follows that amplitude brightening seen in Wapanucka sinkholes requires softening (reduced AI) of the overlying shale, perhaps indicating hypogene karst hydrothermal activity not active in the Viola interval. We acknowledge that amplitude pattern behavior is only indirect evidence of hydrothermal activity.

In map view, sinkholes are circular to elliptical features that occur in all the fault blocks for the Viola Limestone (Figures ¹¹, ¹⁵, ¹⁶, and ¹⁷). The sinkholes in HG2 and HG3 are adjacent to the north-bounding faults of these blocks, and not as numerous to G1 and H1. There is an alignment of sinkholes with the major faults and lineaments. Visual inspection of the mapped intervals on the time (Figure ¹¹), variance (Figure ¹⁵), curvature (Figure ¹⁶), and amplitude (Figure ¹⁷) maps reveal that the Viola Limestone has the greatest sinkhole development of the studied horizons. Sags are poorly developed in the Caney Shale (Figures ¹¹, ¹⁵, ¹⁶, and ¹⁷). The Jefferson Sandstone shows no visible sags on the time or amplitude maps (Figures ¹¹ and ¹⁷), although subtle sags may be indicated on the variance and curvature maps (Figures ¹⁵ and ¹⁶). On the Wapanucka horizon, no sinkholes are evident on the time structure or curvature (Figures ¹¹ and ¹⁶) but are visible on the variance and amplitude (Figures ¹⁵ and ¹⁷). The Wapanucka sinkholes appear to be in the same location as pipe features that show no visible continuation into the Wapanucka on the vertical seismic sections (Figures ⁹ and ¹⁰). Note that these pipes are not seen everywhere in the survey area and are more prominent on the H1 block.

Modern sinkholes with diameters greater than 100 m have been documented in the Yucatan Peninsula, Mexico; Devil’s Sinkhole, Texas and Southern China (Palmer, 2007); and Papua New Guinea, Madagascar and Puerto Rico (Waltham, 2005). We observe that sag diameters are smaller in the Caney compared to sinkholes in the Viola. We relate this to the narrowing of the pipe features toward the top of the pipe. The sag diameters in the Caney range from 93 to 305 m (304–1000 ft) and 195 to 606 m (640–1990 ft) in the Viola. The depression reliefs measured within the Caney sags are approximately 11–33 m (36–110 ft) and 28–49 m (93–160 ft) in the Viola. As previously stated, we observed no sags or sinkholes with measurable time relief in either the Jefferson or Wapanucka.

The pipes originate within the carbonate section (Viola Limestone and below) implying regions of paleokarst, with no evidence of bright spots associated with collapsed paleocave sediments, for example, in the Tarim Basin, China (Zeng et al., 2011a, 2011b). The scale of the pipe features is 150–520 m (500–1700 ft) in diameter, spaced 460–2130 m (1500–7000 ft) apart, and a vertical extent of 213–490 m (700–1600 ft). We acknowledge that the pipe vertical extent reported here may be considered a minimum due to low seismic data quality below the Viola. Similar pipe features described in the Fort Worth Basin (Hardage et al., 1996; Sullivan et al., 2006; McDonnell et al., 2007) have a vertical extent of 760–1100 m (2500–3610 ft), in the Persian Gulf (Burberry et al., 2016) an extent of 1490–2100 m (4900–6900 ft), and in the Pearl River Mouth Basin of China (Sun et al., 2013) an extent of 100–1000 m (330–3300 ft) is reported.

Waltham and Fookes (2003) propose an engineering classification for karst recognizing juvenile, youthful, mature, complex, and extreme categories. We use the term “mature” and “immature” to classify paleokarst features: Mature karst exhibits large sinkholes and collapse features commonly found in both temperate regions, and the wet tropics, whereas immature designates juvenile and youthful karst. Juvenile karst is formed in impure carbonates, or at deserts and periglacial zones with rare sinkholes, and youthful karst formed in temperate regions has small sinkhole features.

The higher distribution and greater development of sinkholes in the Viola with time structure relief suggest that the Viola Limestone is a more mature karst system than the immature karst of the Wapanucka. We interpret that these Viola sinkholes to be dissolution/collapse sinkholes, or cockpit karst as found in a tropical environment (Kumbalek, 2015) that formed by the lowering of the Viola Limestone surface (Waltham et al., 2005). Factors promoting karst maturity during Viola time may include clean, pure, high-strength limestone, possible long exposure, and biogenic soil gas interacting with meteoric water to form a more aggressive fluid.

We do not expect to observe sinkhole formation in shale or sandstone formations because paleokarst is mainly associated with the chemical dissolution of limestone. The presence of sags in the clastic material is likely due to the collapse of underlying carbonate sediments followed by infill and compaction. Hydrocarbon generation in the Woodford and Sylvan Shales may have provided the generation of sulfuric acid, which further enhanced dissolution and rock collapse (Sykes et al., 1997a) — factors that may have also created the pipe features. Therefore, we propose that the Viola Limestone and sinkholes and pipes are indicative of a mature paleokarst system.

We interpret the Wapanucka sinkholes to have formed during a period when there was subaerial exposure of the Wapanucka Limestone. Dannenberg (1952) proposes a major uplift known as the postlower Dornik Orogeny, which occurred before Atoka deposition during the final deposition of the Wapanucka Limestone. Seismic amplitude analysis given earlier suggests a hydrothermal alteration of the Wapanucka Limestone and overlying Pennsylvanian shales in sinkhole features. Therefore, we propose that Wapanucka sinkholes represent immature paleokarst, with hydrothermal rock property alteration with no measurable seismic relief (Aboaba and Liner, 2018).

We observe that the Wapanucka sinkhole features are curiously in the same spatial location as the deeper Viola sinkholes. McDonnell et al. (2007) report a similar phenomenon in the Pennsylvanian Marble Falls and Ordovician Ellenburger karst in the Fort Worth Basin. The lower section of the Marble Falls Formation of Central Texas, which is thought to be of Morrowan age and, thus, possibly comparable to the Wapanucka (Strimple and Nassichuk, 1965). Although we do not see any extension of the pipes cutting through the Wapanucka, the pipes may have induced subseismic faults or fractures serving as fluid pathways leading to alteration of the Wapanucka Limestone and overlying Pennsylvanian shale.

Visual examination of the interpreted Marble Falls and Ellenburger intervals on seismic sections in McDonnell et al. (2007) and Qi et al. (2014) reveal a significant structural low or sag, with pipe features extending beyond the Marble Falls into the Lower Atoka Runaway Formation. McDonnell et al. (2007) note that if there had been a paleokarst occurrence in the Marble Falls, it might have followed pathways developed by earlier the Ellenburger collapse, an explanation perhaps applicable to our data (Aboaba and Liner, 2018).

We have studied four seismic horizons for evidence of paleokarst; the Ordovician Viola Limestone, Mississippian Caney Shale, Pennsylvanian Jefferson Sandstone, and the Pennsylvanian Wapanucka Limestone. In areas of good seismic data quality, probable karst collapse pipe features are observed with a vertical extent up to 490 m (1610 ft). In map view, the pipes have a diameter of 150–520 m (500–700 ft) and are spaced 460–2130 m (1510–7000 ft) apart (Aboaba and Liner, 2018). The collapse pipes extend below the Viola into the Simpson Group and upward cutting across the Caney, but not extending to the top of the Jefferson. The pipes are characterized by high variance and are coincident with Viola sinkholes and Caney sags that show measurable relief. The collapse features below the Viola may actually be karsting of deeper features of the Simpson Group, or they may be velocity pushdown effects due to the decreased porosity and increased velocity of the Viola and the increased thickness infill of lower velocity Caney Shale.

A calibrated Gassmann and amplitude analysis for the Viola implies a drop in acoustic impedance corresponding to a porosity increase of approximately 10% inside sinkholes relative to adjacent rock. We interpret the Viola sinkholes and pipe features to be indicative of a mature epigene paleokarst system formed by subareal exposure and dissolution by meteoric waters.

The Wapanucka Limestone shows no measurable relief in sinkhole features that are observed on seismic amplitude. The Wapanucka sinkholes are seen on the H1 horst block and are spatially coincident with deeper Viola sinkholes. Calibrated Gassmann and seismic amplitude analysis of the Wapanucka Limestone shows that observed sinkhole amplitude cannot be reconciled with any plausible alteration of the limestone alone. We conclude that Wapanucka sinkholes represent immature hypogene paleokarst formed by limited subareal exposure and later hydrothermal alteration of the Wapanucka Limestone and overlying shale.

This study provides an interpretive framework for identifying mature and immature, epigene and hypogene paleokarst, from seismic and well data, which may be applicable to similar subsurface carbonate settings worldwide.

We thank Devon Energy for providing seismic and well data, Schlumberger and CGG for software support, and the Python library Petropy used for plotting wireline logs. This paper benefitted from an early review by J. Blackstock. We thank the anonymous reviewers whose positive feedback further improved the quality of this work.

Data associated with this research are confidential and cannot be released.

Olanrewaju Aboaba received a B.Tech. in applied geophysics from the Federal University of Technology, Akure, Nigeria, and an M.S. in geophysics from the University of Oklahoma, Norman, USA. He is currently a Ph.D. candidate in geosciences at the University of Arkansas, Fayetteville, USA. He has worked with Degeconek and BP. He is a member of SEG, AAPG, NABG, and the Nigerian Association of Petroleum Explorationists. His research interests include seismic interpretation, seismic attribute analysis, and rock physics.

Christopher Liner received a B.S. in geology from the University of Arkansas, an M.S. in geophysics from the University of Tulsa, and a Ph.D. in geophysics from the Colorado School of Mines working with the Center for Wave Phenomena. He served as SEG president (2014–2015), editor of Geophysics (1999–2001), and the 2012 SEG distinguished instructor (Elements of Seismic Dispersion). He is known to many as the author of Elements of 3D seismology (3rd edition) and The art and science of seismic interpretation (coauthored with T. L. McGilvery). He has authored many technical papers and scientific meeting abstracts and the “What’s new in exploration” column in World Oil Magazine (2010). His background is academic and industrial. Eleven years of business experience includes Western Geophysical, Conoco, Golden Geophysical, and Saudi Aramco. He has held faculty positions at the universities of Tulsa (1990–2004), Houston (2008–2012), and, currently, Arkansas. He now serves as the inaugural Maurice F. Storm Endowed Chair in Petroleum Geology in the Department of Geosciences as well as department chair, with a research interest in carbonates, advanced seismic interpretation methods, and machine learning. He is a member of SEG, AAPG, AGU, the Seismological Society of America, an honorary member of the Geophysical Society of Houston, and a corresponding member of the European Academy of Sciences.