Sandbox wave equation V1D
Understanding the 1D damped acoustic wave equation through an interactive sandbox.
Most seismic workflows treat wave physics as a black box. We run migrations, inversions, and synthetics, but the governing equation itself, the engine of all propagation, is usually hidden behind layers of abstraction.
This sandbox for the wave equation reverses that habit. It strips the equation down to its simplest useful form, one spatial dimension, and lets you watch every term at work. The goal is not realism for its own sake, but clarity, the clarity we need in these times of black boxes and automatic processes. When you can see each physical contribution act independently, wave propagation stops being magic and becomes an intelligible, structured process.
Code and reproducibility
The full Python implementation accompanying this article is openly available (for free). The code is released as a stable, versioned reference and can be cloned, executed, and modified freely.
GitHub repository:
https://github.com/maxseg2021/wave-sandbox
Stable release (Version 1D):
Sandbox wave equation V1D – Release v1.0
This article discusses exactly the physics and numerical implementation contained in the v1.0 release. Future extensions (e.g. 2D propagation) will be published as separate, versioned releases.
The governing equation
The sandbox is built around the one dimensional damped acoustic wave equation
∂²u(x,t)/∂t² + 2γ ∂u(x,t)/∂t = c(x)² ∂²u(x,t)/∂x² + s(t) δ(x − xₛ)
where u(x,t) is the wavefield, x is position, and t is time. The velocity c(x) may vary spatially to represent layering. The source wavelet s(t) injects energy at position xₛ, represented mathematically by the Dirac delta (δ).
Each term has a clear physical meaning.
The second time derivative represents inertia. It is the tendency of a disturbed medium to keep moving.
The spatial curvature term describes elastic restoring forces. Differences between neighbouring points pull the wave forward or backward. The squared velocity c(x)² controls how fast this interaction propagates.
The first order time derivative introduces damping, which removes energy from the system at a rate proportional to particle velocity. Without this term, the wave would oscillate indefinitely.
The source term is the only place where energy enters the system. Numerically, the delta function is implemented as a localized forcing at the source grid point.
Attenuation and the role of Q
Intrinsic attenuation is controlled through the quality or attenuation factor Q. Instead of implementing a full frequency dependent attenuation model, the sandbox uses the standard relationship
γ = π f₀ / Q
where f₀ is the dominant frequency of the source.
This captures an essential physical fact: high frequency waves attenuate faster than low frequency waves in the same medium. A high Q implies weak attenuation and slow decay. A low Q implies strong energy loss. Because γ depends on frequency, changing the source spectrum immediately changes the rate at which energy is dissipated.
From equation to code
The equation is discretized using second order finite differences in space and time. For grid point i and time step n, the update rule is
uᶰ⁺¹ᵢ = (2 − 2γΔt) uᶰᵢ − (1 − 2γΔt) uᶰ⁻¹ᵢ
+ cᵢ² (Δt²/Δx²) (uᶰᵢ₊₁ − 2uᶰᵢ + uᶰᵢ₋₁) + sᶰ δᵢ,ᵢₛ
The first two terms carry inertia and damping from previous time steps. The spatial derivative term propagates curvature through the medium. The last term injects the source wavelet at the source index.
Numerical stability is enforced using the Courant condition
Δt ≤ 0.9 Δx / max(c(x))
ensuring that information does not propagate faster than the grid allows.
Source scaling and energy
Source strength is parameterized in terms of source energy, expressed in the sandbox as an equivalent source mass. Amplitude (A) scales with the square root of energy (E):
A ∝ √E
reflecting the physical relationship E ∝ A². This avoids unrealistic linear scaling and ensures that doubling energy does not double displacement.
Different source types (Ricker, Gaussian, vibroseis) inject different frequency content, but they all enter the equation through the same physical mechanism (available in the phyton code).
Geometrical spreading
Because the wave equation itself is strictly one dimensional, geometrical spreading is not built into the propagation. Instead, a post propagation correction is applied at the receiver:
1D spreading: constant amplitude
2D spreading: amplitude ∝ 1/√r
3D spreading: amplitude ∝ 1/r
This separation is deliberate. It prevents artificial amplitude growth inside faster layers and keeps propagation physics independent from geometric energy dilution. Spreading affects detectability and signal to noise, not travel times or wave shapes.
Noise and detectability
Noise is added only at the receiver
d(t) = u(xᵣ,t) + n(t)
with n(t) drawn from a Gaussian distribution with a specified RMS level.
This distinction is crucial. In the equation, waves never disappear. Experimentally, they disappear when their amplitude drops below the noise floor. By tying detection thresholds automatically to the noise level, the sandbox reproduces the intuitive behaviour that weak sources die earlier without altering the underlying physics.
Layering and reflections
Velocity layering is introduced through c(x). Reflections and transmissions emerge naturally from the spatial derivative term. No reflection coefficients are imposed. The wave equation enforces continuity and boundary conditions on its own.
This is one of the most instructive aspects of this sandbox. A single interface is enough to demonstrate reflection, transmission, and interference without any additional assumptions.
Why start in one dimension
This sandbox is intentionally minimal. With one spatial dimension there is no moveout, no hyperbola, and no lateral ambiguity. In that sense it is analogous to a zero offset sonic or a vertical seismic experiment. Source, receiver, time, and rock.
This simplicity is the point. Before extending the model to two dimensions, where geometrical spreading becomes intrinsic and diffractions appear naturally, the one dimensional engine must be understood thoroughly.
The purpose of this sandbox is mastery through simplification. When every wiggle in a trace can be traced back to a specific term in a simple equation, wave physics stops being a black box and becomes a tool you can reason about.
Version two will extend this same logic to two spatial dimensions. But the foundation is here. Just as an example a Pyodide extraction from the code (run the model then play):
Run the 1D wavefield in your browser for two layers (Pyodide)
Comments
Post a Comment