Seismic Wave Simulation Workflow

SeisGo’s simulation module provides a 1-D finite-difference acoustic wave solver and a layered velocity model builder, useful for synthetic testing and method validation.

Building a layered velocity model

build_vmodel creates a fine-grid, linearly increasing velocity model with optional per-layer perturbations:

from seisgo import simulation
import numpy as np

z, v, rho = simulation.build_vmodel(
    zmax=500,          # maximum depth / model length (m or km)
    dz=1.0,            # grid spacing
    nlayer=5,          # number of velocity layers
    vmin=1500,         # minimum velocity (units consistent with dz)
    vmax=3500,         # maximum velocity
    rhomin=1800,       # minimum density
    rhomax=2500,       # maximum density
    zmin=0,            # starting depth (default 0)
    layer_dv=None,     # optional per-layer fractional perturbation
)

Adding velocity anomalies

layer_dv is an array of fractional velocity perturbations applied multiplicatively to each layer:

# Make layer 3 (0-indexed) 10% slower
layer_dv = np.zeros(5)
layer_dv[2] = -0.10

z, v, rho = simulation.build_vmodel(
    zmax=500, dz=1.0, nlayer=5,
    vmin=1500, vmax=3500,
    rhomin=1800, rhomax=2500,
    layer_dv=layer_dv,
)

Running a 1-D finite-difference simulation

fd1d_dx4dt4 solves the 1-D acoustic wave equation with:

4th-order spatial accuracy O(Δx⁴)
4th-order temporal accuracy O(Δt⁴) via the Adams–Bashforth scheme

from seisgo import simulation
import numpy as np

# Model setup
x = np.arange(0, 500, 1.0)    # spatial grid (m)
vmodel = v                      # from build_vmodel above
rho_model = rho

tout, seis = simulation.fd1d_dx4dt4(
    x=x,
    dt=1e-4,           # time step (s) — must satisfy CFL condition
    tmax=0.3,          # total simulation time (s)
    vmodel=vmodel,
    rho=rho_model,
    xsrc=50,           # source grid index
    xrcv=250,          # receiver grid index
    stf_freq=30,       # source central frequency (Hz) for Ricker wavelet
    stf_shift=None,    # time shift (default: 3/stf_freq)
    stf_type="ricker", # "ricker" or "gaussian"
    t_interval=1,      # output time sub-sampling factor
)

CFL stability condition

The solver automatically checks that dt satisfies the Courant–Friedrichs–Lewy condition:

$$dt \leq 0.5 \cdot \frac{\Delta x_{\max}}{v_{\max}}$$

A ValueError is raised if the provided dt is too large.

Source time functions

`stf_type`	`stf_freq` meaning
`"ricker"`	Central frequency (Hz)
`"gaussian"`	Width parameter σ (s)

Both wavelets are available from seisgo.utils:

from seisgo import utils

t_ricker, w_ricker   = utils.ricker(dt=1e-4, f=30, t0=0.05)
t_gauss,  w_gauss    = utils.gaussian(dt=1e-4, width=0.02, shift=0.05)

Plotting the synthetic seismogram

import matplotlib.pyplot as plt

plt.figure(figsize=(9, 3))
plt.plot(tout, seis, "k", lw=0.8)
plt.xlabel("Time (s)")
plt.ylabel("Pressure amplitude")
plt.title("Synthetic seismogram — 1-D FD simulation")
plt.tight_layout()
plt.show()

Full example: synthetic CorrData for method testing

from seisgo import simulation, utils
import numpy as np

# Build model
x = np.arange(0, 1000, 2.0)
z, v, rho = simulation.build_vmodel(1000, 2.0, 6, 1500, 3200, 1800, 2600)

# Simulate with source at one end, receiver at mid-point
tout, seis = simulation.fd1d_dx4dt4(
    x, dt=2e-4, tmax=0.5, vmodel=v, rho=rho,
    xsrc=10, xrcv=250, stf_freq=20, stf_type="ricker",
)

print(f"Simulated {len(tout)} time samples, dt={tout[1]-tout[0]:.2e} s")