Wave equation in 30 minutes

Processing Prerequisites

Learning objectives

  • State the 1D acoustic wave equation and identify what each term means physically
  • Predict where reflections, transmissions, and wavefronts appear in a layered earth
  • Explain the finite-difference scheme and the CFL stability condition
  • Bridge from the PDE to every imaging algorithm that follows: migration, FWI, modeling

Everything in seismic starts with the wave equation. The source injects energy; the equation tells the wavefield how to evolve; the receiver records whatever arrives. Migration, FWI, and all physics-based modeling are just different ways to solve (or invert) this equation.

1. The 1D acoustic wave equation

  ∂2u∂t2  =  c2(x) ∂2u∂x2  \boxed{\;\frac{\partial^{2} u}{\partial t^{2}} \;=\; c^{2}(x)\,\frac{\partial^{2} u}{\partial x^{2}}\;}

Here u(x, t) is the wavefield (pressure, or displacement, or any scalar quantity the wave is carrying), and c(x) is the local propagation velocity. The equation says: the acceleration of the wavefield at any point equals velocity squared times its local curvature in space. Where the wave is sharply curved (an edge of a pulse), it accelerates; where it is flat, it does not.

2. Plane-wave solutions

In a homogeneous medium (c constant), the general solution is

u(x,t)=f(x−ct)+g(x+ct)u(x, t) = f(x - ct) + g(x + ct)

for any twice-differentiable f and g. These are travelling waves: one moving right at speed c, one moving left. The shape is set by the source; the speed is set by the medium.

3. Reflection, transmission, and interfaces

When a wave hits an interface where c changes, energy splits: some continues forward (transmission), some bounces back (reflection). For normal incidence on a step in impedance Z = ρc, the reflection coefficient is

R=Z2−Z1Z2+Z1R = \frac{Z_2 - Z_1}{Z_2 + Z_1}

This one number is the entire point of reflection seismology. Drag the source around in the widget below and switch between velocity profiles to see reflections emerge.

The 1D wave equation∂²u/∂t² = c² ∂²u/∂x²t=0t=Δtt=2Δtpropagation: cInteractive figure — enable JavaScript to adjust c, source position, and BCs.

The bottom panel is your first shot record. Each row is the wavefield at one instant in time; time grows downward. In a homogeneous medium, the pulse draws two straight lines — one to the left, one to the right, each at slope 1/c. Add an interface and you see a new branch peel off: the reflected wave starts at the interface and draws a second line back toward the surface. This is exactly the picture you will be interpreting on every shot gather in Parts 1 and 2.

4. Finite-difference discretization

Solving the PDE on a computer means discretizing. Replace derivatives with differences:

∂ttu≈un+1−2un+un−1(Δt)2,∂xxu≈ui+1−2ui+ui−1(Δx)2\partial_{tt} u \approx \frac{u^{n+1} - 2u^{n} + u^{n-1}}{(\Delta t)^{2}},\qquad \partial_{xx} u \approx \frac{u_{i+1} - 2u_{i} + u_{i-1}}{(\Delta x)^{2}}

Plugging in and solving for the future time step:

uin+1=2uin−uin−1+(ci ΔtΔx)2(ui+1n−2uin+ui−1n)u_i^{n+1} = 2u_i^{n} - u_i^{n-1} + \left(\frac{c_i\,\Delta t}{\Delta x}\right)^{2} (u_{i+1}^{n} - 2u_i^{n} + u_{i-1}^{n})

That single line, run in a loop, is what the widget does. Every 2D and 3D modeler is the same idea in more indices; the engineering is keeping the arithmetic intensity up and the memory access patterns cache-friendly.

5. The CFL condition

For this explicit scheme to be stable, the Courant-Friedrichs-Lewy condition must hold:

CFL=cmax⁡ ΔtΔx  <  1\text{CFL} = \frac{c_{\max}\,\Delta t}{\Delta x} \;\lt\; 1

Intuition: information in the wave equation propagates at speed c, so a time step of Δt carries the wave a distance c·Δt. If that distance exceeds one cell Δx, the numerical update cannot keep up with the physics, and amplitudes explode. Every time you see “CFL violation” in a processing log, this is why.

The widget quietly picks Δt to keep CFL near 0.45 for safety. The info strip shows the resulting CFL.

6. Boundary conditions

Real earth is infinite; computer grids are finite. If you do nothing, waves hit the edge of the grid and reflect back in, contaminating the solution. Three fixes:

  • Reflecting (Dirichlet or Neumann): the wave bounces off. Useful only when the boundary is a real reflector (e.g. free surface at the top of marine data).
  • Absorbing / sponge layer: taper the amplitude near the boundary to bleed energy away. The widget uses a small sponge (you can see very faint damping near the edges).
  • PML (Perfectly Matched Layer): the gold standard in production modeling — a specially designed absorbing layer that damps waves with minimal reflection across frequencies.

7. From 1D to real seismic

The real world is 3D with three degrees of freedom per point (elastic) and variable density. The PDE generalizes; the computational cost scales with the number of grid points cubed and the number of time steps:

  • 2D acoustic: the foundation of 2D RTM. Cheap enough for classroom use.
  • 3D acoustic: what most 3D RTM production codes solve.
  • 3D elastic: three coupled PDEs for displacement vector components; needed for AVO-faithful modeling.
  • 3D anisotropic elastic: adds TTI or orthorhombic parameters; the state of the art in advanced FWI.

Every one of those is a generalization of the four lines of code in the widget’s update rule. The physics is simple; the engineering is heroic.

8. Why this closes Part 0

Look at everything Part 0 has put in your hands:

  • §0.2 complex numbers — the language of frequency-domain modeling.
  • §0.3 convolution — how source wavelet + reflectivity make the trace.
  • §0.4 Fourier — frequency-domain wave propagation; convolution theorem.
  • §0.5 sampling — the grid and Nyquist of the numerical simulator.
  • §0.6 Z-transform — the language of discrete wave-equation stencils.
  • §0.7 linear algebra — the matrix form of every finite-difference step.
  • §0.8 random variables — the noise in the observed traces the simulator is matching.
  • §0.9 optimization — the outer loop that updates c(x) to match data.
  • §0.10 wave equation — the physics the simulator itself solves.

Those nine topics, threaded together, are FWI. They are also migration. They are also inversion. You now have every piece.

**The one sentence to remember**

The wave equation says the wavefield’s acceleration equals velocity-squared times its spatial curvature; finite differences turn that into a loop; CFL tells you the largest stable time step; and every processing algorithm past §0.10 is an application of this one PDE.

What comes next

Part 1 switches from mathematics to the field. We walk through land acquisition, marine acquisition, SEG-Y headers, noise characterization, and survey design sanity — the data the processing algorithms later will operate on. Part 0 is complete. See you in §1.1.

References

  • Claerbout, J. F. (1985). Imaging the Earth’s Interior. Blackwell.
  • Yilmaz, Ö. (2001). Seismic Data Analysis (2 vols.). SEG.
  • Tarantola, A. (1984). Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49, 1259.
  • Virieux, J., Operto, S. (2009). An overview of full-waveform inversion in exploration geophysics. Geophysics, 74, WCC1.

This page is prerendered for SEO and accessibility. The interactive widgets above hydrate on JavaScript load.