lasio / welly          bruges                    segyio / obspy
[Well Log Prep]     --> [Rock Physics]         --> [Synthetics / Tie]
  |                      |                          |
  Load LAS/DLIS          Elastic moduli             Wavelet extraction
  QC & despike           Gassmann fluid sub         Reflectivity series
  Resample curves        AVO intercept/gradient     Convolve synthetic
  Extract Vp, Vs, rho    Backus averaging           Well-seismic tie

import lasio
import numpy as np
from welly import Well

# Load well with sonic, density, and shear sonic
las = lasio.read('well.las')
df = las.df().dropna()

# Extract elastic logs
depth = df.index.values
vp = 1e6 / df['DT'].values       # P-wave velocity (m/s) from sonic (us/ft)
vs = 1e6 / df['DTS'].values      # S-wave velocity (m/s) from shear sonic
rho = df['RHOB'].values * 1000   # Density (kg/m3) from g/cc

# QC: check ranges
assert np.all(vp > 1500) and np.all(vp < 7000), "Vp out of range"
assert np.all(vs > 500) and np.all(vs < 4000), "Vs out of range"
assert np.all(rho > 1500) and np.all(rho < 3200), "Density out of range"

# If no shear sonic, estimate from Vp
# Castagna mudrock line: Vs = 0.8621 * Vp - 1172.4 (m/s)
if 'DTS' not in df.columns:
    vs = 0.8621 * vp - 1172.4
    vs = np.maximum(vs, 300)  # Floor for shallow sediments

import bruges

# Elastic moduli from velocities
K = bruges.rockphysics.moduli.bulk(vp=vp, vs=vs, rho=rho)   # Bulk modulus
G = bruges.rockphysics.moduli.shear(vs=vs, rho=rho)          # Shear modulus
E = bruges.rockphysics.moduli.youngs(K=K, G=G)               # Young's modulus
nu = bruges.rockphysics.moduli.poissons(vp=vp, vs=vs)        # Poisson's ratio
AI = bruges.rockphysics.moduli.impedance(vp=vp, rho=rho)     # Acoustic impedance
SI = bruges.rockphysics.moduli.impedance(vp=vs, rho=rho)     # Shear impedance

# Vp/Vs ratio (key AVO indicator)
vp_vs = vp / vs

# Gassmann fluid substitution
# Replace brine with gas in reservoir interval
phi = df['NPHI'].values  # Porosity

# Mineral and fluid properties
K_mineral = 36.6e9   # Quartz bulk modulus (Pa)
K_brine = 2.6e9      # Brine bulk modulus (Pa)
rho_brine = 1050     # Brine density (kg/m3)
K_gas = 0.02e9       # Gas bulk modulus (Pa)
rho_gas = 100        # Gas density (kg/m3)

# Dry rock modulus from saturated (reverse Gassmann)
K_sat = K.copy()
K_dry = bruges.rockphysics.fluidsub.vrh(
    volumes=[1-phi, phi],
    moduli=[K_mineral, K_brine]
)[0]

# Forward Gassmann: substitute gas for brine
K_sat_gas = bruges.rockphysics.fluidsub.gassmann(
    k_sat=K_sat, k_fl=K_brine, k_min=K_mineral,
    phi=phi, k_fl2=K_gas
)

# Updated density with gas
rho_gas_sat = rho - phi * rho_brine + phi * rho_gas

# Updated velocities
vp_gas = np.sqrt((K_sat_gas + 4/3 * G) / rho_gas_sat)
vs_gas = np.sqrt(G / rho_gas_sat)

# AVO intercept and gradient (Shuey approximation)
# For a single interface between layers i and i+1
for i in range(len(vp) - 1):
    rc = bruges.reflection.shuey(
        vp1=vp[i], vs1=vs[i], rho1=rho[i],
        vp2=vp[i+1], vs2=vs[i+1], rho2=rho[i+1],
        theta=np.arange(0, 40, 1)
    )

# AVO classification from intercept (R0) and gradient (G)
# Class I:   R0 > 0, G < 0  (hard sand, dim with offset)
# Class II:  R0 ~ 0, G < 0  (near-zero, polarity reversal)
# Class III: R0 < 0, G < 0  (soft sand, bright with offset)
# Class IV:  R0 < 0, G > 0  (very soft, dim with offset)

# Zoeppritz exact for full offset range
rc_exact = bruges.reflection.zoeppritz(
    vp1=vp[i], vs1=vs[i], rho1=rho[i],
    vp2=vp[i+1], vs2=vs[i+1], rho2=rho[i+1],
    theta=np.arange(0, 50, 1)
)

# Create reflectivity series
rc_series = bruges.reflection.reflectivity(vp, rho)

# Create wavelet
duration = 0.128  # seconds
dt = 0.002        # sample rate (2ms)
wavelet = bruges.filters.ricker(duration=duration, dt=dt, f=25)

# Convolve to create synthetic
synthetic = np.convolve(rc_series, wavelet, mode='same')

# If tying to seismic, extract wavelet from seismic trace
import segyio
with segyio.open('seismic.sgy') as f:
    near_trace = f.trace[0]  # Nearest trace to well
    # Extract statistical wavelet from trace
    # Or use bruges.filters for analytic wavelets

- [ ] Load well logs (Vp, Vs, Rho, porosity) with lasio
- [ ] QC logs: check ranges, despike, fill gaps
- [ ] If no Vs log: estimate from Castagna or Greenberg-Castagna
- [ ] Calculate elastic moduli and impedances with bruges
- [ ] Run Gassmann fluid substitution (brine to gas/oil)
- [ ] Compare Vp, Vs, density, impedance before/after fluid sub
- [ ] Compute AVO response at target interface (Shuey or Zoeppritz)
- [ ] Classify AVO response (Class I-IV)
- [ ] Generate synthetic seismograms for both fluid scenarios
- [ ] Plot AVO crossplot (intercept vs gradient)

- [ ] Load well logs and seismic trace at well location
- [ ] Create time-depth relationship from check shots or sonic
- [ ] Convert logs to time domain
- [ ] Extract wavelet from seismic (statistical or deterministic)
- [ ] Generate synthetic seismogram from reflectivity * wavelet
- [ ] Cross-correlate synthetic with seismic trace
- [ ] Adjust stretch/squeeze to optimize tie
- [ ] Report correlation coefficient

- [ ] Load thin-bed well logs at fine sampling (0.5 ft)
- [ ] Define averaging window (e.g., quarter wavelength at target frequency)
- [ ] Apply Backus averaging to get effective anisotropic elastic properties
- [ ] Compare fine-scale vs upscaled reflectivity
- [ ] Assess thin-bed tuning effects