Drilling and completion of horizontal wells are the largest expenses in unconventional reservoir plays, where the cost of drilling a well is proportional to the time it takes to reach the target objective. Accordingly, the faster the desired penetration depth and offset is achieved, the lower the cost of the drilling process. The rate of penetration (ROP) is measured in all wells, but rarely examined by geoscientists. ROP depends on many factors, but the primary factors are weight on the drill bit, drill bit rotation speed, drilling fluid-flow rate, and the characteristics of the formation being drilled (Bourgoyne et al., 1986). In this study, all wells were drilled within a 2-year period using similar drilling parameters, allowing investigation of the formation characteristics on the ROP.

Various approaches have been applied to estimate ROP. One of the main challenges for ROP estimation is the variability in the interplay between the rock and drilling speed (Farrokh et al., 2012). A “drill-off test” is a method primarily used to determine an optimum ROP for a set of conditions; however, a limitation of the drill-off test is that this process produces a static weight only valid for limited conditions during the test. The drill-off test does not work well under more complex geological conditions (King and Pinckard, 2000). Gong et al. (2007) utilized numerical simulations to investigate how rock properties affected penetration rates in Tunnel Boring Machines (TBM) and found that an increase in rock brittleness caused an increase in penetration rate. Later, a numerical model was created to model penetration rate for TBM’s by Gong and Zhao (2009), who found that an increase in compression strength decreased ROP and an increase in volumetric joint count increased ROP.

In addition to well logs and cores, seismic attributes are widely used to predict lithological and petrophysical properties of reservoirs. For example, curvature anomalies commonly indicate an increase in rock strain, which in turn can be used to infer fractures (White et al., 2012). Impedance inversion is currently the most direct seismic-based estimate of rock properties. Seismic-impedance inversion results have been used to predict fault zones, potential fractures, and lithology in the Mississippian limestone (Dowdell et al., 2013; Roy et al., 2013; Verma et al., 2013; Lindzey et al., 2015). Young’s modulus, E, and Poisson’s ratio, $υ$⁠, calculated from bulk density, $ρ$⁠, compressional velocity, V_p, and shear velocity, V_s logs can be used to estimate rock brittleness (Harris et al., 2011).

Drilling and borehole measurements such as ROP are usually not linearly related to volumetric seismic attributes, such that the use of multilinear regression is limited. Artificial neural networks (ANN) are commonly used to link attributes to properties such as gamma-ray response (Verma et al., 2013), total organic carbon (TOC; Verma et al., 2016), and well production (Da Silva et al., 2012). The proximal support vector machine (PSVM) method is a more recent innovation that has been successfully used to predict brittleness (Zhang et al., 2015). PSVM utilizes pattern recognition and classifies points by mapping them to a higher dimension before assigning them to categories. PSVM has been applied in seismic facies recognition (e.g., channels, mass-transport complexes, etc.) (Zhao et al., 2015) and lithofacies classification (Zhao et al., 2014). Zhao et al. (2014) used PSVM to categorize shale and limestone on well logs with training inputs of gamma-ray and sonic logs.

With the recent onset of unconventional techniques such as horizontal drilling and hydraulic fracturing, the Mississippian limestone has seen a growth in activities. Where operators once targeted structural traps with vertical wells, now they target stratigraphic traps with horizontal wells (Lindzey et al., 2015). Such horizontal wells require a better understanding of the variability within the Mississippian limestone to increase the success and efficiency of precisely targeted directional wells. Throughout this study, a workflow is presented to establish a relationship between seismic attributes and rock mechanical properties with ROP to optimize well placement and decrease the drilling cost.

The Mississippian limestone is a broad informal term that refers to dominantly carbonate deposits of the midcontinent (Parham and Northcutt, 1993). The main depositional environment represented in north-central Oklahoma is associated with the east–west trending ramp margin of the Burlington shelf of a starved basin environment (Costello et al., 2013). The thickness of the Mississippian limestone ranges from $350ft(106.7m)to700ft(213.4m)$ north to south over the study area (Costello et al., 2013).

Mississippian limestone in the study area was deposited in a southward prograding system near the shelf margin during the Osagean and Meramecian (Costello et al., 2013). This environment has resulted in commonly acknowledged facies within the Mississippian carbonates, ranging from shale, chert conglomerate, tripolitic chert, dense chert, altered chert-rich limestone, dense limestone, to shale-rich limestone (Lindzey et al., 2015). In the study area, tripolitic chert is most prevalent in the Upper Mississippian zones and rapidly decreases in abundance at depth greater than 150 ft (45.7 m) below the pre-Pennsylvanian unconformity (Lindzey et al., 2015).

During the Early Mississippian, warm oxygenated water covered much of the ramp in the study area. Sponge-microbe bioherms formed elongate mounds below storm wavebase and produced abundant $SiO2$ spicules, which led to formation of spicule-rich wackestones and packstones (Lindzey et al., 2015). Limestone and cherty limestone rich in marine fauna were the dominant sediments deposited at this time (Parham and Northcutt, 1993).

Regional uplift occurred during the Pennsylvanian, creating the Pennsylvanian unconformity that overlies most of the Mississippian in the midcontinent (Parham and Northcutt, 1993). The uplift not only removed large sections of rock but also reworked and caused alteration at the top of the Mississippian section and created detrital deposits of reworked Mississippian-aged rocks (Rogers, 2001). These altered sections of rocks comprise highly porous tripolitic chert and very dense glass-like chert. The leaching due to meteoric waters during relative sea-level fall has led to karstification and the formation of caverns and solution-channel features (Parham and Northcutt, 1993).

In the study area, diagenesis left intensely altered Mississippian limestone after deposition, and one of the most prominent of these diagenetic features is silica replacement (Lindzey et al., 2015). Water washed through the pores and redistributed the siliceous volcanic ash and some macrofossils, which left extensive microscale porosity (Lindzey et al., 2015). The dissolved silica precipitated in pore space and partially or completely replaced some carbonate fossils (Lindzey et al., 2015). Pore sediment structures are not well preserved due to the strong diagenetic overprint. Chert nodules are present, especially in highly reworked and bioturbated zones. Fractures are often filled with silica or calcite (Costello et al., 2013).

Molds, fractures, channels, and especially vugs are the most prominent pore type observed in the Mississippian interval of the study area (Lindzey et al., 2015). Vuggy porosity is often associated with tripolite but also exists in the other dominant facies. In many places where silica replacement took place, extensive secondary porosity formed in the shape of vugs (Rogers, 2001). Moldic porosity is also common, especially in packstone and grainstone facies that exhibit more skeletal grains. Moldic porosity develops by dissolution of sponge spicules (Montgomery et al., 1998). Both fracture and channel porosity exist but are less abundant compared to the other pore types (Lindzey et al., 2015).

In 2010, Chesapeake Energy acquired a $400mi2(1036km2)$ 3-D seismic survey in Woods County, Oklahoma (Figure 1A). The seismic processing workflow included refraction statics, velocity analysis, residual statics, prestack time migration, frequency–space–time (FXY) predictive noise rejection, and Ormsby filtering. The overall data quality is excellent. The signal-to-noise ratio (S/N) is relatively high and the wavelet amplitude appears continuous throughout the Mississippian target. The data set includes digital well logs and mud logs for 83 wells consisting of 52 horizontal and 31 vertical wells. For the ROP analysis, only horizontal wells were used. These data consisted of 52 gamma-ray logs, 51 mud logs, and 18 of them are open-hole logs.

The wells in the area of interest were drilled by the same operator in a similar time period; therefore, we assume consistency between the wells regarding weight on bit, mud type, and bit type. This study evaluates the impact of geological properties on the ROP. The work flow contains three steps: training, validation, and classification (Figure 2). Prestack inversion and seismic-attribute volumes were generated for the Mississippian limestone and converted to depth. Geomechanical rock properties (from seismic inversion) and seismic-attribute values were interpolated and then extracted along each wellbore every 2 ft (0.61 m) corresponding to the well-coordinate system from the mud logs. The mud logs give ROP in units of min/ft, which is an inverse velocity. We define the inverse velocity to be the cost of penetration (COP). The mean and standard deviation of COP for the 51 horizontal wells resulted in two categories: high and low COP with average values of $27and2.5min/ft(89and8.2min/m)$⁠, respectively (Figure 3). Each coordinate location is assigned a COP category and a set of values including seismic attributes and geomechanical rock properties. The category and values for $30%$ of the wells were used as inputs to train the model. The remaining $70%$ of the wells were used to validate the model. When an optimal accuracy is reached, the model is used to classify the entire data set where wells have yet to be drilled and no COP data are available.

Formation tops for the Lansing, Mississippian, and Woodford units were interpreted on the time-migrated seismic data in the time domain and on well logs in the depth domain. A conversion velocity model $(V0(x,y))$ was built using commercial software PETREL (© Schlumberger), where velocity, V₀, is defined at the top of the Lansing datum, $Z0(x,y)$⁠. Depth, Z, is calculated by adding the depth below the Lansing, $ΔZ=V0Z(x,y)×[t−t0(x,y)]$ to the datum. The well tops were used as a correction factor in the creation of the velocity model. Well data were assigned more weights than the seismic data. We followed the recommended settings to build the velocity model, such that a moving-average method was used as an interpolation approach for creating the new depth surfaces and an inverse-distance-squared algorithm was used to compute the inverse distance during the interpolation processes. Because the seismic horizons honored the faults in the study area, the velocity model is computed taking faults into consideration.

Geometric attributes for this data set were generated using software AASPI developed at the University of Oklahoma. The attributes generated included: most positive curvature, k₁, most negative curvature, k₂, curvedness, $k12+k22$⁠, shape index, $s=2πATAN(k2+k1k2−k1)$⁠, coherent energy, and coherence. These attributes were chosen because of their ability to delineate the structural complexity in the area of interest. The sampling interval of these attributes is the same as the original seismic data volume, $110×110ft(33.5×33.5m)$⁠. To match the mud log coordinate spacing, linear interpolation was used to generate values at $2ft(0.61m)$ intervals.

Geomechanical rock properties were derived from prestack inversion results using commercial software Hampson Russell (©CGG GeoSoftware). Data preconditioning steps, prior to a prestack seismic inversion included phase shift, bandpass filtering $(10–15—110–120Hz)$⁠, parabolic Radon transform, and trim statics.

Exploratory data analysis consisted of evaluation of two different families of volumetric attributes as input to PSVM classification: geometric attributes and geomechanical rock properties with the goal of determining which attributes are most sensitive to COP in the heterogeneous Mississippian limestone.

Geometric attributes are used to aid in the interpretation of folds and faults. Based on the TBM observation by Gong et al. (2007), we hypothesize that COP is affected by faults and fractures. Therefore, we examined the correlation of the structural attributes coherence, dip magnitude, curvature, and curvedness to our two well clusters (Figure 6). The attribute histograms indicate little to no separation for coherence and dip magnitude; however, curvature and curvedness exhibit measurable separation. Figure 4D indicates that low curvedness correlates to low COP.

TBM analysis by Gong et al. (2007) also suggested that mechanical properties play a significant role in the variation of COP. Using prestack seismic inversion we computed porosity, lambda–rho, $λρ$⁠, mu–rho, $μρ$⁠, and P-impedance volumes to analyze the Mississippian Limestone (Figure 5). The P-impedance measures the product of density and seismic P velocity. $λρ$ and $μρ$ are used to estimate lithology and geomechanical behavior such as the brittleness index (Perez-Altamar and Marfurt, 2014). Figure 4B, C, E shows the high degree of separation for these rock properties. Low COP is related to low porosity, high $λρ$⁠, high $μρ$⁠, and high P-impedance values. Conversely, high porosity, low $λρ$⁠, low $μρ$⁠, and low P-impedance values are indicative of high COP. These differences were used to train the PSVM model and classify COP data based on the geomechanical rock properties within the Mississippian interval in the study area (Figure 10).

The rectangular frame separating the dark gray circle from the light gray circle in Figure 7B is called a discriminator. Note that many of the measurements cannot be separated in Figure 7A. Because Gong and Zhao (2009) found that increased brittleness improved TBM performance, we examine brittleness as a means to predict COP. Perez-Altamar and Marfurt (2014) used geomechanical properties to predict brittleness index for shale plays in the U.S.A. We display a crossplot in Figure 9 where each sample was color-coded by COP and plotted in a 2-D space. Then we manually defined high COP (red), low COP (green), and mixed COP (yellow) polygons to define a 3-cluster template. A crossplot of $λρ$ and $μρ$llustrates the limitations of manually picking in Figure 9A illustrates the limitations of manually picking clusters in two-attribute space, where 50% of the voxels fall into the mixed COP (yellow) class. Figure 9B, a crossplot of $ρ$ and $Vp/Vs$⁠, further shows this problem with the handpicked clusters where an even larger number of voxels falls into the mixed (yellow) class. Figure 8 suggests improved class separation when using three attributes. However, interactive cluster definition in 3-D is significantly more challenging than in 2-D.

Visualization and interactive visualization with more than three attributes is intractable. PSVM addresses this problem in two ways. First, it projects the data, in this case, two attributes defining a 2-D space that cannot be separated by a linear discriminator, into a higher 3-dimensional space (Figure 7) where separation by a planar discriminator is possible. Second, because the discriminator generation is machine driven rather than interpreter driven, one can introduce more than three input attributes. We used the five attributes, curvedness, $λρ$⁠, $μρ$⁠, porosity, and P-impedance found to exhibit good histogram separation in all exploratory data analysis steps (Figure 3). The PSVM method allows us to create a classification model based on a set of training input. As the dimensionality of the input increases, the model becomes more accurate at classifying COP within the data set. For instance, during the validation process, we found the model to be sensitive to porosity. Before porosity was introduced to the model, the accuracy was $88.9%.$ When porosity was added as a new degree of dimensionality, the accuracy increased to $94%.$ This allowed for the creation of an optimal model with five degrees of dimensionality for COP classification across the study area.

A comparison of the histograms (Figure 4) shows that the generated PSVM model is more sensitive to geomechanical rock properties than geometric attributes. Indeed, strain (measured by curvature) is only one component necessary to generate natural fractures. Stearns (2015) found fractures measured in horizontal image logs through a different Mississippi limestone survey were highly correlated to gamma ray (lithology) response and only less connected to curvature, if at all. Nevertheless, this is not to say structural attributes such as curvature have no effect on the model. We observed that higher COP values are linked with higher curvedness, which indicated that it is harder to drill through the formation with higher structural complexities. Studies have found that large curvature values are related with natural fractures, which may or not be cemented (Bourgoyne et al., 1986; Hunt et al., 2011). Such heterogeneities may slow the drilling progress. Porosity is another good indicator of microstructures associated with fracture geometry. Low porosity observed in low-COP wells may seem counter-intuitive at first; however, woodworkers observe that there are few bit problems when drilling through oak, but that the bit often gets stuck or even breaks when drilling relatively “soft” pine (Neher, 1993). Again using the woodworker’s analogy, one uses different saw blades for different woods. The bits used in this survey may have been chosen to deal effectively with the very hard chert.

COP is a major factor affecting the time spent drilling a well and is directly related to the overall cost of the drilling process. This is the first study that links COP to seismic data and seismic-related attributes. Clustering five attributes using a PSVM classification method, we were able to correlate COP with seismic attributes and geomechanical rock properties and obtain a confidence of $94%$⁠. Low COP was observed in wells encountering low porosity, high $λρ$⁠, high $μρ$⁠, low curvedness, and high P-impedance. High COP was observed in wells encountering high porosity, low $λρ$⁠, low $μρ$⁠, high curvedness, and low P-impedance. By using this workflow, we can use COP of previously drilled wells with 3-D seismic data to predict COP over the study area. Whereas one may still wish to drill a specific target objective, we claim that this statistical analysis technique will provide a more accurate cost estimate and help choose the appropriate drilling equipment.

We thank the sponsors of the “Mississippi Lime” and Attribute Assisted Seismic Processing & Interpretation (AASPI) Consortia at the University of Oklahoma: Anadarko, Arcis, BGP, BHP Billiton, Lumina Geophysical, Chesapeake Energy, Chevron, ConocoPhillips, Devon Energy, EnerVest, ExxonMobil, Geophysical Insights, Institute of Petroleum, Marathon Oil, Newfield Exploration, Occidental Petroleum, Petrobras, Pioneer Natural Resources, QEP, Remark Energy Technology, Schlumberger, Shell, SM Energy, Southwestern Energy, and Tiptop Energy (Sinopec). We especially thank Chesapeake Energy for providing the 3-D seismic and well data. Petrel 2015 software was graciously provided by Schlumberger. Seismic attributes and the proximal Support Vector Machine (PSVM) results were generated using AASPI software. We would like to acknowledge CGG GeoSoftware for its donation of its software Hampson Russell, that was used for seismic inversion analysis and porosity estimation. We would like to thank our colleague Mr. Abdulmohsen Al Ali for generating seismic prestack inversion.