Drilling Analytics and Optimization

schedule15 min readfitness_center10 exercises

Why This Chapter Exists

Drilling a well is one of the most expensive operations in the petroleum industry. An onshore well might cost $5 million. A deepwater well can exceed$ 150 million. Every day on the rig costs between $100,000 and$ 1,000,000 depending on the environment. When the rate of penetration drops, when the bit stalls, when the drillstring vibrates in ways it should not, the cost clock keeps ticking.

Drilling generates enormous volumes of real-time data: weight on bit, rotary speed, torque, rate of penetration, mud properties, pressures at multiple points in the system, vibration measurements, and more --- all recorded every few seconds. The problem is that most of this data is watched in real time by a driller staring at gauges, then archived and rarely analysed systematically.

Python changes that. This chapter teaches you to ingest drilling data, compute the engineering metrics that reveal inefficiency, detect common drilling dysfunctions from data patterns, and build models that predict rate of penetration. The goal is to turn drilling data from an archive into a decision-making tool.

infoWhat You Will Learn

Load, clean, and structure real-time drilling data
Calculate rate of penetration and drilling efficiency metrics
Compute Mechanical Specific Energy and interpret what it reveals
Detect drilling dysfunctions from data signatures
Build ROP prediction models using engineering features
Perform well time and cost analysis

Drilling Data

Drilling data arrives in time-indexed records, typically at intervals of one to ten seconds during active drilling. The primary measurements are:

Bit Depth / Hole Depth --- how deep the well currently is (feet or metres)
Rate of Penetration (ROP) --- how fast the bit is advancing through rock (ft/hr)
Weight on Bit (WOB) --- the force applied to the bit face (thousands of pounds, or klbf)
Rotary Speed (RPM) --- how fast the drillstring is turning (revolutions per minute)
Torque (TQ) --- the rotational force on the drillstring (ft-lbs)
Standpipe Pressure (SPP) --- the pressure of the mud being pumped downhole (psi)
Mud Weight In / Out --- the density of mud going in and coming back (ppg)
Flow Rate --- volume of mud being circulated (gal/min)

Not all of this data is useful at all times. When the driller is tripping pipe (pulling the drillstring out of the hole or running it back in), the ROP is zero and WOB is meaningless. When circulating without drilling, there is no rock being cut. Filtering the data to active drilling intervals is the first and most important preprocessing step.

main.py

import numpy as np
import pandas as pd

# Generate realistic drilling data for a 10,000 ft well
np.random.seed(42)
n_points = 5000

# Depth increases with variable ROP
depth_start = 500
rop_base = np.piecewise(
    np.linspace(0, 1, n_points),
    [np.linspace(0, 1, n_points) < 0.3,
     (np.linspace(0, 1, n_points) >= 0.3) & (np.linspace(0, 1, n_points) < 0.6),
     np.linspace(0, 1, n_points) >= 0.6],
    [lambda x: 120 - 80*x,        # soft formation, fast drilling
     lambda x: 45 + 10*np.sin(20*x),  # medium formation, variable
     lambda x: 25 - 10*(x-0.6)]    # hard formation, slow
)
rop = rop_base + np.random.normal(0, 5, n_points)
rop = np.clip(rop, 5, 200)

# Cumulative depth from ROP (assuming 2-minute measurement intervals)
time_interval_hr = 2 / 60  # 2 minutes in hours
depth = depth_start + np.cumsum(rop * time_interval_hr)

# Drilling parameters correlated with formation
wob = np.where(depth < 4000, 15 + np.random.normal(0, 2, n_points),
       np.where(depth < 7000, 25 + np.random.normal(0, 3, n_points),
                35 + np.random.normal(0, 3, n_points)))
wob = np.clip(wob, 5, 50)

rpm = np.where(depth < 4000, 140 + np.random.normal(0, 10, n_points),
       np.where(depth < 7000, 120 + np.random.normal(0, 8, n_points),
                100 + np.random.normal(0, 8, n_points)))
rpm = np.clip(rpm, 40, 180)

torque = 0.8 * wob * (depth / 1000) + np.random.normal(0, 500, n_points)
torque = np.clip(torque, 500, 30000)

spp = 2000 + 0.15 * depth + np.random.normal(0, 100, n_points)

# Inject a dysfunction zone (stick-slip vibration between 6000-6500 ft)
stick_slip_mask = (depth > 6000) & (depth < 6500)
rpm[stick_slip_mask] = rpm[stick_slip_mask] + 40 * np.sin(np.linspace(0, 50, stick_slip_mask.sum()))
torque[stick_slip_mask] = torque[stick_slip_mask] * 1.5

drilling_data = pd.DataFrame({
    "Depth (ft)": np.round(depth, 1),
    "ROP (ft/hr)": np.round(rop, 1),
    "WOB (klbf)": np.round(wob, 1),
    "RPM": np.round(rpm, 0).astype(int),
    "Torque (ft-lbs)": np.round(torque, 0).astype(int),
    "SPP (psi)": np.round(spp, 0).astype(int),
})

print("DRILLING DATA SUMMARY")
print("=" * 55)
print(f"  Records:       {len(drilling_data):,}")
print(f"  Depth range:   {drilling_data['Depth (ft)'].min():,.0f} – {drilling_data['Depth (ft)'].max():,.0f} ft")
print(f"  Avg ROP:       {drilling_data['ROP (ft/hr)'].mean():.1f} ft/hr")
print(f"  Avg WOB:       {drilling_data['WOB (klbf)'].mean():.1f} klbf")
print(f"  Avg RPM:       {drilling_data['RPM'].mean():.0f}")
print()
print(drilling_data.describe().round(1).to_string())

Mechanical Specific Energy

Mechanical Specific Energy (MSE) is the single most important drilling efficiency metric. It represents the energy required to destroy a unit volume of rock:

MSE = \frac{480 \times T \times RPM}{D_b^2 \times ROP} + \frac{4 \times WOB}{\pi \times D_b^2}

where:

$T$ = torque (ft-lbs)
$RPM$ = rotary speed (rev/min)
$D_b$ = bit diameter (inches)
$ROP$ = rate of penetration (ft/hr)
$WOB$ = weight on bit (lbf)

MSE has units of psi (energy per unit volume of rock). In a perfectly efficient system, MSE would equal the rock's confined compressive strength (CCS). In practice, MSE is always higher than CCS because energy is wasted to friction, vibration, and inefficient cutting.

The ratio of CCS to MSE is the mechanical efficiency. An efficiency of 30–40% is typical. When efficiency drops below 20%, something is wrong --- the bit may be dull, the well may have vibration problems, or the drilling parameters may be poorly chosen.

main.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def compute_mse(torque_ftlbs, rpm, wob_lbf, rop_fthr, bit_diameter_in):
    """
    Calculate Mechanical Specific Energy.

    Parameters
    ----------
    torque_ftlbs : array-like
        Torque (ft-lbs).
    rpm : array-like
        Rotary speed (rev/min).
    wob_lbf : array-like
        Weight on bit (lbf). Note: input in lbf, not klbf.
    rop_fthr : array-like
        Rate of penetration (ft/hr).
    bit_diameter_in : float
        Bit diameter (inches).

    Returns
    -------
    numpy.ndarray
        MSE values (psi).
    """
    torque = np.asarray(torque_ftlbs, dtype=float)
    r = np.asarray(rpm, dtype=float)
    w = np.asarray(wob_lbf, dtype=float)
    rop = np.asarray(rop_fthr, dtype=float)
    db = bit_diameter_in

    # Avoid division by zero
    rop_safe = np.where(rop > 0, rop, 1e-6)

    rotary_term = 480 * torque * r / (db**2 * rop_safe)
    wob_term = 4 * w / (np.pi * db**2)

    return rotary_term + wob_term


# Calculate MSE for the drilling data
bit_diameter = 8.5  # inches
drilling_data["MSE (psi)"] = compute_mse(
    drilling_data["Torque (ft-lbs)"],
    drilling_data["RPM"],
    drilling_data["WOB (klbf)"] * 1000,  # convert klbf to lbf
    drilling_data["ROP (ft/hr)"],
    bit_diameter
)

# Estimated rock strength for efficiency calculation
# (simplified: increases with depth)
ccs_estimate = 5000 + 1.5 * drilling_data["Depth (ft)"]
drilling_data["Efficiency (%)"] = np.round(ccs_estimate / drilling_data["MSE (psi)"] * 100, 1)
drilling_data["Efficiency (%)"] = drilling_data["Efficiency (%)"].clip(0, 100)

# Plot MSE vs depth
fig, axes = plt.subplots(1, 3, figsize=(14, 8), sharey=True)

# ROP
axes[0].plot(drilling_data["ROP (ft/hr)"], drilling_data["Depth (ft)"], 'k-', linewidth=0.3, alpha=0.5)
axes[0].set_xlabel("ROP (ft/hr)", fontsize=11)
axes[0].set_ylabel("Depth (ft)", fontsize=11)
axes[0].set_title("Rate of Penetration", fontsize=12)
axes[0].invert_yaxis()
axes[0].grid(True, alpha=0.3)

# MSE
mse_clipped = drilling_data["MSE (psi)"].clip(0, 150000)
axes[1].plot(mse_clipped, drilling_data["Depth (ft)"], 'b-', linewidth=0.3, alpha=0.5)
axes[1].set_xlabel("MSE (psi)", fontsize=11)
axes[1].set_title("Mechanical Specific Energy", fontsize=12)
axes[1].grid(True, alpha=0.3)

# Efficiency
axes[2].plot(drilling_data["Efficiency (%)"], drilling_data["Depth (ft)"], 'g-', linewidth=0.3, alpha=0.5)
axes[2].axvline(x=30, color='r', linestyle='--', linewidth=1, label="30% threshold")
axes[2].set_xlabel("Mechanical Efficiency (%)", fontsize=11)
axes[2].set_title("Drilling Efficiency", fontsize=12)
axes[2].legend(fontsize=9)
axes[2].grid(True, alpha=0.3)

plt.suptitle("Drilling Performance vs. Depth", fontsize=13, y=1.01)
plt.tight_layout()
plt.show()

# Summary statistics by depth interval
drilling_data["Depth Bin"] = pd.cut(drilling_data["Depth (ft)"], bins=10)
summary = drilling_data.groupby("Depth Bin", observed=True).agg({
    "ROP (ft/hr)": "mean",
    "WOB (klbf)": "mean",
    "RPM": "mean",
    "MSE (psi)": "mean",
    "Efficiency (%)": "mean",
}).round(1)

print("\nDRILLING PERFORMANCE BY DEPTH INTERVAL")
print(summary.to_string())

Drilling Dysfunction Detection

Drilling dysfunctions waste time, damage equipment, and increase cost. The most common dysfunctions leave distinct signatures in the drilling data:

Stick-slip vibration occurs when the bit alternates between sticking (zero RPM at the bit) and spinning free (very high RPM). The surface RPM appears normal, but downhole the bit experiences extreme torsional oscillation. In the data, this manifests as high torque variability and inconsistent ROP at constant surface RPM.

Bit bounce is axial vibration where the bit repeatedly lifts off and impacts the rock. It causes low ROP relative to WOB and erratic weight measurements.

Whirl is lateral vibration where the bit orbits the borehole rather than rotating concentrically. It causes an enlarged, out-of-gauge hole and rapid bit wear.

main.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def detect_stick_slip(torque, rpm, window=50, torque_cv_threshold=0.3):
    """
    Detect stick-slip vibration from surface measurements.

    Stick-slip is indicated by high coefficient of variation
    in torque at relatively stable surface RPM.

    Parameters
    ----------
    torque : array-like
        Torque measurements (ft-lbs).
    rpm : array-like
        Surface RPM measurements.
    window : int
        Rolling window size.
    torque_cv_threshold : float
        Coefficient of variation threshold for flagging.

    Returns
    -------
    pandas.Series
        Boolean mask where True indicates stick-slip.
    """
    tq = pd.Series(torque)
    r = pd.Series(rpm)

    tq_mean = tq.rolling(window).mean()
    tq_std = tq.rolling(window).std()
    rpm_std = r.rolling(window).std()

    # Coefficient of variation in torque
    tq_cv = tq_std / tq_mean.replace(0, np.nan)

    # RPM should be relatively stable (driller holding steady)
    rpm_stable = rpm_std < 15

    return (tq_cv > torque_cv_threshold) & rpm_stable


# Apply detection
drilling_data["Stick-Slip Flag"] = detect_stick_slip(
    drilling_data["Torque (ft-lbs)"],
    drilling_data["RPM"]
)

ss_pct = drilling_data["Stick-Slip Flag"].mean() * 100
ss_depth_min = drilling_data[drilling_data["Stick-Slip Flag"]]["Depth (ft)"].min()
ss_depth_max = drilling_data[drilling_data["Stick-Slip Flag"]]["Depth (ft)"].max()

print("DYSFUNCTION DETECTION REPORT")
print("=" * 50)
print(f"  Stick-slip occurrences: {ss_pct:.1f}% of drilling time")
if ss_pct > 0:
    print(f"  Affected depth range:   {ss_depth_min:,.0f} – {ss_depth_max:,.0f} ft")
    avg_rop_ss = drilling_data[drilling_data["Stick-Slip Flag"]]["ROP (ft/hr)"].mean()
    avg_rop_normal = drilling_data[~drilling_data["Stick-Slip Flag"]]["ROP (ft/hr)"].mean()
    print(f"  Avg ROP during stick-slip: {avg_rop_ss:.1f} ft/hr")
    print(f"  Avg ROP normal drilling:   {avg_rop_normal:.1f} ft/hr")
    print(f"  ROP reduction:             {(1 - avg_rop_ss/avg_rop_normal)*100:.1f}%")

# Plot torque and RPM in the affected zone
mask = (drilling_data["Depth (ft)"] > 5500) & (drilling_data["Depth (ft)"] < 7000)
zone = drilling_data[mask].reset_index(drop=True)

fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)

axes[0].plot(zone.index, zone["RPM"], 'k-', linewidth=0.5)
axes[0].set_ylabel("RPM", fontsize=10)
axes[0].set_title("Drilling Parameters — Stick-Slip Zone (5,500–7,000 ft)", fontsize=12)
axes[0].grid(True, alpha=0.3)

axes[1].plot(zone.index, zone["Torque (ft-lbs)"], 'b-', linewidth=0.5)
axes[1].set_ylabel("Torque (ft-lbs)", fontsize=10)
axes[1].grid(True, alpha=0.3)

# Highlight stick-slip flags
ss_zone = zone[zone["Stick-Slip Flag"]]
axes[2].plot(zone.index, zone["ROP (ft/hr)"], 'k-', linewidth=0.5, label="ROP")
axes[2].scatter(ss_zone.index, ss_zone["ROP (ft/hr)"], c='red', s=5, zorder=5, label="Stick-Slip")
axes[2].set_ylabel("ROP (ft/hr)", fontsize=10)
axes[2].set_xlabel("Sample Index", fontsize=10)
axes[2].legend(fontsize=9)
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

ROP Prediction

Rate of penetration depends on controllable parameters (WOB, RPM, flow rate, mud weight) and uncontrollable parameters (rock strength, formation type, pore pressure). A prediction model that captures these relationships allows the drilling engineer to choose the controllable parameters that maximize ROP for the given formation.

main.py

import numpy as np
import pandas as pd
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, r2_score
import matplotlib.pyplot as plt

# Features for ROP prediction
features = ["WOB (klbf)", "RPM", "Torque (ft-lbs)", "SPP (psi)", "Depth (ft)"]
target = "ROP (ft/hr)"

X = drilling_data[features].values
y = drilling_data[target].values

# Train/test split — use depth-based split (not random)
# to simulate predicting ROP in a new formation
split_depth = 7000
train_mask = drilling_data["Depth (ft)"] < split_depth
test_mask = drilling_data["Depth (ft)"] >= split_depth

X_train, X_test = X[train_mask], X[test_mask]
y_train, y_test = y[train_mask], y[test_mask]

print(f"Training samples: {len(X_train):,} (above {split_depth:,} ft)")
print(f"Testing samples:  {len(X_test):,} (below {split_depth:,} ft)")

# Train model
model = GradientBoostingRegressor(
    n_estimators=200,
    max_depth=5,
    learning_rate=0.1,
    random_state=42
)
model.fit(X_train, y_train)

# Evaluate
y_pred_train = model.predict(X_train)
y_pred_test = model.predict(X_test)

print(f"\nTraining MAE:  {mean_absolute_error(y_train, y_pred_train):.2f} ft/hr")
print(f"Training R²:   {r2_score(y_train, y_pred_train):.4f}")
print(f"Testing MAE:   {mean_absolute_error(y_test, y_pred_test):.2f} ft/hr")
print(f"Testing R²:    {r2_score(y_test, y_pred_test):.4f}")

# Feature importance
importance = pd.DataFrame({
    "Feature": features,
    "Importance": model.feature_importances_
}).sort_values("Importance", ascending=True)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Predicted vs actual
axes[0].scatter(y_test, y_pred_test, s=3, alpha=0.3, color='#2E86AB')
axes[0].plot([0, 100], [0, 100], 'k--', linewidth=1)
axes[0].set_xlabel("Actual ROP (ft/hr)", fontsize=11)
axes[0].set_ylabel("Predicted ROP (ft/hr)", fontsize=11)
axes[0].set_title("ROP Prediction — Test Set", fontsize=12)
axes[0].grid(True, alpha=0.3)

# Feature importance
axes[1].barh(importance["Feature"], importance["Importance"], color='#2E86AB')
axes[1].set_xlabel("Importance", fontsize=11)
axes[1].set_title("Feature Importance", fontsize=12)
axes[1].grid(True, axis='x', alpha=0.3)

plt.tight_layout()
plt.show()

Well Time and Cost Analysis

Every drilling operation is ultimately measured against two numbers: time and cost. A time-depth curve plots cumulative depth against elapsed time and is the standard way to compare drilling performance between wells, between rigs, or against the plan.

main.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Simulate time-depth data for two wells
# Well A: actual performance
# Well B: offset well (benchmark)

# Well A — includes a 2-day stuck pipe event at 6,500 ft
depths_a = np.array([0, 500, 1500, 3000, 4500, 5500, 6500, 6500, 6500, 7500, 8500, 9200])
days_a = np.array([0, 0.5, 1.5, 3, 5.5, 8, 11, 12, 13, 16, 20, 24])
cost_a = days_a * 280000  # $280k/day spread rate

# Well B — smoother, no NPT
depths_b = np.array([0, 500, 1500, 3000, 4500, 5500, 6500, 7500, 8500, 9200])
days_b = np.array([0, 0.5, 1.5, 3, 5, 7.5, 10, 13, 16.5, 20])
cost_b = days_b * 280000

# AFE (Authorization for Expenditure) — the budget
depths_plan = np.array([0, 500, 1500, 3000, 4500, 5500, 6500, 7500, 8500, 9200])
days_plan = np.array([0, 0.5, 1.5, 2.5, 4.5, 6.5, 9, 12, 15, 18])

fig, axes = plt.subplots(1, 2, figsize=(13, 7))

# Time-Depth plot
axes[0].plot(days_a, depths_a, 'r-o', markersize=4, linewidth=2, label="Well A (actual)")
axes[0].plot(days_b, depths_b, 'b-s', markersize=4, linewidth=2, label="Well B (offset)")
axes[0].plot(days_plan, depths_plan, 'k--', linewidth=1.5, label="AFE Plan")
# Annotate stuck pipe
axes[0].annotate("Stuck pipe\n(2 days NPT)", xy=(12, 6500), xytext=(16, 5500),
                 fontsize=9, color='red',
                 arrowprops=dict(arrowstyle='->', color='red', lw=1.2))
axes[0].set_xlabel("Days", fontsize=11)
axes[0].set_ylabel("Depth (ft)", fontsize=11)
axes[0].set_title("Time-Depth Curve", fontsize=12)
axes[0].invert_yaxis()
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)

# Cost curve
axes[1].plot(depths_a, cost_a / 1e6, 'r-o', markersize=4, linewidth=2, label="Well A")
axes[1].plot(depths_b, cost_b / 1e6, 'b-s', markersize=4, linewidth=2, label="Well B")
axes[1].plot(depths_plan, days_plan * 280000 / 1e6, 'k--', linewidth=1.5, label="AFE Budget")
axes[1].set_xlabel("Depth (ft)", fontsize=11)
axes[1].set_ylabel("Cumulative Cost ($M)", fontsize=11)
axes[1].set_title("Cost vs. Depth", fontsize=12)
axes[1].legend(fontsize=9)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Cost summary
print("DRILLING COST SUMMARY")
print("=" * 45)
print(f"  Well A total days:  {days_a[-1]:.0f}")
print(f"  Well A total cost:  ${cost_a[-1]/1e6:.1f}M")
print(f"  Well B total days:  {days_b[-1]:.0f}")
print(f"  Well B total cost:  ${cost_b[-1]/1e6:.1f}M")
print(f"  AFE budget:         ${days_plan[-1] * 0.28:.1f}M")
print(f"  Well A overrun:     ${(cost_a[-1] - days_plan[-1]*280000)/1e6:+.1f}M")
print(f"  NPT cost (stuck):  ${2 * 280000 / 1e6:.1f}M")

Exercises

fitness_center

Exercise 14.1Practice

-- Drilling Data Cleanup

Load a drilling dataset and filter it to active drilling intervals only (ROP > 0, WOB > 5 klbf, RPM > 30). What percentage of the raw data represents ...

arrow_forward