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

$\boxed{\;\frac{\partial^{2} u}{\partial t^{2}} \;=\; c^{2}(z)\,\frac{\partial^{2} u}{\partial z^{2}}\;}$

Here $u(z, t)$ is the wavefield (pressure, or displacement, or any scalar quantity the wave is carrying) at depth $z$ and time $t$ , and $c(z)$ 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, at the crest of a pulse, it accelerates hardest; where it is straight, at the inflection points on its flanks, it does not accelerate at all. The boxed form assumes constant density. With density varying it becomes $\\frac{1}{\\rho c^{2}}\\,\\partial\_{tt} p = \\partial\_z\\!\\left(\\frac{1}{\\rho}\\,\\partial\_z p\\right)$ for the pressure $p$ , which is what the figure below solves.

2. Plane-wave solutions

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

u(z, t) = f(z - ct) + g(z + ct)

for any twice-differentiable $f$ and $g$ . These are travelling waves: $f$ moving down at speed $c$ , $g$ moving up. 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 the rock changes, it splits: some continues forward (transmission), some bounces back (reflection). For normal incidence on a step in acoustic impedance $Z = \\rho c$ , the reflection coefficient is

R = \\frac{Z\_2 - Z\_1}{Z\_2 + Z\_1}

and the pressure carried on is $T = 1 + R$ . It is the impedance contrast, not the velocity contrast, that sets $R$ . This one number is the entire point of reflection seismology. In the figure below a pulse leaves a source near the top of a two-layer earth: drag the time to follow it, then change the velocities and the density below the interface and watch what comes back.

At the starting setting the lower layer is faster and denser, $Z\_2 = 7.20$ MRayl against $Z\_1 = 4.40$ , so the impedances predict $R = +0.241$ . The loop returns +0.239, measured on the trace in (e), 381 ms after the direct wave, and 1.239 of the pulse goes on down. The space-time diagram (d) stacks every snapshot of the wavefield side by side, depth down and time across. In a uniform layer the pulse draws straight lines of slope $1/V$ , one going down and one going up; at the interface a new line peels off and runs back toward the surface. The trace recorded at the receiver, plate (e), is one row of that diagram: it is your first shot record. Set $V\_2$ to 2500 m/s and $\\rho\_2$ to 1.76 g/cc and the impedances match: the line in (d) still bends to the new slope, yet nothing comes back. Stretch the picture to two dimensions, record along the surface, and each straight reflection line becomes the hyperbola you will read 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:

\\partial\_{tt} u \\approx \\frac{u^{n+1} - 2u^{n} + u^{n-1}}{(\\Delta t)^{2}},\\qquad \\partial\_{zz} u \\approx \\frac{u\_{i+1} - 2u\_{i} + u\_{i-1}}{(\\Delta z)^{2}}

Plugging in and solving for the future time step:

u\_i^{n+1} = 2u\_i^{n} - u\_i^{n-1} + \\left(\\frac{c\_i\\,\\Delta t}{\\Delta z}\\right)^{2} (u\_{i+1}^{n} - 2u\_i^{n} + u\_{i-1}^{n})

That single line, run in a loop, is what the figure does. Its version carries the density as well, and with constant density it is exactly this line. 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:

\\text{CFL} = \\frac{c\_{\\max}\\,\\Delta t}{\\Delta z} \\;\\le\\; 1 \\qquad (\\le 1/\\sqrt{d}\\ \\text{in}\\ d\\ \\text{dimensions})

Intuition: information in the wave equation propagates at speed $c$ , so a time step of $\\Delta t$ carries the wave a distance $c\\,\\Delta t$ . If that distance exceeds one cell $\\Delta z$ , 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 figure sets $\\Delta t = \\mathrm{CFL}\\,\\Delta z/V\_{\\max}$ from its Courant number, 0.50 at the start, which gives steps of 0.417 ms on its 2.5 m grid. Push the Courant number just past 1 and the run explodes within milliseconds of the pulse entering the faster layer: at 1.05 it blows up at 275 ms. Stability is not accuracy, though. The grid also needs enough cells in the shortest wavelength, about 8 or more for this scheme; with 10 m cells and a 30 Hz pulse it has 2.7, and the pulse smears into a trailing ripple.

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, such as the free surface at the top of marine data, where the pressure is zero and $R = -1$ .
Absorbing / sponge layer: taper the amplitude near the boundary, or impose a one-way wave equation there, so the energy leaves the grid. The figure absorbs with Mur’s one-way condition, which is exact in 1D at CFL 1; switch its bottom edge to Reflecting and the edge sends the wave straight back like a false interface.
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 ( $N^{3}$ in 3D) times the number of time steps, and because CFL ties $\\Delta t$ to the cell size, halving the cell size costs about $2^{4} = 16$ times more:

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 one-line update the figure runs. The physics is simple; the engineering is heroic.

8. Why this closes Part 0

Look at everything Part 0 has put in your hands:

Section 0.2 complex numbers, the language of frequency-domain modeling.
Section 0.3 convolution, how source wavelet + reflectivity make the trace.
Section 0.4 Fourier, frequency-domain wave propagation; convolution theorem.
Section 0.5 sampling, the grid and Nyquist of the numerical simulator.
Section 0.6 Z-transform, the language of discrete wave-equation stencils.
Section 0.7 linear algebra, the matrix form of every finite-difference step.
Section 0.8 random variables, the noise in the observed traces the simulator is matching.
Section 0.9 optimization, the outer loop that updates $c(z)$ to match data.
Section 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 imaging and inversion algorithm past Section 0.10, from migration to FWI, 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 Section 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.