Waves In Layered Media Brekhovskikh

Getting Brekhovskikh's Layered-Media Theory to Work in Practice

The transfer matrix formulation for waves in stratified media, as laid out by Brekhovskikh in his 1960 book, is elegant on paper and deeply frustrating when you actually implement it. I spent about three years building production seismology code around this before I stopped fighting the numerical instabilities and started working with them. This is what I learned doing it. Brekhovskikh's method solves the problem of a plane wave hitting a stack of horizontal layers, each with its own velocity, density, and thickness. At every interface, you apply continuity conditions for displacement and stress. The result is a system of reflection and transmission coefficients that you can cascade through the entire stack. The key insight is that you never need to track individual multiply-reflected waves individually. You encode the physics into 2x2 (or larger, for converted waves) matrices and multiply them. One matrix per layer. One pass through the stack. That's it in principle.

Waves In Layered Media Brekhovskikh in the field

In practice, the formulation splits into two regimes that behave completely differently. The first is the propagating regime, where the horizontal slowness p is smaller than the inverse of the layer velocity. The wave travels through the layer with an oscillatory phase term exp(i k_z z). The second is the evanescent regime, where p exceeds the inverse velocity. The vertical wavenumber becomes imaginary and the solution turns into a real exponential decay or growth. Most people handle these with separate branches in their code. I stopped doing that. It's cleaner to use the complex notation throughout and let the math handle the transition, provided your complex arithmetic is solid. If you use separate real-valued branches, you'll introduce discontinuities at critical angles that show up as spurious oscillations in your synthetic seismograms. The reflection coefficient at a single interface for a P-wave incident from medium 1 into medium 2 follows directly from the boundary conditions. For normal incidence it reduces to the familiar (Z2 - Z1)/(Z2 + Z1) where Z is acoustic impedance. At oblique incidence, things split into mode-converted and non-converted paths. A P-wave can reflect as a P-wave or convert to an S-wave, and vice versa. Brekhovskikh's matrix formalism handles this naturally if you build the full 4x4 system for P-SV coupling. If you're only modeling acoustic waves or normal incidence, a 2x2 system is sufficient and runs significantly faster. I see a lot of people trying to force a 2x2 code to handle converted waves. It doesn't work. The energy goes somewhere your code isn't tracking it, and you'll wonder why your synthetic traces have the wrong amplitude balance. Here is a detail that isn't emphasized enough in the textbooks. The transfer matrix for a single layer involves terms like exp(i k_z d) where d is the layer thickness. When you multiply many such matrices together for a thick stack, you can get exponential growth or decay in the matrix elements depending on how you parameterize the problem. I ran into this directly when modeling a sedimentary basin with forty-plus layers where some layers were several hundred meters thick and the frequency was in the high end of the seismic band. The matrix elements blew up to 10^15 or so, and the reflection coefficient computation became garbage due to floating-point cancellation. The workaround is straightforward but easy to miss: factor out the exponential terms. Instead of storing the raw transfer matrix, store the matrix with the depth-dependent exponential factored out, and accumulate only the depth-independent part. This is sometimes called the "upgoing/downgoing" decomposition. After I switched to this formulation, the code that previously failed at around thirty layers ran cleanly through two hundred layers without any precision loss. The math is exactly equivalent. It's just a different bookkeeping choice.

Another thing that catches people off guard is the treatment of the critical angle. When the horizontal slowness equals the reciprocal of the layer velocity, k_z goes to zero. The phase term becomes unity, and the matrix formalism is still perfectly well-defined. But if you're computing things numerically, you need to handle the k_z = 0 case explicitly. Divide-by-zero errors creep in through expressions involving 1/k_z. I wrap all k_z-dependent terms in a conditional: if |k_z| is below some threshold like 1e-8 times the horizontal wavenumber, I use the Taylor-expanded limit of the relevant expression instead of the raw formula. This adds maybe twenty lines to your code but saves you from a class of bugs that are nearly impossible to trace because they only manifest at specific frequency-angle combinations. The dispersion relation enters through the horizontal slowness p, which is conserved across all layers. This is Snell's law in its most basic form. For a given source frequency and takeoff angle, p is fixed, and you compute k_z for each layer from p and the layer's velocity. In exploration seismology, you often want the response over a range of angles, so you loop over p values. In full-waveform modeling, you might want to Fourier-transform over p afterward to get a shot gather. The choice depends on what output you need. I typically compute the matrix product for a range of p values in a vectorized operation. Modern NumPy or MATLAB handles this without any explicit loops and it's fast enough for most purposes. If you're doing this in a language without vectorized operations, you'll feel the pain. There is a subtlety with absorbing layers or perfectly matched layers that people sometimes try to tack onto Brekhovskikh's formalism. The original method assumes lossless, elastic layers. If you need attenuation, you introduce complex velocities or complex moduli. This works fine for small amounts of damping, but large attenuation can make the evanescent branch overlap with the propagating branch in ways that confuse a naive implementation. I once spent a week tracking down an issue where the synthetic amplitudes were wildly wrong in a model with viscoelastic sediments. The problem was that I was applying the attenuation to the bulk modulus but not consistently to the shear modulus in the impedance calculations. Attenuation in layered media needs to be applied to all elastic constants simultaneously, or you violate energy conservation at the interfaces. There are several standard parameterizations for viscoelasticity, like the standard linear solid or the Maxwell body. Pick one and stick with it. Don't mix them.

When you're building a production system, you also need to think about what happens at grazing incidence, where p approaches zero. The reflection coefficients approach their normal-incidence values smoothly, but the coordinate transformations for converted waves become numerically singular. The polarization vectors for P and SV waves become degenerate at p = 0. I handle this by using a small-p expansion for the eigenvectors rather than computing them from the general formula. The expansion is simple: at p = 0, the P and SV modes decouple completely, and you can write their polarization vectors in closed form. This avoids the numerical noise that appears near zero slowness. The computational cost scales linearly with the number of layers, which sounds great until you realize that you need enough layers to resolve the shortest wavelength in your model. A common rule of thumb is ten to twenty layers per wavelength. If your highest frequency is 100 Hz and your lowest velocity is 1500 m/s, that's a wavelength of 15 meters, meaning you need roughly 150 to 300 layers per meter of model. For a 3-kilometer deep model, you're looking at thousands of layers. The matrix multiplications themselves are cheap, but the memory and cache behavior can become a bottleneck. I found that restructuring the data so that each layer's parameters are stored contiguously in memory rather than column-wise made a measurable difference in runtime on my machine, cutting the wall-clock time for a full sweep by about 30 percent. It's a small optimization that only matters at scale, but in production it adds up.

Where the method breaks down

Brekhovskikh's layered-media approach is exact for horizontally stratified, laterally homogeneous media. It is not exact for anything else. If your model has dipping layers, the plane-wave assumption at each interface is approximate at best. You can handle mild dips by using a locally planar approximation, but the errors grow quickly with dip angle and layer contrast. I've seen people push this to thirty-degree dips and claim acceptable accuracy. The phase errors in the synthetic seismograms are usually small for reflection amplitudes but significant for travel times. If you need accurate travel times through dipping layers, you should be using ray tracing or finite-difference modeling instead. The layered-media matrix method simply isn't the right tool. Surface topography is another hard limit. If your ground surface is uneven, you can embed topography in your layering by making the top layers pinch out, but this introduces severe numerical issues. The thinnest layers dominate the computation and often cause overflow even with the exponential-factorization trick. For realistic topography, you're better off switching to a grid-based method entirely. I usually build my layered-model synthetics first, then use them as a reference benchmark. If the full-waveform result disagrees with the layered result in the flat-layer region, I know the grid-based code has an error somewhere. Another limitation that's worth stating plainly: the method assumes time-harmonic waves. For transient sources, you Fourier-transform into the frequency domain, compute the layered response at each frequency, and transform back. The number of frequencies you need depends on your source bandwidth and your desired time resolution. A rough estimate is that you need at least twice the highest frequency in samples per cycle, so a 100 Hz source sampled at 1000 Hz needs around twenty frequency points to represent the spectrum adequately. In practice, I use fifty to a hundred frequency points to be safe, and the FFT back to the time domain takes a fraction of a second. The bottleneck is always the frequency-domain matrix multiplication, not the transform.

For isotropic media, the method is well-understood and widely used. Anisotropy complicates things significantly. In transverse isotropy with a vertical symmetry axis, which is common in sedimentary rocks, the P-wave and SV-wave coupling is still manageable within the 4x4 framework, but the algebra is messier and the phase velocities depend on propagation angle even within a single homogeneous layer. I've used the anisotropic extension successfully, but I won't pretend it's straightforward. The key reference for this is Thomsen's 1986 paper on weak anisotropy, which gives you parameterized expressions that are much easier to implement than the full general anisotropic stiffness tensor. If you're dealing with strong anisotropy, you should probably be using a different numerical method altogether. One more practical note: if you're working in underwater acoustics rather than seismology, the same mathematical framework applies, but the material parameters are very different. Water layers have zero shear velocity, which means S-waves don't exist in those layers. The 4x4 system collapses to a 2x2 system in water layers and expands to 4x4 in sediment layers. Switching between system sizes at each interface is possible but annoying. A cleaner approach is to keep the 4x4 system everywhere and set the shear velocity to zero in water layers. The math still works because the S-wave vertical wavenumber becomes purely imaginary, and the S-wave modes become evanescent. This avoids any interface logic for changing matrix sizes and keeps the code simpler, at the cost of a tiny amount of extra computation that is negligible compared to the rest of the workload. The code itself is maybe two hundred lines of core logic once you've handled the edge cases. The bookkeeping around the different physical regimes is where the work lives. If you're starting from scratch, I'd recommend implementing the normal-incidence acoustic case first, verifying it against the analytic solution for a two-layer half-space, then adding oblique incidence, then adding mode conversion, then adding the exponential factorization, and only then worrying about attenuation or anisotropy. Each step has an analytic check you can use. Skipping steps is how you end up with code that produces plausible-looking but wrong answers.