Part II: Petroleum Data Engineering

Chapter 7

Well Log Analysis with Python

schedule15 min readfitness_center10 exercises

Why This Chapter Exists

A well log is the closest thing petroleum engineers have to seeing underground. When a well is drilled, logging tools are lowered into the borehole and pulled back to surface, measuring the physical properties of the rock at every half-foot of depth. The resulting curves — gamma ray, resistivity, density, neutron porosity — encode the answers to the three questions that drive every reservoir decision:

What kind of rock is down there? (Sand, shale, limestone, or something else.)
Does it contain hydrocarbons? (Oil, gas, or just water.)
How much can we produce? (Porosity, permeability, net pay thickness.)

These questions determine whether a company spends $80 million completing a well or walks away. They determine reserve bookings worth billions. They shape field development plans that last decades.

A petrophysicist answers these questions by interpreting log data — applying physics-based equations to convert raw measurements into rock and fluid properties. This chapter teaches you to do that interpretation in Python, from raw LAS file to a complete petrophysical evaluation with net pay determination.

infoWhat You'll Learn

Understand what each logging tool measures physically and why it matters
Calculate shale volume from gamma ray logs
Compute porosity from density and neutron measurements
Determine water saturation using Archie's equation
Apply cutoffs to identify net pay zones
Build a complete petrophysical evaluation pipeline

What Logging Tools Actually Measure

Before writing any equations, you need to understand what each tool is responding to. This is not optional knowledge — if you do not understand the physics behind the measurement, you cannot interpret the results correctly, and no amount of Python will fix that.

Gamma Ray (GR)

Every rock contains trace amounts of radioactive elements — potassium, thorium, and uranium. Shales contain significantly more of these elements than sandstones or limestones. The gamma ray tool measures the natural radioactivity of the formation.

Why it matters: The gamma ray log is the primary tool for distinguishing reservoir rock (sand, with low GR) from non-reservoir rock (shale, with high GR). This distinction is the first step in every petrophysical analysis. If you cannot separate sand from shale, nothing else you calculate is reliable.

Typical values:

Clean sand: 15–45 GAPI
Shaly sand: 45–90 GAPI
Shale: 90–150+ GAPI

Deep Resistivity (RT)

Electrical resistivity measures how strongly the rock resists the flow of electric current. This depends on what fluid is in the pore space.

Why it matters: Water conducts electricity. Oil and gas do not. A rock saturated with salt water has low resistivity. The same rock saturated with oil has high resistivity. By measuring resistivity and knowing the porosity, you can calculate what fraction of the pore space contains water versus hydrocarbons. That fraction — the water saturation — is the single most important number in reservoir evaluation.

Typical values:

Shale: 1–10 ohm·m
Water-bearing sand: 2–20 ohm·m
Oil-bearing sand: 20–2000+ ohm·m

Bulk Density (RHOB)

The density tool emits gamma rays into the formation and measures how many return. Dense rock absorbs more gamma rays, so fewer return. The measurement gives the bulk density of the rock — the combined density of the rock matrix and the fluids in its pores.

Why it matters: If you know the density of the rock matrix (about 2.65 g/cc for quartz sandstone) and the density of the fluid in the pores (about 1.0 g/cc for water, 0.85 g/cc for oil), you can calculate the porosity — the fraction of the rock that is pore space. Porosity determines how much fluid the rock can hold.

Typical values:

Pure sandstone (zero porosity): 2.65 g/cc
Porous sandstone (20% porosity): 2.32 g/cc
Shale: 2.40–2.70 g/cc
Limestone: 2.71 g/cc

Neutron Porosity (NPHI)

The neutron tool emits fast neutrons into the formation. Hydrogen atoms slow neutrons down most effectively, and hydrogen is primarily found in water and hydrocarbons occupying the pore space. The tool measures the concentration of hydrogen, which is approximately proportional to porosity.

Why it matters: Neutron porosity provides a second, independent estimate of porosity. When density porosity and neutron porosity agree, you have confidence in the result. When they disagree, that disagreement itself carries information — it can indicate the presence of gas (which has less hydrogen per unit volume than oil or water) or the presence of shale (which contains bound water in its clay structure).

Typical values:

Clean sand with water: porosity matches true porosity
Gas-bearing sand: neutron porosity reads lower than density porosity (the "gas crossover effect")
Shale: neutron porosity reads higher than true porosity due to bound water in clays

Setting Up the Data

We will use a synthetic but realistic well log dataset that demonstrates the key features a petrophysicist looks for: clean sand zones, shale intervals, and a transition from water-bearing to oil-bearing reservoir.

main.py

import numpy as np
import pandas as pd

np.random.seed(101)

# Build a synthetic well log that tells a geological story.
# Depths from 7000 to 8000 ft, half-foot sampling.
depth = np.arange(7000, 8000.5, 0.5)
n = len(depth)

# Define lithology zones (this is what we're trying to recover from the logs)
# 7000-7200: shale
# 7200-7400: water-bearing sand
# 7400-7600: oil-bearing sand (the pay zone)
# 7600-7700: shale
# 7700-7850: water-bearing sand
# 7850-8000: shale

gr = np.full(n, 110.0)    # start with shale baseline
rt = np.full(n, 3.0)
rhob = np.full(n, 2.50)
nphi = np.full(n, 0.30)

def set_zone(depth_arr, d_start, d_end, gr_val, rt_val, rhob_val, nphi_val,
             gr_arr, rt_arr, rhob_arr, nphi_arr, noise=True):
    mask = (depth_arr >= d_start) & (depth_arr < d_end)
    n_pts = mask.sum()
    gr_arr[mask] = gr_val + (np.random.normal(0, 5, n_pts) if noise else 0)
    rt_arr[mask] = rt_val + (np.random.normal(0, rt_val*0.1, n_pts) if noise else 0)
    rhob_arr[mask] = rhob_val + (np.random.normal(0, 0.02, n_pts) if noise else 0)
    nphi_arr[mask] = nphi_val + (np.random.normal(0, 0.01, n_pts) if noise else 0)

# Shale zones
set_zone(depth, 7000, 7200, 115, 2.5, 2.52, 0.32, gr, rt, rhob, nphi)
set_zone(depth, 7600, 7700, 120, 2.0, 2.55, 0.34, gr, rt, rhob, nphi)
set_zone(depth, 7850, 8000, 110, 3.0, 2.50, 0.30, gr, rt, rhob, nphi)

# Water-bearing sand
set_zone(depth, 7200, 7400, 35, 8.0, 2.28, 0.18, gr, rt, rhob, nphi)
set_zone(depth, 7700, 7850, 30, 10.0, 2.25, 0.17, gr, rt, rhob, nphi)

# Oil-bearing sand (THE pay zone)
set_zone(depth, 7400, 7600, 28, 85.0, 2.20, 0.15, gr, rt, rhob, nphi)

# Clip to physical bounds
gr = np.clip(gr, 0, 200)
rt = np.clip(rt, 0.1, 5000)
rhob = np.clip(rhob, 1.5, 3.0)
nphi = np.clip(nphi, -0.05, 0.60)

logs = pd.DataFrame({
    "DEPTH": depth,
    "GR": np.round(gr, 2),
    "RT": np.round(rt, 2),
    "RHOB": np.round(rhob, 3),
    "NPHI": np.round(nphi, 4)
})

print(f"Well log: {len(logs)} samples, {logs['DEPTH'].min():.0f} to {logs['DEPTH'].max():.0f} ft")
print()
print(logs.describe().round(3).to_string())

Step 1: Shale Volume (Vshale)

The first calculation in any petrophysical analysis is determining how much shale is present at each depth. Shale degrades everything — it reduces effective porosity, traps water in clay minerals, and lowers permeability. You need to quantify it before you can reliably calculate anything else.

The simplest and most common approach is the linear gamma ray method:

V_{shale} = \frac{GR - GR_{clean}}{GR_{shale} - GR_{clean}}

where:

$GR$ = the gamma ray reading at a given depth
$GR_{clean}$ = the gamma ray reading in the cleanest sand (the lowest GR value in a known sand zone)
$GR_{shale}$ = the gamma ray reading in a pure shale (the highest GR value in a known shale zone)

This gives a value between 0 (pure sand) and 1 (pure shale). Values above 1 or below 0 are possible due to noise and should be clipped.

main.py

# Determine GR endpoints from the data.
# In practice, you would pick these from a histogram or
# from known clean sand and shale intervals.

gr_clean = logs["GR"].quantile(0.05)   # 5th percentile ≈ cleanest sand
gr_shale = logs["GR"].quantile(0.95)   # 95th percentile ≈ pure shale

print(f"GR endpoints:")
print(f"  GR_clean (sand baseline): {gr_clean:.1f} GAPI")
print(f"  GR_shale (shale baseline): {gr_shale:.1f} GAPI")

# Calculate linear Vshale
logs["VSHALE"] = (logs["GR"] - gr_clean) / (gr_shale - gr_clean)
logs["VSHALE"] = logs["VSHALE"].clip(0, 1)

# Summary by zone
print()
print("Average Vshale by depth zone:")
for label, d1, d2 in [("7000-7200 (shale)", 7000, 7200),
                       ("7200-7400 (wet sand)", 7200, 7400),
                       ("7400-7600 (pay sand)", 7400, 7600),
                       ("7600-7700 (shale)", 7600, 7700),
                       ("7700-7850 (wet sand)", 7700, 7850)]:
    zone = logs[(logs["DEPTH"] >= d1) & (logs["DEPTH"] < d2)]
    print(f"  {label}: Vshale = {zone['VSHALE'].mean():.3f}")

The shale volumes confirm what we expect: shale zones near 1.0, sand zones near 0.0, with some scatter from measurement noise. If these numbers did not make physical sense — for example, if a known sand zone showed a Vshale of 0.6 — you would revisit your GR endpoints or consider using a non-linear Vshale model (Larionov, Clavier, or Steiber).

Step 2: Porosity

Porosity — the fraction of rock volume that is pore space — determines how much fluid a reservoir can hold. Two independent estimates come from the density and neutron logs.

Density Porosity

\phi_{density} = \frac{\rho_{matrix} - \rho_{bulk}}{\rho_{matrix} - \rho_{fluid}}

For a quartz sandstone with water in the pores:

$\rho_{matrix}$ = 2.65 g/cc (quartz grain density)
$\rho_{fluid}$ = 1.0 g/cc (formation water density)

Effective Porosity

Raw porosity calculated from logs includes the pore space occupied by bound water in shale. Since that water is immobile and cannot be produced, it does not contribute to producible reserves. Effective porosity corrects for this:

\phi_{effective} = \phi_{total} \times (1 - V_{shale})

main.py

# Matrix and fluid parameters
rho_matrix = 2.65   # g/cc, quartz sandstone
rho_fluid = 1.00    # g/cc, formation water

# Density porosity
logs["PHID"] = (rho_matrix - logs["RHOB"]) / (rho_matrix - rho_fluid)
logs["PHID"] = logs["PHID"].clip(0, 0.45)

# Neutron porosity is already in the NPHI column (fractional units)

# Average of density and neutron porosity (common practice in clean sands)
logs["PHIA"] = (logs["PHID"] + logs["NPHI"]) / 2

# Effective porosity (corrected for shale)
logs["PHIE"] = logs["PHIA"] * (1 - logs["VSHALE"])
logs["PHIE"] = logs["PHIE"].clip(0, 0.40)

print("Porosity summary by zone:")
for label, d1, d2 in [("7200-7400 (wet sand)", 7200, 7400),
                       ("7400-7600 (pay sand)", 7400, 7600),
                       ("7700-7850 (wet sand)", 7700, 7850)]:
    zone = logs[(logs["DEPTH"] >= d1) & (logs["DEPTH"] < d2)]
    print(f"  {label}:")
    print(f"    Density porosity:    {zone['PHID'].mean():.3f}")
    print(f"    Neutron porosity:    {zone['NPHI'].mean():.3f}")
    print(f"    Effective porosity:  {zone['PHIE'].mean():.3f}")

In the pay zone (7400–7600 ft), effective porosity should be in the range of 0.15–0.25 for a producible sandstone reservoir. Values below 0.08 generally indicate rock too tight to produce economically. Values above 0.35 are uncommon in consolidated sandstones and may indicate a data or calculation error.

Step 3: Water Saturation (Archie's Equation)

This is the critical calculation. Water saturation ( $S_w$ ) tells you what fraction of the pore space contains water. The remainder ( $1 - S_w$ ) is hydrocarbons.

Archie's equation, published by Gus Archie in 1942, relates water saturation to resistivity and porosity:

S_w = \left( \frac{a \cdot R_w}{\phi^m \cdot R_t} \right)^{1/n}

where:

$S_w$ = water saturation (fraction)
$R_w$ = resistivity of formation water (ohm·m) — measured from water samples or calculated from SP logs
$R_t$ = true formation resistivity from the deep resistivity log (ohm·m)
$\phi$ = porosity (fraction)
$a$ = tortuosity factor (typically 0.62–1.0)
$m$ = cementation exponent (typically 1.7–2.2)
$n$ = saturation exponent (typically 1.8–2.2)

Why this equation works: In a water-saturated rock, the only conductive pathway for electricity is through the brine in the pore space. As oil replaces water, the conductive pathways become fewer and more tortuous, and resistivity increases. Archie's equation quantifies this relationship.

Why the parameters matter: The exponents $m$ and $n$ are not universal constants — they depend on the rock type, the pore geometry, and the wetting characteristics of the formation. Using incorrect values produces wrong water saturations, which produces wrong reserve estimates. In practice, these parameters are calibrated against core data when available.

main.py

# Archie parameters — these would normally come from
# core analysis or regional knowledge
a = 0.81       # tortuosity factor (Humble formula for sandstone)
m = 2.0        # cementation exponent
n = 2.0        # saturation exponent
Rw = 0.04      # formation water resistivity at reservoir temperature (ohm·m)

# Calculate water saturation using Archie's equation
# Only calculate where porosity is meaningful (> 0.01)
logs["SW"] = np.nan
valid = logs["PHIE"] > 0.01

logs.loc[valid, "SW"] = (
    (a * Rw) / (logs.loc[valid, "PHIE"]**m * logs.loc[valid, "RT"])
) ** (1/n)

# Clip to physical bounds
logs["SW"] = logs["SW"].clip(0, 1)

# Hydrocarbon saturation
logs["SHC"] = 1 - logs["SW"]

print("Water saturation by zone:")
for label, d1, d2 in [("7200-7400 (wet sand)", 7200, 7400),
                       ("7400-7600 (pay sand)", 7400, 7600),
                       ("7700-7850 (wet sand)", 7700, 7850)]:
    zone = logs[(logs["DEPTH"] >= d1) & (logs["DEPTH"] < d2)]
    sw_mean = zone["SW"].mean()
    shc_mean = zone["SHC"].mean()
    print(f"  {label}:")
    print(f"    Sw = {sw_mean:.3f}  (water)")
    print(f"    Shc = {shc_mean:.3f}  (hydrocarbons)")

The results should show what the resistivity log already told us: the wet sand zones have high water saturation (Sw near 1.0), and the pay zone has low water saturation (Sw near 0.2–0.4), meaning 60–80% of the pore space is filled with oil.

If the water saturation in a "known" oil zone comes out above 0.6, something is wrong — either the Archie parameters are incorrect, the porosity calculation has errors, or the zone is not actually oil-bearing. This is why petrophysics requires engineering judgment, not just equations.

Step 4: Net Pay Determination

Not all rock with hydrocarbons is worth producing. Net pay is the thickness of reservoir rock that meets minimum quality standards. The industry-standard approach is to apply cutoffs — threshold values that a depth interval must pass to count as pay.

Common cutoffs for a clastic (sandstone) reservoir:

Parameter	Cutoff	Why
Vshale	< 0.40	Too much shale means low permeability
Effective porosity	> 0.08	Too little pore space means too little oil
Water saturation	< 0.60	More than 60% water means oil won't flow preferentially

These cutoffs are not universal. They depend on the reservoir rock type, the fluid properties, and the economics of the field. A deepwater well with high development costs might require tighter cutoffs than an onshore well.

main.py

# Net pay cutoffs
VSHALE_CUTOFF = 0.40
PHIE_CUTOFF = 0.08
SW_CUTOFF = 0.60

# Apply cutoffs
logs["FLAG_VSHALE"] = logs["VSHALE"] < VSHALE_CUTOFF
logs["FLAG_PHIE"] = logs["PHIE"] > PHIE_CUTOFF
logs["FLAG_SW"] = logs["SW"] < SW_CUTOFF

# A depth is net pay only if ALL three criteria are met
logs["NET_PAY"] = (
    logs["FLAG_VSHALE"] &
    logs["FLAG_PHIE"] &
    logs["FLAG_SW"]
)

# Calculate net pay thickness
step = 0.5  # ft per sample
net_pay_ft = logs["NET_PAY"].sum() * step
gross_interval = logs["DEPTH"].max() - logs["DEPTH"].min()

print("NET PAY SUMMARY")
print(f"  Gross interval:  {gross_interval:.0f} ft "
      f"({logs['DEPTH'].min():.0f} - {logs['DEPTH'].max():.0f} ft)")
print(f"  Net pay:         {net_pay_ft:.1f} ft")
print(f"  Net-to-gross:    {net_pay_ft/gross_interval:.2%}")
print()

# Average properties in net pay
pay = logs[logs["NET_PAY"]]
print(f"  Average properties in net pay:")
print(f"    Porosity (PHIE):  {pay['PHIE'].mean():.3f}")
print(f"    Water sat (Sw):   {pay['SW'].mean():.3f}")
print(f"    HC sat (1-Sw):    {pay['SHC'].mean():.3f}")
print(f"    Vshale:           {pay['VSHALE'].mean():.3f}")
print()

# Identify pay zones as continuous intervals
pay_mask = logs["NET_PAY"].values
zone_starts = []
zone_ends = []
in_zone = False
for i in range(len(pay_mask)):
    if pay_mask[i] and not in_zone:
        zone_starts.append(logs["DEPTH"].iloc[i])
        in_zone = True
    elif not pay_mask[i] and in_zone:
        zone_ends.append(logs["DEPTH"].iloc[i - 1])
        in_zone = False
if in_zone:
    zone_ends.append(logs["DEPTH"].iloc[-1])

print("Pay zones identified:")
for i, (top, base) in enumerate(zip(zone_starts, zone_ends)):
    thickness = base - top
    print(f"  Zone {i+1}: {top:.1f} - {base:.1f} ft ({thickness:.1f} ft)")

The Complete Petrophysical Pipeline

Bringing all four steps together into a single, reusable function:

main.py

def petrophysical_evaluation(df, gr_col="GR", rt_col="RT",
                              rhob_col="RHOB", nphi_col="NPHI",
                              depth_col="DEPTH",
                              rho_matrix=2.65, rho_fluid=1.0,
                              a=0.81, m=2.0, n=2.0, Rw=0.04,
                              vsh_cutoff=0.40, phie_cutoff=0.08,
                              sw_cutoff=0.60):
    """
    Complete petrophysical evaluation from raw log data.

    Parameters
    ----------
    df : DataFrame with log curves
    Archie parameters: a, m, n, Rw
    Matrix/fluid densities: rho_matrix, rho_fluid
    Cutoffs: vsh_cutoff, phie_cutoff, sw_cutoff

    Returns
    -------
    DataFrame with all computed properties and net pay flags.
    Dictionary with summary statistics.
    """
    result = df.copy()

    # GR endpoints
    gr_clean = result[gr_col].quantile(0.05)
    gr_shale = result[gr_col].quantile(0.95)

    # Step 1: Vshale
    result["VSHALE"] = ((result[gr_col] - gr_clean) /
                        (gr_shale - gr_clean)).clip(0, 1)

    # Step 2: Porosity
    result["PHID"] = ((rho_matrix - result[rhob_col]) /
                      (rho_matrix - rho_fluid)).clip(0, 0.45)
    result["PHIA"] = (result["PHID"] + result[nphi_col]) / 2
    result["PHIE"] = (result["PHIA"] * (1 - result["VSHALE"])).clip(0, 0.40)

    # Step 3: Water saturation (Archie)
    result["SW"] = np.nan
    valid = result["PHIE"] > 0.01
    result.loc[valid, "SW"] = (
        (a * Rw) / (result.loc[valid, "PHIE"]**m * result.loc[valid, rt_col])
    ) ** (1/n)
    result["SW"] = result["SW"].clip(0, 1)
    result["SHC"] = 1 - result["SW"]

    # Step 4: Net pay
    result["NET_PAY"] = (
        (result["VSHALE"] < vsh_cutoff) &
        (result["PHIE"] > phie_cutoff) &
        (result["SW"] < sw_cutoff)
    )

    # Summary
    step = abs(result[depth_col].diff().median())
    net_pay_ft = result["NET_PAY"].sum() * step
    pay = result[result["NET_PAY"]]

    summary = {
        "gross_interval_ft": result[depth_col].max() - result[depth_col].min(),
        "net_pay_ft": net_pay_ft,
        "net_to_gross": net_pay_ft / (result[depth_col].max() - result[depth_col].min()),
        "avg_phie": pay["PHIE"].mean() if len(pay) > 0 else 0,
        "avg_sw": pay["SW"].mean() if len(pay) > 0 else 0,
        "avg_vshale": pay["VSHALE"].mean() if len(pay) > 0 else 0,
        "gr_clean": gr_clean,
        "gr_shale": gr_shale,
    }

    return result, summary


# Run the complete evaluation
evaluated, summary = petrophysical_evaluation(logs)

print("PETROPHYSICAL EVALUATION — COMPLETE RESULTS")
print("=" * 50)
for key, val in summary.items():
    if isinstance(val, float):
        print(f"  {key:25s}: {val:.4f}")
    else:
        print(f"  {key:25s}: {val}")

This function encapsulates the entire workflow. For a new well, you load the LAS file (using the pipeline from Chapter 6), pass the DataFrame to petrophysical_evaluation(), and receive a complete analysis in seconds. Parameters can be adjusted per well or per field without rewriting any logic.

Summary

This chapter covered the core workflow of petrophysical analysis:

Gamma ray logs measure natural radioactivity and distinguish sand from shale. Shale volume (Vshale) is calculated from the GR response relative to clean sand and pure shale baselines.
Density and neutron logs provide independent porosity estimates. Effective porosity corrects for shale content, giving the porosity available for hydrocarbon storage.
Archie's equation converts resistivity and porosity into water saturation — the fraction of pore space filled with water. The remainder is hydrocarbons. The equation's parameters (a, m, n, Rw) must be calibrated to the specific reservoir.
Net pay cutoffs (Vshale, porosity, water saturation) determine which intervals are commercially producible. These cutoffs are engineering judgments, not universal constants.
A complete pipeline automates the entire workflow, making it repeatable and auditable.

The next chapter uses reservoir fluid properties — bubble point pressure, formation volume factor, viscosity, and gas solubility — to understand how the hydrocarbons we just identified will behave when we produce them.

Exercises

fitness_center

Exercise 7.1Practice

GR Endpoint Sensitivity

Using the well log dataset from this chapter, calculate Vshale using three different GR_clean / GR_shale pairs: 5th/95th percentiles (as we used)Minim...

arrow_forward