Exploration Seismology: Full Waveform Inversion and Subsurface Acoustic Velocity Models

Exploration Seismology: Full Waveform Inversion and Subsurface Acoustic Velocity Models

Mapping subterranean geological formations for geothermal reservoir discovery, aquifer mapping, and sub-sea carbon sequestration ($CO_2$ storage monitoring) relies on Exploration Seismology. Compressed airguns tow behind survey vessels, firing acoustic impulse bursts down into the seabed; hydrophone sensor arrays record reflected pressure wavelets over several seconds.

1. 2D Acoustic Wave Equation Formulation

In isotropic acoustic media with heterogeneous velocity distribution $v(\mathbf{x})$, the propagation of pressure fields $p(\mathbf{x}, t)$ generated by source $s(\mathbf{x}, t)$ satisfies:

$$\frac{1}{v^2(\mathbf{x})} \frac{\partial^2 p}{\partial t^2} - \nabla^2 p = s(\mathbf{x}, t)$$

At geological boundaries where acoustic impedance $Z = \rho \cdot v$ shifts (e.g., transition from shale to sandstone), reflection coefficient $R$ dictates wave energy partitioning:

$$R = \frac{\rho_2 v_2 - \rho_1 v_1}{\rho_2 v_2 + \rho_1 v_1}$$

       Marine Reflection Seismology Acquisition Setup
       Vessel (Airgun Array)
       |---- Source Shot
       v 
       ~~~~~~~~~~~~~~~~~~~~ Sea Surface ~~~~~~~~~~~~~~~~~~~~~
       ~~~~~~~~~~~~~~~~~~~~ Sea Floor ~~~~~~~~~~~~~~~~~~~~~~~
              |             /                           | Downward   / Reflected     \ Ocean Bottom Nodes
              V Wavefront /  Acoustic Ray   \ (OBN Seismic Receivers)
       ------------------------------------------------------ Strata 1 (v1)
       \\\\\\\\\\\\\\\\\\\\\\\\\\\
       ------------------------------------------------------ Strata 2 (v2 > v1)

2. Full Waveform Inversion (FWI) via Adjoint State Method

Traditional seismic imaging merely tracks arrival travel times (tomography). Full Waveform Inversion (FWI) solves a non-linear inverse problem minimizing the complete $L_2$ norm misfit between observed field waveforms $d_{\text{obs}}$ and numerically simulated wavefields $d_{\text{syn}}(v)$:

$$\chi(v) = \frac{1}{2} \sum_{\text{shots}} \int_0^T | p_{\text{syn}}(\mathbf{x}r, t; v) - d{\text{obs}}(\mathbf{x}_r, t) |^2 dt$$

Calculating gradients with respect to millions of subterranean velocity grid points without astronomical computational cost leverages the Adjoint State Method:

1. Forward propagate the synthetic wavefield forward in time: $p(\mathbf{x}, t)$.
2. Back-propagate the data residuals (adjoint wavefield) reverse in time: $\lambda(\mathbf{x}, t)$.
3. Compute the velocity gradient via zero-lag cross-correlation:

$$\frac{\partial \chi}{\partial v(\mathbf{x})} = -\frac{2}{v^3(\mathbf{x})} \int_0^T \lambda(\mathbf{x}, t) \frac{\partial^2 p(\mathbf{x}, t)}{\partial t^2} dt$$

import numpy as np

def compute_seismic_synthetic_trace(velocity_profile, density_profile, ricker_freq=25.0, dt=0.001, duration_s=1.0):
    '''
    1D convolutional seismic synthetic trace generation.
    '''
    impedance = velocity_profile * density_profile
    # Reflection coefficient series
    rc = (impedance[1:] - impedance[:-1]) / (impedance[1:] + impedance[:-1])

    # Ricker wavelet generator
    time = np.arange(-0.06, 0.06, dt)
    tau = np.pi * ricker_freq * time
    ricker = (1.0 - 2.0 * (tau**2)) * np.exp(-(tau**2))

    # Convolution
    synthetic_trace = np.convolve(rc, ricker, mode='same')
    return synthetic_trace

Full Waveform Inversion processes petabytes of field telemetry to render sub-kilometer crustal velocity volumes, providing geological certainty for clean geothermal energy extraction.