Absorbing Boundary Conditions: Definition, FD Modeling, and PML

Absorbing boundary conditions (ABCs) are numerical algorithms applied along the outer edges of a finite computational domain in seismic wave simulation to suppress artificial reflections that would otherwise propagate back into the model interior and contaminate the wavefield solution. Because finite-difference (FD) and finite-element (FE) grids are necessarily bounded in space, wave energy that reaches a grid boundary has nowhere to go; without an appropriate termination strategy, the wavefield reflects off that boundary exactly as it would off a physical impedance contrast. ABCs mimic the behavior of an infinite, reflection-free medium so that outgoing waves simply exit the domain without returning. They are indispensable in acoustic and elastic modeling, full-waveform inversion (FWI), reverse-time migration (RTM), and any other workflow in which a synthetic wavefield must replicate propagation through a geologically realistic, effectively unbounded Earth.

Key Takeaways

• ABCs prevent artificial boundary reflections in finite computational domains, preserving solution accuracy in seismic modeling and inversion workflows.
• The Perfectly Matched Layer (PML), introduced by Jean-Pierre Berenger in 1994, is the current gold-standard ABC: it absorbs waves of any frequency and angle of incidence with near-zero reflection, outperforming all earlier paraxial and sponge-layer approaches.
• PML thickness typically ranges from 10 to 30 grid points; thinner layers save memory but risk numerical instability, while thicker layers improve absorption at the cost of added computation.
• Stability of PML formulations in anisotropic (VTI/TTI) elastic media is an active research area; standard PML can become unstable in strongly anisotropic models, driving adoption of Convolutional PML (CPML) and Multi-Axial PML (M-APML) variants.
• In open-source codes such as Devito and SeisFlows, PML is configured via a small set of parameters (damping profile, layer thickness, reflection coefficient target), making implementation accessible without deriving the boundary equations from scratch.

How Absorbing Boundary Conditions Work

In a standard FD seismic simulation, the computational domain is a rectangular (2D) or cuboid (3D) grid of nodes at which particle velocity or pressure is updated at each time step according to the elastic or acoustic wave equation. When a propagating wavefront reaches the edge of this grid, the boundary conditions imposed there determine what happens to the wave energy. A free-surface boundary (zero normal stress) reflects waves back with opposite polarity, which is physically correct at the Earth-air interface but disastrous on the other five faces of the model where no physical reflector exists. A rigid (Dirichlet) boundary reflects waves with the same polarity. Both are numerically simple but geophysically wrong for the sides and bottom of a synthetic model.

The earliest practical ABCs were paraxial approximations, developed by Reynolds (1978) and Clayton and Engquist (1977, 1980). These one-way wave equation operators allow waves traveling roughly perpendicular to the boundary to pass through with minimal reflection, but their absorption efficiency degrades rapidly for waves arriving at oblique incidence angles above roughly 30 degrees. Because a realistic seismic model contains energy at many angles, paraxial ABCs leave residual reflections that smear into the image. Sponge layers (also called absorbing boundary zones or damping zones) address this by attaching a thick buffer strip around the model interior within which an exponential or quadratic amplitude damping function multiplies the wavefield at each time step. Sponge layers work adequately but require a large number of grid points (often 40-80) to achieve strong attenuation, adding significant memory and compute cost.

The Perfectly Matched Layer solved the oblique-incidence problem analytically. Berenger showed that by splitting each field component and introducing a complex frequency-shifted stretching variable in the boundary layer, the governing wave equations could be transformed so that any outgoing wave, regardless of frequency or angle, decays exponentially within the layer without any reflection at the interior-PML interface. The theoretical reflection coefficient is exactly zero for a continuous PML, and in practice (after discretization) it is typically below -80 dB with a layer only 10-20 grid points thick. The computational overhead of PML over a sponge layer of equivalent effectiveness is modest: PML needs far fewer grid points to achieve the same absorption, so total memory consumption is usually lower. For industrial-scale 3D FWI and RTM, where model domains can exceed 500x500x200 grid points at fine sampling, this difference is material. The Convolutional PML (CPML) variant, formulated by Komatitsch and Martin (2007), is particularly robust because it handles evanescent waves and low-frequency components better than the original split-field formulation, and it integrates cleanly into second-order velocity-stress FD stencils.

PML Formulation and Design Parameters

Within the PML layer, spatial derivatives in the direction normal to the boundary are replaced by a complex stretched coordinate:

d/dx → (1 / s_x) d/dx,   where s_x = 1 + d(x) / (iω + α(x))

Here, d(x) is a positive, spatially varying damping coefficient (typically a polynomial profile peaking at the outer edge of the PML), ω is angular frequency, and α(x) is an optional attenuation shift that improves absorption of evanescent and very low-frequency waves. The target reflection coefficient R determines how steeply d(x) must rise:

d_max = -(n+1) V_P ln(R) / (2 L_PML)

where n is the polynomial order (typically 2 or 3), V_P is P-wave velocity at that boundary, and L_PML is the physical thickness of the PML in meters. A common design target is R = 10^-3 (reflection coefficient -60 dB). With V_P = 3,000 m/s, n = 2, and L_PML = 300 m (10 grid points at 30 m spacing), d_max works out to approximately 207 s^-1. In the time domain, the CPML auxiliary memory variable evolves as a first-order ODE alongside the main wavefield variables, requiring storage of one additional array per spatial direction per wavefield component inside the PML layer. For a 3D elastic simulation, this adds roughly 30-40% memory within the PML zone; because the zone is thin relative to the full model, total memory increase is small.

PML stability in anisotropic media is the main practical challenge. In tilted transverse isotropic (TTI) media, certain combinations of Thomsen parameters (δ, ε, θ) produce group-velocity vectors that are not parallel to the phase-velocity vector, allowing wave energy to re-enter the interior from the PML. The Multi-Axial PML (M-APML) and Anisotropic PML (APML) formulations address this by applying damping in multiple coordinate directions simultaneously, though at added implementation complexity. In practice, most production FWI codes either use CPML with empirical stability checks or revert to a thick sponge layer at model boundaries where TTI instability is anticipated.

Applications in Seismic Workflows

ABCs appear in every forward-modeling step of modern quantitative seismic workflows. In full-waveform inversion, the acoustic or elastic wave equation is solved forward in time for each source, then the data residual is back-propagated as an adjoint wavefield; the gradient of the objective function is the zero-lag cross-correlation of forward and adjoint fields. Any boundary reflection in either wavefield masquerades as a physical reflector in the gradient, degrading convergence. FWI practitioners therefore use CPML with R targets of 10^-4 or better, and typically extend the PML to all six faces of the 3D model, including the free surface when marine data are modeled without sea-surface multiples.

In reverse-time migration, the source wavefield is forward-propagated with ABCs, saved to disk (or recomputed via checkpointing), and the receiver wavefield is back-propagated with the same ABCs operating in time-reversed mode. Imperfect ABCs create "rabbit ears" artifacts at the image edges: false high-amplitude zones that degrade interpretation, particularly on deep targets near the model boundary. A 20-point CPML border reduces these artifacts to below the noise floor on typical model sizes. Industry RTM implementations in packages such as Schlumberger's Omega, TGS's proprietary FD engine, and open-source frameworks like Devito all implement CPML by default.

Acoustic versus elastic formulations have different ABC requirements. Acoustic modeling treats the subsurface as a fluid (only P-waves), and PML requires only one auxiliary variable per direction. Elastic modeling includes both P- and S-wave modes, and since P- and S-wave velocities differ, the PML stretching parameters must accommodate both simultaneously; a single damping profile that absorbs P-waves well may under-absorb S-waves if V_S is much slower. CPML handles this by applying the same formal stretching to all stress and velocity components; d_max is commonly set based on P-wave velocity (the faster mode), accepting slightly more residual reflection on S-waves as a practical compromise.

It is important to note that the term "absorbing boundary conditions" also appears in reservoir simulation, where it means something entirely different: a boundary condition on a flow model that allows fluid to leave the domain without reflecting (i.e., a no-flow or constant-pressure outer boundary). The two usages are technically analogous but operate in completely different physical contexts and should not be confused.