Exerciseschevron_rightChapter 7chevron_right7.8
fitness_center

Exercise 7.8

LAS-Driven Petrophysics on OD-007

Level 3
Chapter 7: Well Log Analysis
descriptionProblem

In the field you almost never get a tidy numpy array; you get a LAS file. This one is OD-007 on OML 58, four curves (GR RT RHOB NPHI) logged on a 5-ft step. One NPHI sample in the upper shale is dead and was written as the LAS null -999.25. The section is classic Niger-Delta stacking, top to bottom: an upper shale, a wet sand, an oil-pay sand beneath it, and a lower shale. Stacked deltaic bodies, each its own reservoir. Your job is to pick the one pay sand out of that stack.

Write evaluate_las(las_text) that runs the full base-case pipeline on the LAS and returns a summary dict:

  1. Read the LAS with lasio.read(io.StringIO(las_text)) and build las.df()

(the -999.25 null becomes NaN automatically; never hard-replace it).

  1. Vshale (linear): gr_clean / gr_shale are the 5th / 95th percentiles

of GR; IGR = (GR - gr_clean)/(gr_shale - gr_clean); VSHALE = clip(IGR, 0, 1).

  1. Porosity: PHID = (2.65 - RHOB)/(2.65 - 1.0) clipped 0..0.45;

PHIA = (PHID + NPHI)/2; PHIE = (PHIA * (1 - VSHALE)) clipped 0..0.40.

  1. Archie Sw (base a=0.81, m=2, n=2, Rw=0.04):

SW = ((a*Rw)/(PHIEm * RT))(1/n) clipped 0..1, only where PHIE > 0.01.

  1. Net pay (base cutoffs): VSHALE < 0.40 AND PHIE > 0.08 AND SW < 0.60.

net_pay_ft = (passing samples) * step_ft, where step_ft is the LAS depth step (here 5.0 ft).

Return {well, n_curves, net_pay_ft, phie_pay, sw_pay} where well is the LAS well name, n_curves is how many log curves came in (excluding depth), net_pay_ft is the pay footage, and phie_pay / sw_pay are the mean PHIE and mean SW over the net-pay samples only.

The discipline: the dead NPHI sample must stay NaN so it can't sneak into a pay count, and your endpoints come from the data, not a guess.

lightbulbHints (0/3)

Stuck? Reveal hints one at a time — they progress from nudge to near-solution.

codeYour solution
main.py
visibilityReveal reference solutionexpand_more

Try solving it yourself first — the hints walk you through it. The solution below is one valid approach; yours may differ and still be correct.

import lasio
import io
import numpy as np
import pandas as pd

LAS_TEXT = """~VERSION INFORMATION
 VERS.   2.0 : CWLS LOG ASCII STANDARD - VERSION 2.0
 WRAP.   NO  : ONE LINE PER DEPTH STEP
~WELL INFORMATION
 WELL.  OD-007 : Well Name
 FLD.   OML 58 : Field Name
 NULL.  -999.2500 : NULL VALUE
 STEP.FT 5.0 : Depth Step
~CURVE INFORMATION
 DEPT.FT   : Depth
 GR  .GAPI : Gamma Ray
 RT  .OHMM : Deep Resistivity
 RHOB.G/CC : Bulk Density
 NPHI.V/V  : Neutron Porosity
~A  DEPT       GR       RT        RHOB     NPHI
  8000.000   106.50     3.760    2.542     0.325
  8005.000   108.00     3.880    2.547     0.330
  8010.000   109.50     4.000    2.552     0.335
  8015.000   106.50     4.120    2.557  -999.250
  8020.000   108.00     4.240    2.542     0.330
  8025.000   109.50     3.760    2.547     0.335
  8030.000   106.50     3.880    2.552     0.325
  8035.000   108.00     4.000    2.557     0.330
  8040.000   109.50     4.120    2.542     0.335
  8045.000   106.50     4.240    2.547     0.325
  8050.000    38.00     2.068    2.312     0.220
  8055.000    39.50     2.134    2.317     0.225
  8060.000    36.50     2.200    2.303     0.215
  8065.000    38.00     2.266    2.308     0.220
  8070.000    39.50     2.332    2.312     0.225
  8075.000    36.50     2.068    2.317     0.215
  8080.000    38.00     2.134    2.303     0.220
  8085.000    39.50     2.200    2.308     0.225
  8090.000    36.50     2.266    2.312     0.215
  8095.000    30.00    40.280    2.257     0.160
  8100.000    31.50    35.720    2.243     0.165
  8105.000    28.50    36.860    2.248     0.155
  8110.000    30.00    38.000    2.252     0.160
  8115.000    31.50    39.140    2.257     0.165
  8120.000    28.50    40.280    2.243     0.155
  8125.000    30.00    35.720    2.248     0.160
  8130.000    31.50    36.860    2.252     0.165
  8135.000    28.50    38.000    2.257     0.155
  8140.000    30.00    39.140    2.243     0.160
  8145.000    31.50    40.280    2.248     0.165
  8150.000    28.50    35.720    2.252     0.155
  8155.000   112.00     3.880    2.567     0.340
  8160.000   113.50     4.000    2.553     0.345
  8165.000   110.50     4.120    2.558     0.335
  8170.000   112.00     4.240    2.562     0.340
  8175.000   113.50     3.760    2.567     0.345
  8180.000   110.50     3.880    2.553     0.335
  8185.000   112.00     4.000    2.558     0.340
  8190.000   113.50     4.120    2.562     0.345
  8195.000   110.50     4.240    2.567     0.335
"""

STEP_FT = 5.0
A, M, N, RW = 0.81, 2.0, 2.0, 0.04
RHO_MA, RHO_FL = 2.65, 1.0


def evaluate_las(las_text):
    """Run the base-case petrophysical pipeline on a LAS string.

    Return {well, n_curves, net_pay_ft, phie_pay, sw_pay}.
    """
    las = lasio.read(io.StringIO(las_text))
    df = las.df()  # depth index; -999.25 already NaN. NPHI sample stays NaN.

    well = las.well["WELL"].value
    n_curves = df.shape[1]  # GR, RT, RHOB, NPHI -> 4

    # Linear Vshale from data-driven endpoints.
    gr_clean = df["GR"].quantile(0.05)
    gr_shale = df["GR"].quantile(0.95)
    igr = (df["GR"] - gr_clean) / (gr_shale - gr_clean)
    vshale_frac = igr.clip(0.0, 1.0)

    # Porosity: density porosity averaged with neutron, shale-corrected.
    phid = ((RHO_MA - df["RHOB"]) / (RHO_MA - RHO_FL)).clip(0.0, 0.45)
    phia = (phid + df["NPHI"]) / 2.0
    phie_frac = (phia * (1.0 - vshale_frac)).clip(0.0, 0.40)

    # Archie water saturation (base parameters), only where porosity is real.
    sw_frac = (((A * RW) / (phie_frac ** M * df["RT"])) ** (1.0 / N)).clip(0.0, 1.0)
    sw_frac[phie_frac <= 0.01] = 1.0

    # Net pay (base cutoffs). NaN anywhere in a row makes the comparison False,
    # so the dead NPHI sample can never be counted as pay.
    pay = (vshale_frac < 0.40) & (phie_frac > 0.08) & (sw_frac < 0.60)
    net_pay_ft = int(pay.sum()) * STEP_FT

    return {
        "well": well,
        "n_curves": n_curves,
        "net_pay_ft": net_pay_ft,
        "phie_pay": float(phie_frac[pay].mean()),
        "sw_pay": float(sw_frac[pay].mean()),
    }


result = evaluate_las(LAS_TEXT)
print(result)
print(f"OD-007 net pay: {result['net_pay_ft']} ft at Sw={result['sw_pay']:.3f}")

lockCopying code is a Full Access feature.