Mathematical Problems

in Geophysical Wave Propagation

George Papanicolaou

Abstract. We review several aspects of the mathematical theory of wave propagation in random media with particular emphasis on topics of geophysical interest.

Contents

1 Introduction 404 1.1 Why geophysical wave propagation . . . . . . . . . . . . . . . . . . 404 1.2 Random media or environments in general . . . . . . . . . . . . . . 404 1.3 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . 406

2 General notions about waves in random media 406 2.1 Scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 406 2.2 Types of waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 407 2.3 Coherent and incoherent fields . . . . . . . . . . . . . . . . . . . . 407 2.4 Localization and transport . . . . . . . . . . . . . . . . . . . . . . . 407 2.5 Nonlinearity and randomness . . . . . . . . . . . . . . . . . . . . . 408 2.6 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . 408 2.7 Parameter estimation and imaging . . . . . . . . . . . . . . . . . . 409

3 The transport regime 409 3.1 Waves to transport . . . . . . . . . . . . . . . . . . . . . . . . . . . 411 3.2 Transport for electromagnetic and elastic waves . . . . . . . . . . . 414 3.3 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 415 3.4 The diffusive regime . . . . . . . . . . . . . . . . . . . . . . . . . . 415 3.5 Parameter identification and inverse problems . . . . . . . . . . . . 415

4 The localization regime 416 4.1 Pulse reflection from randomly layered media . . . . . . . . . . . . 416 4.2 Scale separation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 417 4.3 Localization regime asymptotics . . . . . . . . . . . . . . . . . . . . 418 4.4 Time domain localization . . . . . . . . . . . . . . . . . . . . . . . 420 4.5 Statistical inverse problems . . . . . . . . . . . . . . . . . . . . . . 421 4.6 Reflection and transmission of time harmonic plane elastic waves . 423 4.7 Pulse stabilization and imaging . . . . . . . . . . . . . . . . . . . . 423 1 Introduction

1.1 Why geophysical wave propagation

In deciding what to present as a plenary lecture in applied mathematics and probability at ICM-98, I considered several areas with which I am familiar and decided to focus on geophysical wave propagation for a couple of reasons. One, a technical one, is that inhomogeneities are strong and highly anisotropic so that the modeling and analysis of wave propagation in the earth’s crust is mathematically interesting and quite difficult. The other reason is closely related to what my view of modern applied mathematics is: the creation and development of a mathematical environment for physical, economic, biological or other phenomena. This involves active participation of the mathematician in the quantitative modeling, in the analysis, in the computations, as well as in the interpretation of results and assessment of the effectiveness of the modeling.

The resulting mathematical methodology will be uneven, from routine offthe-shelf toolbox applications to entirely uncharted problems that need new ideas and techniques, and it is up to the mathematician to decide what the right mix of mathematical sophistication and rough heuristics should be. An overly mathematical approach will impede communication with nonmathematical specialists who value results and do not care much for mathematical generality. Accepting the conventional wisdom in a field, and concentrating on technical mathematical issues, is not a good idea either. Geophysical wave propagation is a case in point. It is fair to say that wave localization is virtually unknown to geophysicists. But, as I will try to explain in this lecture, wave localization is quite important in exploration geophysics because, among other things, it influences the resolution of seismic imaging and the effective depth penetration of seismic probes. What is the best way to approach these problems mathematically?

A few years ago, K. Aki, a distinguished seismologist whose ideas about the role of crustal inhomogeneities in seismic wave propagation have been very influential, heard a seminar that I gave on wave localization and asked this question: How can one tell from seismic observations that wave localization has taken place? Electronic wave localization in semiconductors goes back forty years [1], with the strong participation of mathematicians during the last twenty years, so we should be able to say quite a bit, as I will try to explain in this lecture. But Aki’s question is a profound one that leads to the most complex and least understood issue in geophysical wave propagation, the localization-transport transition. It is a pragmatic, operational question which reminds us that great intellectual challenges can have humble, unpretentious origins. I think that it takes a mathematician to answer Aki’s question and perhaps it will be one that does it.

1.2 Random media or environments in general

Mathematical Problems in Geophysical Wave Propagation 405

I will treat the earth’s crust as a random medium, that is, as an elastic medium with density and Lame parameters that are random functions of space. The equations of linear elastic wave propagation become now stochastic partial differential equations. Initial and boundary conditions must also be specified and they could bring in additional randomness, from modeling the rough surface of the earth. At this level of generality the randomness is nothing more than variable coefficients and non-flat boundaries, so general linear PDE methods can deal with everything (symmetric hyperbolic systems). If dissipation is important, and it is in some contexts, it can be put into the equations in different ways. There is no general agreement on how to best model dissipation analytically and this is an interesting issue that I will not address here.

But even this much is somewhat grudgingly accepted by geophysicists. I regularly hear comments like: there is only one earth and it is not changing all that fast, so where is the statistical ensemble of realizations coming from? If stochastic modeling is to be criticized along such lines then why are we modeling the Dow Jones Industrial Average, or some other index or asset price, as a stochastic process? There is only one realization of the DJIA just as there is only one realization of the earth. What the stochastic processes model is uncertainty, lack of information and its consequences when only imperfect and sparse observations or measurements are available, and even desirable. The notion of ‘effective’ medium is very much part of the mathematical physics of the 19th century, of Maxwell, Rayleigh and others, which is why equations with constant coefficients have any relevance at all in modeling. The conceptual barrier seems to come up when one thinks of fluctuations.

It is not an accident, therefore, that in one of the first instances of wave propagation in random media, natural light propagation through a turbulent atmosphere, astrophysicists at the turn of the 20th century did not go to Maxwell’s equations (or the wave equation if the vector nature of light waves can be neglected) but developed a new, phenomenological theory, the radiative transport theory, to interpret observed phenomena. There are a few isolated attempts to consider random media, with fluctuations, during the first half of the 20th century but it is with the advent of radar and sonar during in the forties that random waves emerge as a subject. Keller’s papers in the sixties [2] where very influential because they were the first ones written by a mathematician, who thought about the conceptual foundations and separated heuristics from legitimate calculations. It was also in the sixties that the connection between radiative transport theory and stochastic wave equations was clarified, as I will discuss in section 3.1.

Atmospheric wave propagation, from radio to radar to optical frequencies, and underwater sound propagation, from 20 hertz to kilohertz, were the main applications driving the theory of wave propagation in random media in the seventies and are discussed in Ishimaru’s book [3]. It is interesting to note that the notion of wave localization is nowhere to be found in this book. Random media in seismology appeared first in the mid eighties in a simple version of radiative transport [4]. Transport theory is now just beginning to become mainstream in seismology as is seen from the recent book of Fehler and Sato [5]. But wave localization is not discussed in this book either. A treatment of waves in random media that deals extensively with wave localization is given by Ping Sheng [6].

What is wave localization anyway? I will explain it in some detail in section 4.3 but, roughly, it is when random inhomogeneities trap wave energy in a finite region and do not allow it to spread as it would normally. Random media behave then like periodic media that have band-gap spectra, allowing wave propagation in some frequency ranges but not in others. It is remarkable that this happens for random media that are not close to periodic ones at all. Mathematically, it is shown that wave or wave-like operators with stationary (translation invariant) random coefficients in unbounded regions have discrete spectra [7]. Discrete spectra means that the wave energy in each mode initially will remain there for ever, oscillating in time but not propagating out to infinity. In three dimensional wave propagation this can happen only when parameter fluctuations are very large. This is not the case for electromagnetic waves in the atmosphere or sound waves in the ocean. The fluctuations are weak, a few percent, and when they are important they lead to radiative transport, which allows spreading of wave energy in diffusive rather than wave-like manner.

Where then do we see wave localization in classical wave propagation? We see it when wave energy is channeled, by a waveguide, by a transmission line, by an optical fiber, by strong anisotropy due to layering in the lithosphere, etc. We also see it in nearly periodic structures. Waves in an one dimensional random medium will localize, even if the fluctuations in the medium parameters are weak. In geophysical wave propagation and elsewhere (in optical localization) a key issue is the identification of structures, more complicated than simple channeling or periodicity, that tend to enhance the onset of wave localization by random fluctuations. This is the localization-transport transition problem.

1.3 Acknowledgements

I would like to thank several friends and colleagues with whom I worked in the past: R. Burridge, J.B. Keller, W. Kohler and B.S. White with whom I have been associated for a long time, and my more recent collaborators: M. Asch, G. Bal, J. Berryman, J. Bronski, O. Dorn, J.P. Fouque, F. Herrmann, R. Knapp, P. Lewicki, A. Nachbin, M. Postel, L. Ryzhik, Y. Samuelides, P. Sheng, K. Solna and S. Weinryb.

2 General notions about waves in random media

2.1 Scales

There are three basic length scales in wave propagation phenomena:

• The typical wavelength λ

• The typical propagation distance L

• The typical size of the inhomogeneities l

In geophysical wave propagation it is difficult to associate a ‘typical’ scale that is characteristic of the inhomogeneities. The density and local speed of propagation of waves vary on many scales. We may think of l as a typical correlation length. When the standard deviation of the fluctuations is small then the most effective interaction of the waves with the random medium will occur when l ∼ λ, that is, the wavelength is comparable to the correlation length. And this interaction will not be observable unless the propagation distance is large (L >> λ). If propagation distances are short, a few wavelengths or correlation lengths, then effective medium theories will work fine. There will be a deterministic propagation speed (for scalar waves), the effective speed, with which energy will propagate as if the medium were deterministic. The effective medium theory will be valid also when the wavelength is long compared to propagation distances, even if the correlation length is short (l << L) and the fluctuations have a large standard deviation. This is the homogenization limit.

Of course this rough way of thinking with scales does not capture the effect of a waveguide geometry, or the effective dimension of the propagation phenomenon. But thinking with scales is very useful and, with some experience, it can become a very good heuristic tool.

2.2 Types of waves

It is classical waves, solutions of the wave equation or more general symmetric hyperbolic systems, that we want consider, rather than electronic waves which are solutions of the Schrödinger equation. The waves are vector fields in general, as with electromagnetic waves which are solutions of Maxwell’s equations or elastic waves where the elastic displacement field is a solution of the elastic wave equations. Mode conversion, the transfer of energy from compression to shear waves for example, is an important effect in random media. So is polarization, which is associated with vector waves all of whose components travel with the same speed. Polarization tends to get lost in a random medium and the way this happens is an important way to make inferences about the nature of the propagation environment.

2.3 Coherent and incoherent fields

When the random fluctuations of the medium parameters are small then the random fluctuations in the solutions will be small, if the propagation distances are not too big. The mean solution, the coherent field, will carry most of the energy. As the waves propagate their fluctuating component, the incoherent field, gets more energy. The total energy is conserved, if there is no dissipation, but the coherent field loses energy and slows down. This behavior of the coherent field is something that can be calculated easily and is well established in the engineering literature.

2.4 Localization and transport

If fluctuations are weak and propagation distances large, most of the wave energy will be incoherent. In seismology, for example, after the first arrival from a disturbance far away the seismogram is dominated by strong fluctuations from multiple scattering. The later part, the coda of a seismogram is mostly incoherent field measurement. It is in this regime that radiative transport is a good approximation. It allows accurate calculation of the envelopes of the seismograms without resolving the detailed fluctuations. A new scale enters the description of propagation phenomena: the mean free path. This is a length scale that gives an indication of the importance of multiple scattering and is much more relevant than the correlation length of the inhomogeneities.

Wave localization is total trapping of the wave energy by scattering from the random inhomogeneities. It is the regime where fluctuation phenomena dominate so we have little intuition for what should happen. For one thing the random fluctuations must be very strong and the structure of the propagating medium must be special (a channeling medium or an ordered, periodic structure). In the lithosphere fluctuations in the speeds of propagation of elastic waves can be as large as 15% and they can be highly anisotropic, with horizontal correlation lengths much larger than vertical ones. Localization manifests itself in fat codas of seismograms, or codas with envelopes that decay slowly. This is a clear indication that there is a lot of multiple scattering going on. Moreover, radiative transport would tend to underestimate the size of the codas indicating that a different analytical theory is needed. What is missing at present is a robust and effective criterion for discriminating between these two situations.

2.5 Nonlinearity and randomness

Nonlinearity and randomness interact significantly only in very special situations, as in soliton propagation in optical fibers or when high intensity laser beams interact with material inhomogeneities. Nonlinearity is rarely an issue in seismic wave propagation except very near sources. In one dimensional wave propagation both nonlinearity and randomness are strongly felt and a long-standing problem is the analysis of their interaction. Is there, for example, wave localization when we have nonlinearities? This is a very difficult question that cannot be answered by a yes or no. The phenomena depend sensitively on the exact setup of the problem: the form of the nonlinearity, the various scales associated with the inhomogeneities and the propagation phenomenon, and the form of the excitation [8, 9, 10, 11, 12, 13].

2.6 Numerical simulations

At the dawn of the 21st century, when computational power is doubling every two years or so, and computational cost is dropping to the point where a good laptop computer today is more powerful than the Cray I supercomputer of the late seventies, why is anybody interested in analytical methods? We have the computational power to simulate anything we want and we have the ability to make detailed and extensive measurements, which in seismology result in huge data sets. What could mathematical analysis contribute in this context?

Being skeptical about the utility of mathematical analysis and believing that we can compute or simulate everything we need may appear naive to a mathematician but it is increasingly the dominant view in many fields, in geophysical wave propagation for example.

The fact is that if we want to understand the behavior of seismic codas we cannot rely on direct numerical simulations. If the typical wavelength is of the order of 3-5 km and we want to calculate a synthetic seismogram 1000 km from the source we need a spatial grid that has at least five points per wavelength, and more if we want to simulate accurately random fluctuations in the parameters. In a realistic three dimensional setup it is impossible to generate numerical solutions that will yield a 3 second synthetic seismogram with millisecond resolution. What is even more important to realize is that we should not really want to do this because with radiative transport theory seismic coda envelopes can be calculated. What is holding up realistic numerical computations is not computing power but analysis: we do not have good enough transport theoretic boundary conditions on the earth’s surface and at interfaces. The mean free path may be as large as 20-30 km and Monte Carlo methods can give reasonably accurate solutions using a high-end workstation. Transport theory does what is called ‘sub-grid’ modeling in computational fluid dynamics. We do not have to resolve the small scale inhomogeneities if we can do some analysis, which is in fact difficult but doable.

2.7 Parameter estimation and imaging

Imaging of the earth’s interior is a challenge that will be with us for a very long time because the inhomogeneities are so strong. In exploration seismology, where seismic probing can generate huge data sets, the issue is not so much good algorithms for imaging but low complexity algorithms. Efficient compression of geophysical data sets is perhaps the most urgent problem that exploration seismology faces at present.

It appears at first that this has nothing to do with waves. Wavelets or other tools for compression from signal processing come to mind, and they are being used. If noise effects are ignored and if the typical wave length of a probing pulse is 100-150 m (for shorter wavelengths dissipation effects are much stronger), we cannot expect image resolution better that 25 m or so at a depth of a few kilometers. And if noise and multiple scattering are to be taken into consideration it is not at all clear what the achievable resolutions are without some compensation. Noise reducing methods (stacking) that are used in imaging are not so effective. Much more needs to be done analytically here. Imaging itself, without noise, is based on variants of a backward wave propagation method (migration) that has now a substantial theoretical basis [14, 15, 16].

The best compression method is to go from the seismic data to the image itself, of course, so good compression has to be adapted to the specific data set and its structure. But there must be interesting algorithms, yet to be found, that are somewhere between know-nothing methods like wavelet decomposition and thresholding, and know-all full imaging.

3 The transport regime

Radiative transport is a phenomenological theory that was introduced to describe the propagation of light intensity through the Earth’s atmosphere. It has been applied successfully to many other problems of wave propagation in a complex medium. In its simplest form, let a(t, x, k) denote the angularly resolved energy density defined for all wave vectors k, position x and time t. Because of interaction with the inhomogeneous medium through which it propagates, a wave with wave

vector k may be scattered into any other direction k̂′, where k̂ = \frac{k}{|k|}. The

transport equation gives the energy balance

\frac{∂a(t, x, k)}{∂t} + ∇_{k}ω(x, k) · ∇_{x}a(t, x, k) − ∇_{x}ω(x, k) · ∇_{k}a(t, x, k)  (1)
= ∫_{R^{n}} σ(x, k, k′)a(t, x, k′)dk′ − Σ(x, k)a(t, x, k).

Here n is the dimension of space (n = 2 or 3), ω(x, k) is the local frequency at position x of the wave with wave vector k, the differential scattering cross-section σ(x, k, k′) is the rate at which energy with wave vector k′ is converted to wave energy with wave vector k at position x, and

∫ σ(x, k′, k)dk′ = Σ(x, k)  (2)

is the total scattering cross-section. The function σ(x, k, k′) is nonnegative and usually symmetric in k and k′. The left side of (1) is the total time derivative of a(t, x, k) at a point moving along a trajectory in phase space (x, k) and may be written as a Liouville equation

\frac{∂a}{∂t} = \{ω, a\} ,  (3)

where \{f, g\} = ∑_{i=1}^{n} (\frac{∂f}{∂x_{i}} \frac{∂g}{∂k_{i}} − \frac{∂f}{∂k_{i}} \frac{∂g}{∂x_{i}} ) is the Poisson bracket. The right side of

(1) represents the effects of scattering.

The transport equation (1) is conservative when (2) holds because then

∫ ∫ a(t, x, k)dxdk = const

independent of time. Absorption may be accounted for easily by letting the total scattering cross-section be the sum of two terms

Σ(x, k) = Σ_{sc}(x, k) + Σ_{ab}(x, k)

where Σ_{sc}(x, k) is the total scattering cross-section given by (2) and Σ_{ab}(x, k) is the absorption rate.

The radiative transport equation (1) was derived from the microscopic equations in the sixties and seventies by many authors (see [17] for references). A nice overview of these methods and results is presented in a recent review [18]. We have recently considered scattering of high frequency waves in a random medium [17] and established validity of the radiative transport theory for scalar and vector waves, including mode conversion and polarization in the following regime:

• Distances of propagation L are much larger than the wave length λ, • The medium parameters vary on the scale comparable to the wave length,

• The mismatch between the inhomogeneities and the background medium is small,

• Absorption is small.

This regime arises in many physically important situations. In seismic wave propagation, teleseismic events can be modeled by radiative transport equations [4, 5].

3.1 Waves to transport

Transport equations for the phase space wave energy densities are constructed [17, 19, 20] as follows. We assume here that the space domain is R^{3} (n = 3) and deal with acoustic waves. The acoustic equations for the velocity v and pressure p are

ρ\frac{∂v}{∂t} + ∇p = 0  (4)

κ\frac{∂p}{∂t} + ∇ · v = 0.

This system may be written in a general form of a symmetric hyperbolic system (with convention of summation over repeated indices):

A(x)\frac{∂u}{∂t} + D^{j} \frac{∂u}{∂x^{j}} = 0,  (5)

where u = (v, p), and x ∈ R^{n}. The matrix A(x) = diag(ρ, ρ, ρ, κ) is symmetric and positive definite and the matrices D^{j} are symmetric and independent of x and t. We consider high frequency solutions of (5). Physically this means that the typical wave length λ of the initial data is much smaller than the overall

propagation distance L with ε = \frac{λ}{L} ≪ 1. The spatial energy density for the solutions of (5) is given by

E(t, x) = \frac{ρv^{2}}{2} + \frac{κp^{2}}{2} = \frac{1}{2}(A(x)u(t, x) · u(t, x)) = \frac{1}{2}A_{ij}(x)u_{i}(t, x)ū_{j}(t, x) (6)

and the flux F(x) by

F_{i}(t, x) = pv = \frac{1}{2}(D^{i}u(t, x) · u(t, x)).  (7)

We have the energy conservation law

\frac{∂E}{∂t} + ∇ · F = 0.  (8)

and thus the total energy is conserved:

\frac{d}{dt} ∫ E(t, x)dx = 0.  (9)

The high frequency limit ε → 0 of the energy density E(t, x) is described in terms of the Wigner transform, which is defined by

W_{ε}(t, x, k) = (\frac{1}{2π} )^{n} ∫ e^{ik·y}u_{ε}(t, x − εy/2)u_{ε}^{∗}(x + εy/2)dy,  (10)

where u_{ε}(t, x) is the solution of (5). The Wigner transform W_{ε} is a 4 × 4 Hermitian matrix. Its limit as ε → 0 is called the Wigner distribution and is denoted by W (t, x, k). The limit Wigner matrix is not only Hermitian but also positive definite. The limit energy density and flux are expressed in terms of W (t, x, k) by

E(t, x) = \frac{1}{2} ∫ Tr(A(x)W (t, x, k))dk

and

F_{i}(t, x, k) = \frac{1}{2} ∫ Tr(D^{i}W (t, x, k))dk.

The limit Wigner distribution may be decomposed over different wave modes in a way that generalizes the plane wave decomposition in a homogeneous medium. The dispersion matrix of the system (5) is defined by

L(x, k) = A^{−1}(x)k_{i}D^{i} = ( \begin{matrix}0 & 0 & 0 & k_{1}/ρ \\ 0 & 0 & 0 & k_{2}/ρ \\ 0 & 0 & 0 & k_{3}/ρ \\ k_{1}/κ & k_{2}/κ & k_{3}/κ & 0 \end{matrix} ) .  (11)

It has one double eigenvalue ω_{1} = ω_{2} = 0 and two simple eigenvalues

ω_{f} = v|k| , ω_{b} = −v|k| ,  (12)

where |k| = \sqrt{k_{1}^{2} + k_{2}^{2} + k_{3}^{2}} and v is the sound speed

v = \frac{1}{\sqrt{κρ}}.  (13)

The corresponding basis of eigenvectors is

b_{1} = \frac{1}{\sqrt{ρ}}(z^{(1)}(k), 0)^{t}, b_{2} = \frac{1}{\sqrt{ρ}}(z^{(2)}(k), 0)^{t},

b_{f} = (\frac{k̂}{\sqrt{2ρ}}, \frac{1}{\sqrt{2κ}})^{t}, b_{b} = (\frac{k̂}{\sqrt{2ρ}}, −\frac{1}{\sqrt{2κ}})^{t},  (14)

where the vectors k̂, z^{(1)}(k) and z^{(2)}(k) form an orthonormal triplet. The eigenvectors b_{1}(k) and b_{2}(k) correspond to transverse advection modes, orthogonal to the direction of propagation. These modes do not propagate because ω_{1,2} = 0. The eigenvectors b_{f}(k) and b_{b}(k) represent forward and backward acoustic waves, which are longitudinal , and which propagate with the sound speed v given by (13).

The limit Wigner distribution matrix W (t, x, k) has the form [17]:

W (t, x, k) =∑_{τ=1}^{2} W_{ij}^{τ}(t, x, k)b^{i}(k)b^{j∗}(k)
+a_{f}(t, x, k)b_{f}(k)b_{f}^{∗}(k) + a_{b}(t, x, k)b_{b}(k)b_{b}^{∗}(k).  (15)

The first term corresponds to the non-propagating modes and may be set to zero here without any loss of generality. The last two terms correspond to forward and backward propagating sound waves. The scalar functions a_{f,b} are related by a_{f}(t, x, k) = a_{b}(t, x, −k), and a_{f} satisfies the Liouville equation

\frac{∂a}{∂t} + ∇_{k}ω · ∇_{x}a − ∇_{x}ω · ∇_{k}a = 0.  (16)

They may be interpreted as phase space energy densities since they are nonnegative (because the matrix W (t, x, k) is non-negative) and

E(x) = \frac{1}{2} ∫ dk[a_{f}(t, x, k) + a_{b}(t, x, k)] = ∫ dka_{f}(t, x, k).

The flux is given by

F = \frac{v}{2} ∫ dk[k̂a_{f}(t, x, k) − k̂a_{b}(t, x, k)] = v ∫ dkk̂a_{f}(t, x, k).  (17)

The radiative transport equation (1) arises when the density ρ and compressibility κ are random and oscillating on the scale of the wave length, so we assume they have the form

ρ → ρ(1 + \sqrt{ε}ρ_{1}(\frac{x}{ε})), κ → κ(1 + \sqrt{ε}κ_{1}(\frac{x}{ε})).

The random processes ρ_{1} and κ_{1} are mean zero space homogeneous with power spectral densities R̂_{ρρ}, R̂_{κκ}, and cross spectral density R̂_{κρ}. The limit ε → 0 is the high frequency limit since the parameter ε is the ratio of wave length to propagation distance. In (3.1) we take the ratio of correlation length to propagation distance to be of order ε also, and we take the standard deviation of the fluctuations to be of order \sqrt{ε}. It is in this scaled limit that radiative transport theory emerges. The radiative transport equation for a(t, x, k) = a_{f}(t, x, k) is

\frac{∂a}{∂t}+vk̂ · ∇_{x}a − |k|∇_{x}v · ∇_{k}a = \frac{πv^{2}|k|^{2}}{2} ∫ δ(v|k| − v|k′|)[a(k′) − a(k)]
· {(k̂ · k̂′)^{2}R̂_{ρρ}(k − k′) + 2(k̂ · k̂′)R̂_{ρκ}(k − k′) + R̂_{κκ}(k − k′)} dk′. (18)

This equation is of the form (1). The mean free path is a typical value of the ratio \frac{v}{Σ}, the speed over the total scattering cross-section. It can be thought of as the distance over which scattering by the inhomogeneities is effective. It is a length scale that can be estimated from seismic data while correlation lengths and standard deviations of parameter fluctuations are usually not observable.

The radiative transport equation (1) has been derived from equations governing particular wave motions by various authors, such as Stott [21], Watson et.al. [22], [23], [24], Barabanenkov et.al. [25], Besieris and Tappert [26], Howe [27], Ishimaru [3] and Besieris et. al. [28] with a recent survey presented in [29]. These derivations also determine the functions ω(x, k) and σ(x, k, k′) and show how a is related to the wave field. In [17], (1) and these functions are derived as a special case of a more general theory using the Wigner distribution and symmetric hyperbolic systems.

There is not a lot of mathematical work on the wave-to-transport limit, and most of it is for the Schrödinger equation with random potential. We cite here the work of Martin and Emch [30], of Spohn [31], of Dell’Antonio [32] and the recent extensive study of Ho, Landau and Wilkins [33] as well as [34]. They treat only spatially homogeneous problems but it is known how to extend the analysis to the spatially inhomogeneous case (slow x-dependent initial data and potential) [35]. A really satisfactory mathematical treatment of radiative transport asymptotics from random wave equations is lacking at present.

3.2 Transport for electromagnetic and elastic waves

Transport theory for electromagnetic and elastic waves is interesting because of wave polarization. This is important in astrophysics and is analyzed in great detail in Chandrasekhar’s treatise [36]. Coherence of polarized light persists and must be tracked, leading to a system of transport equations for the Stokes parameters that fix the state of polarization. The derivation of this system from Maxwell’s equations was first done in the early seventies, and using symmetric hyperbolic systems and Wigner distributions in [17], where the earlier papers are cited.

The main reason we wanted a general derivation of transport equations for general waves was so that we could deal with elastic waves. One can, of course, write down phenomenological equations for the transport of elastic wave energy and this was done often in the last 10-15 years [5]. The problem is that shear waves were treated like acoustic waves and the role of polarization was not accounted for correctly in the geophysics literature, even though the similarity with electromagnetic waves (Chandrasekhar’s work) should be clear. In [17] it is shown that elastic wave transport is like E&M for shear and like acoustics for compressional waves, and the two wave modes are coupled.

A simple but interesting consequence of the general derivation is the symmetry (self-adjointness) of the transport equations. This implies immediately that the only equilibrium phase space energy densities are the uniform ones (over the support of the energy surface). The spatial energy densities for the compressional P waves and the shear S waves must be in a fixed ratio to each other, which turns out to be

E_{P} = \frac{v_{S}^{3}}{2v_{P}^{3}} E_{S}

Here v_{S} is the shear speed (about 3km/sec) and v_{P} is the compressional speed (about 5km/sec). This makes the P wave energy about one tenth of the S wave energy deep in the coda of seismograms, assuming surface effects are not important so that the free space theory can be used. This is independent of what the source is and of the details of the scattering medium, as long as there is effective scattering. The asymptotic energy law is well known empirically but it did not have an explanation from first principles so we presented it in detail and related it to the seismological literature in [37]. It turns out that this kind of long time P-to-S energy equilibrium was known in connection with remote sensing with ultrasound [41].

3.3 Boundary conditions

Finding appropriate transport theoretic boundary conditions for wave propagation in the transport regime is perhaps the most pressing issue both theoretically and for the applications, in geophysics, in electromagnetics, in ultrasound and elsewhere. The problems are analytically difficult as can be seen from [20] where the relatively simple case of inhomogeneous, slowly varying deterministic media with a flat interface is considered and transport theoretic boundary conditions are derived in the high frequency limit.

There is a lot of physical and applied literature on scattering form random rough surfaces [38, 39]. The issue is to determine what is appropriate as a boundary or interface condition for radiative transport equations. As with polarization, interfaces are a source of coherence in an otherwise incoherent scattering process. So they must be treated carefully to avoid oversimplifications. In [40] we consider acoustic reflection and transmission by a flat interface and derive transport theoretic boundary conditions, but a lot more has to be done here, including the derivation of boundary conditions for E&M and elastic wave transport.

3.4 The diffusive regime

It is well known, primarily from studies that originated in neutron scattering and reactor theory, that when the propagation distance in the transport regime is large compared to the mean free path a simpler diffusion theory emerges. In some seismic propagation problems the mean free path is 20-30 km but propagation is over 1000 km and more. So it is quite clear that a diffusion approximation for the transport equations is called for. We know how to do this when there are no boundaries present [17], even with polarization for E&M and elastic waves.

The problem is that the crustal wave guide is 30-40 km deep and it is not clear how to use the diffusion approximation, or even how to decide if it should be used at all. But the mathematical problem of finding asymptotic boundary and interface conditions in the diffusive regime is interesting, quite delicate analytically and potentially very useful [42]. In radar scattering from clutter, the diffusive transport theory is very likely the most appropriate one to use for wavelengths in the 10 cm to 1 m range, for example.

3.5 Parameter identification and inverse problems

Parameter identification for radiative transport has received relatively little attention in geophysics [5]. In light propagation through the atmosphere the situation is, of course, very different if only because the measurements that can be made are very different. The recent activity in diffusive tomography [43] should eventually find applications in geophysics as well, but there are many difficult problems that must be settled along the way, such as getting the right transport theoretic boundary conditions.

4 The localization regime

I will review briefly reflection of acoustic plane wave pulses normally incident on a randomly layered half space, z < 0, with z the direction of the layering [45, 46]. A good reference for deterministic wave propagation in layered media is Brekhovskikh’s book [44]. It is in randomly layered media that wave localization is dominant. I will describe it in the time domain, for pulses, because this is the most interesting case in geophysical wave propagation, in reflection seismology and elsewhere. It is also not treated much in the mathematical or physical literature specialized to localization problems, and the simple intuition that most specialists have for time harmonic, one dimensional wave localization is not quite adequate for pulses. This was pointed out some time ago [50].

Radiative transport theory is not, of course, valid for randomly layered media. This was also considered long ago in connection with wave guides and optical fibers [47]. But it is not well understood in applied fields, even today, and papers appear occasionally that attempt to ‘derive’ radiative transport equations for propagation in layered media. I do not mean here three dimensional radiative transport in plane parallel structures. I mean random layering. If radiative transport were valid in this case, then the differential scattering cross-section would be singular, concentrated in only two (in the simplest case) directions, up and down or forwards and backwards propagation.

In the long paper [48] we deal in detail with the point source case, that is, the propagation of an acoustic pulse generated by a point source over a layered random medium. Here I will describe only the reflection of acoustic plane wave pulses.

4.1 Pulse reflection from randomly layered media

The acoustic pressure p(t, z) and velocity u(t, z) satisfy the continuity and momentum equations

\frac{1}{K} p_{t} + u_{z} = 0

ρu_{t} + p_{z} = 0  (19)

Here ρ is the material density and K the bulk modulus. As in [48] we assume for simplicity that the density has no random variation

ρ(z) = \{ \begin{matrix}ρ_{0}, & z > 0, \\ ρ_{1}, & z < 0 \end{matrix}  (20)

with ρ_{0} and ρ_{1} constants. For the bulk modulus we assume that

K^{−1}(z) = \{ \begin{matrix}K_{0}^{−1}, & z > 0, \\ K_{1}^{−1}(z) (1 + ν(\frac{z}{ε^{2}} )) , & z < 0 \end{matrix}  (21)

with K_{0} a constant, K_{1}(z) a smooth deterministic function of z and ν(s) a bounded stationary random function with mean zero, representing the fluctuations in K^{−1}. Note that they vary on the scale ϵ^{2}, where ϵ is a small parameter. If z is measured in kilometers and the fluctuations vary on the scale of a few meters then a value of ϵ around 0.05 captures the scale separation we wish to model. We assume that the random function ν(s) has a correlation length of order one so that the correlation length of ν(z/ϵ^{2}) is of order ϵ^{2} in kilometers (about 2.5 meters for ϵ = 0.05). The mean sound speed c is given by

c(z) = \{ \begin{matrix}c_{0} = \sqrt{\frac{K_{0}}{ρ_{0}}} , & z > 0 \\ c̄(z) = \sqrt{\frac{K_{1}(z)}{ρ_{1}}}, & z < 0 \end{matrix}  (22)

Note that the fluctuations in the sound speed are not assumed to be small. The estimation of the vertical correlation length of the inhomogeneities in the lithosphere from well-log data is considered in [51]. They found that 2-3 m is a reasonable estimate of the correlation length of the fluctuations in sound speed.

For t < 0 a normally incident plane wave solution in z > 0 has the form

u(t, z) = \frac{1}{\sqrt{ϵ}} \frac{1}{\sqrt{ρ_{0}c_{0}}} f (\frac{t + z/c_{0}}{ϵ} )

p(t, z) = −\frac{1}{\sqrt{ϵ}}\sqrt{ρ_{0}c_{0}}f (\frac{t + z/c_{0}}{ϵ} )  (23)

Here f is the pulse shape function which is assumed to vanish for negative arguments and to have support that is of order one in the macroscopic t units that are seconds. With ϵ = 0.05, the pulse width is about 50 msec or, with a speed of 3 km/sec, 150 meters. The multiplicative factor 1/\sqrt{ϵ} in (23) makes the total energy of the incident plane wave pulse independent of ϵ. Continuity of p and u at the interface z = 0 makes (19) and (23) a complete problem. We are interested in p(t, 0) or u(t, 0) for t > 0, the pressure or velocity at the interface, and this involves the solution of a complicated random scattering problem because of the form (21) of K^{−1}.

4.2 Scale separation

The scaling that we have chosen, and the asymptotic limit ε → 0 that we will consider, models well problems in reflection seismology and is quite different from transport theoretic scaling. The main differences are that the fluctuations are not assumed to be small and the typical wavelength of the probing pulse (150 m) is small with respect to the probing depth (5 km, say) but large compared to the correlation length (2-3 m). The parameter ε is then the ratio of the (typical) wave length to propagation depth, as well as the ratio of correlation length to wave length. This is a particularly interesting scaling limit mathematically because it is a high frequency limit with respect to the large scale variations of the medium that we want to detect, but it is a low frequency limit with respect to the fluctuations, whose effect acquires a canonical form independent of details.

Is this model realistic and can it be used effectively? One argument that can be made against it is this: There is no real scale separation in sound speed fluctuations, as one can see from well-logs [52], so this neat way of dealing with fluctuations, background and probing pulse cannot possibly be right, even if it can handle large fluctuations. Another is that perfectly layered random media are an unacceptable idealization.

Regarding scale separation, it is fair to say that the scope of the analytical theory that has been developed, and is described briefly here, is well beyond anything that could be expected from any theory that deals with strong fluctuations in a serious way. Radiative transport theory is more robust because the fluctuations are assumed to be small, and then it is not necessary to have scale separation (correlation lengths and wave lengths are comparable). Moreover, the analytical tools that emerge from the asymptotic scale separation theory are far more flexible and robust than the crude thinking with scales implies. Discontinuities and imperfections that are comparable to the pulse width can be handled by the theory and do not make it unusable. The problem is that the theory is not easy to follow, it is analytically difficult to implement and not nearly enough has been done to test it in situations that push against the scale separation assumptions. The statistical analysis of well-log data that was done in [51], that produced the estimate of 2-3 m for the correlation length of the sound speed fluctuations, is quite thorough, but perhaps more can be done here also.

The modifications to the theory that are needed to account for imperfect layering are far more important than anything missed by scale separation asymptotics. This goes back to the localization-delocalization transition that I have mentioned several times already. It remains a big gap in our understanding of wave propagation in random media.

4.3 Localization regime asymptotics

We will consider the reflected pressure p_{refl}(t, 0), at z = 0 and t > 0, which is the total pressure minus the incident pressure (23). After a time of order ϵ, the duration of the incident pulse, the two are the same. Of particular interest is the two-time reflected average pressure intensity.

I(t, t̄) = \frac{1}{ρ_{0}c_{0}} < p_{refl}(t + \frac{ϵt̄}{2}, 0)p_{refl}(t − \frac{ϵt̄}{2}, 0) >  (24)

with the angular brackets denoting statistical average. The factor 1/ρ_{0}c_{0} is a normalization.

For simplicity, we will assume in the sequel that there is no macroscopic discontinuity at z = 0 so that ρ_{0} = ρ_{1} and K_{0}^{−1} = K_{1}^{−1}(0).

Note that the time offset in (24) is proportional to the pulse width ϵ. The reason for this is that for time offsets of more than a few pulse widths the reflected signals are essentially uncorrelated. Moreover, in the absence of discontinuities in the medium, < p_{refl}(t, 0) > is essentially zero except for a time of order ϵ near t = 0 when the reflection from the interface z = 0 is felt. That is, there is no coherent backscattering. We formulated a scattering problem where the quantity of interest is as directly related to the medium fluctuations as is possible.

Fix a t > 0, not close to zero, and a small ϵ. Since I(t, t̄) is essentially zero for large t̄ we can introduce its (essentially local) Fourier transform

Λ(t, ω)|f̂(ω)|^{2} = ∫ e^{iωt̄}I(t, t̄)dt̄  (25)

in which Λ is the normalized local power spectral density. The normalization is |f̂(ω)|^{2} with f̂(ω) the Fourier transform of the pulse shape function f(t). The two-time intensity function can be written as

I(t, t̄) = \frac{1}{2π} ∫ |f̂(ω)|^{2}Λ(t, ω) e^{−iωt̄}dω  (26)

The main thrust of our theoretical work in [45, 48, 46] is that in the limit ϵ → 0 the local power spectral density can be calculated by solving a system of partial differential equations where

Λ(t, ω) = W^{1}(0, t, ω)  (27)

and the W^{N}(z, t, ω), N ≥ 0 satisfy the equations

\frac{∂W^{N}}{∂z} + \frac{2N}{c̄(z)} \frac{∂W^{N}}{∂t} − \frac{2αω^{2}N^{2}}{c̄^{2}(z)} \{W^{N+1} − 2W^{N} + W^{N−1}\} = 0  (28)

for −L < z ≤ 0, with

W^{N}(−L, t, ω) = δ(t)δ_{N,0}  (29)

Here the mean sound speed c̄(z) is given by (22) and α > 0 is the noise intensity level of the fluctuations

α = \frac{1}{4} ∫_{0}^{∞} < ν(s)ν(0) > ds  (30)

The length L is arbitrary, provided that for any given t > 0 for which we want to calculate Λ(t, ω) it satisfies

L > c_{max} \frac{t}{2}  (31)

with c_{max} the maximum speed c̄(z) in z ≤ 0. Because of the hyperbolic nature of the equations (28) it is easy to see (and explained in the references) that the choice of L satisfying (31) does not affect Λ(t, ω) given by (27).

4.4 Time domain localization

There is no quick and simple way to explain the result (27)-(31) that relates the local power spectral density Λ(t, ω; c̄(.)), the mean sound speed profile c̄(z) and the noise intensity level α. But we will now make several remarks that will help explain the nature of this relationship.

From the definition (24) and (25) it is clear that Λ(t, ω) is a local Fourier transform but it is not necessarily positive as it would have to be if p_{refl}(t, 0), the reflected pressure, were a stationary process in t so that I(t, t̄) were independent of t. However, in the limit ϵ → 0, and hence when ϵ is small, Λ(t, ω) given by (27)- (31) is indeed positive. For a general profile c̄(z) it cannot be computed explicitly but for c̄(z) = c̄, a constant, it has the form

Λ(t, ω) = \frac{\frac{αω^{2}}{c̄}}{(1 + \frac{αω^{2}}{c̄} t)^{2}}  (32)

In terms of the localization length [49] at frequency ω

l(ω) = \frac{c̄^{2}}{2αω^{2}}  (33)

we can write (32) in the form

Λ(t, ω) = \frac{1}{2} \frac{c̄l(ω)}{(l(ω) + \frac{c̄t}{2})^{2}}  (34)

As shown in [49] and the many references cited there, the localization length at frequency ω is a measure of the depth of penetration of a time harmonic wave with this frequency into a randomly layered medium with uniform sound speed c̄ and noise level α for the fluctuations. Wave energy does not penetrate much below this length. If T (L, ω) is the time harmonic transmission coefficient for a randomly layered medium of width L, with ω the frequency of the incident plane wave, then

lim_{L→∞} \frac{1}{L} log |T (L, ω)| = \frac{−1}{l(ω)}

with probability one. This defines the localization length l(ω) > 0, which is always positive for a large class of random media. It cannot be computed explicitly but in the low frequency limit it has the form (33). The lower the frequency the deeper the penetration of the waves into the randomly layered medium.

In the time domain, the normalized local power spectral density of the reflected signal at a fixed time t, Λ(t, ω) in (34), has a maximum ω_{max} = ω_{max}(t) that depends on time. From (34) the maximum is calculated to be

ω_{max} = \sqrt{\frac{c̄}{αt}}  (35)

In a more physical way [49] we can say that the maximum of the local power spectral density at time t occurs for that ω = ω_{max} for which

l(ω_{max}) = \frac{c̄t}{2}  (36)

Thus, for the frequency for which the localization length equals the mean distance traveled into the medium, we have the maximum contribution to the noise spectrum of the reflected signal. This is a stochastic resonance relation that identifies precisely the main source of noise in the reflected signals.

It is because of wave localization and its manifestations in the time domain described above that signals reflected by randomly layered media are so noisy. From (32) we find, by integrating over ω, that the envelope of the root mean square of the reflected pulse is of the form constant×t^{−3/4}. Thus, the fluctuations in the reflected signal decay very slowly, indicating that a great deal of multiple scattering is taking place and that wave localization is dominant.

We can interpret (28) as a hierarchy of equations for moments associated with the scattering problem [48]. The infinite hierarchy (N ≥ 0 in (28)) indicates that the second moment that we are interested in (I of (24) or Λ of (25)) cannot be computed separately from all higher moments (the W^{N}, N > 2 in (28)). This is another manifestation of localization.

When we use the parameters of section 4.2 that are typical in reflection seismology we find that the minimum localization length occurs in the 20-30 Hz regime and is about 15-20 km [51]. This means that random inhomogeneities will effectively prevent probing below this depth because all the wave energy is reflected to the surface by multiple scattering.

What is missing at present is a more general theory that allows us to compute the changes in the one dimensional theory that occur when small three dimensional inhomogeneities are introduced into the model. We need a more general theory that lets us go from localization to transport as the random layering is reduced and isotropic inhomogeneities replace it.

4.5 Statistical inverse problems

I will describe briefly how the mean sound speed profile c̄(z) can be estimated from observations of p_{refl}(t, 0) or

R_{f}(t) = \frac{1}{\sqrt{ρ_{0}c_{0}}} p_{refl}(t, 0)  (37)

in which dependence of the pulse shape function f is indicated. The inversion strategy is based on one more fact about the reflected signal R_{f}(t), in addition to (27)-(31). It is that as ϵ tends to zero R_{f}(t) becomes approximately a Gaussian process. It has not been possible to prove this so far but there are some good heuristic indications that it is true [45] and extensive numerical simulations corroborate it very well [48]. From the Gaussian property of R_{f}(t) we conclude that

\frac{1}{|f̂(ω)|^{2}} ∫ e^{iωt̄}R_{f} (t + \frac{ϵt̄}{2}) R_{f} (t − \frac{ϵt̄}{2}) dt̄= Λ̂(t, ω)  (38)

is approximately, when ϵ is small, an exponential random variable with mean Λ(t, ω) given by (27)-(31), when c̄(z) is known. Moreover, for distinct 0 < t_{1} < t_{2} < ... < t_{N_{t}} and 0 < ω_{1} < ω_{2} < ... < ω_{N_{f}}, where N_{t} and N_{f} are integers, the random variables {Λ̂(t_{j}, ω_{l})} are independent with exponential distribution having mean {Λ(t_{j}, ω_{l})} .

The inversion strategy is now this: Depending on the available data, fix a set of time points {t_{j}} and frequencies points {ω_{l}} as above. For each realization of R_{f}(t) that is available, estimate Λ̂(t_{j}, ω_{l}) from (38). This is actually a very delicate step that must be done carefully as we discuss in [48], Appendix E. Then form

O(c̄) = ∏_{realiz} ∏_{j=1}^{N_{t}} ∏_{l=1}^{N_{f}} \frac{e^{−Λ̂(t_{j},ω_{l})/Λ(t_{j},ω_{l};c̄(.))}}{Λ(t_{j}, ω_{l}; c̄(.))}  (39)

where the first product is over different independent realizations. This is the likelihood functional for the estimates Λ̂, given a known mean speed profile c̄(z). We now choose c̄(z) in order to maximize this functional. This is a rather usual maximum likelihood estimation except that now the maximization must be done over the profiles c̄(.) which in turn determine Λ(t, ω; c̄(.)) in (39) via the partial differential equations (28)-(29) and the relation (27).

The most convenient way to solve the maximization problem for (39), and thus estimate c̄(z), is to assume that it is piece-wise linear over a few macroscopically large layers and then maximize O over a finite set of speeds c̄^{1}, c̄^{2}, ..., c̄^{N_{z}}. These speeds are approximations of c̄(z) at successively larger depths numbered from 1 to N_{z}. Moreover, because of the hyperbolic nature of (28)-(29) the maximization can be done one layer at a time with increasing depth. This avoids the difficult problem of finding the maximum of a complicated function of several variables. Physically this layer peeling process makes sense because there is a direct relation between the sound speed profile up to a certain depth and the smallest time before which the rest of the medium is not felt in the reflected signal R_{f}(t).

Of course we need a lot of independent realizations to get reasonable results and this is unrealistic in a geophysical context. But it is important in principle to make this strategy work and amazingly enough it does [50], [48]). It is amazing because we are trying to determine the smooth, mean speed profile from the reflected signals that are swamped by fluctuations due to multiple scattering. The computational and other implementation details are described in [48]).

Could we do this kind of inversion from extremely noisy reflections if we only had one realization? Yes, if we have reflection measurements at different offsets (distances from the source) on the interface, generated by a point source over a randomly layered medium [53]. This is a very difficult problem that requires a great deal of numerical computation. The inversion is not as good as in the plane wave case (with many realizations) but it is reasonably good and, in any case, it shows that the strategy does work. But improving the results requires very careful attention to a host of implementation issues that can be settled only empirically, by trial and adjustment, at present.

An interesting discussion of reflections from time reversed reflections, their statistical properties and their relation in turn to the hierarchy of moments equations (28) is given by Clouet and Fouque [54]. This work should have important applications in statistical inverse problems of geophysical interest.

Another application of time domain localization asymptotics is to surface water waves over a rough bottom [55].

4.6 Reflection and transmission of time harmonic plane elastic waves

We have described a variety of results for acoustic pulse reflection from randomly layered media, emphasizing time domain effects. For geophysical applications we must also consider elastic waves in randomly layered media. The analytical difficulties in extending the theory that we briefly described above to the elastic case are enormous, mainly because there are two wave modes, P and S waves, that are coupled by the inhomogeneities. In [56] we extended the scale-separation asymptotic theory to time harmonic, obliquely incident elastic plane waves. We calculate in detail mode coupling in reflection and transmission, with various kinds of interfacial discontinuities. It is surprising that so many things can be calculated analytically and in such detail, given the complexity of the problem.

However, despite considerable efforts we have not been able to extend the results to the time domain. The hierarchy of moment equations that we used in the analysis of acoustic pulse reflection does not seem to work for elastic wave pulse reflection. The analysis of reflections for elastic wave pulses generated by a point source, the analog of the analysis carried out in [48] for acoustic waves, seems to be out of reach at present.

4.7 Pulse stabilization and imaging

We have focused mostly on reflection in the time domain because the bulk of the measurements that can be made in geophysics, in nondestructive testing with ultrasound and elsewhere are surface measurements. However, transmission is also important as is the analysis of reflections from imbedded discontinuities in a randomly layered medium. The vicinity of the front of the pulse, or the vicinity of first arrival from the discontinuity, has an interesting structure that can be analyzed in considerable detail. This is called the O’Doherty-Anstey theory because it was first discussed by these two geophysicists in the early seventies [57]. The main point is that if the fluctuations are weak and the pulse is followed with its random speed, then it will appear to stabilize (not fluctuate) and become broader as it advances into the medium. This is discussed in detail in [48] where many other papers are cited.

What if the fluctuations are not weak, and we have scale separation as described above? Do we have an O’Doherty-Anstey theory? This question was answered in [58, 59] by overcoming what was the main obstacle before: finding the right random speed with which to center the advancing pulse. The fact that the advancing pulse spreads and loses energy (to fluctuations in its coda) is not so surprising and is true for general random media, not only layered media, although the fluctuations must be weak. What is surprising, and not generally known or anticipated in the geophysics literature, is that in the case of large fluctuations the centering speed is not the local random speed but a function of it, and the centered pulse stabilizes with probability one with minimal spreading (relative to other centerings).

In [60], Solna shows how this theory can be used to improve the resolution of discontinuity identification in a random medium. He also extends the O’Doherty-Anstey theory to a class of locally layered random media, that is, he allows for slow horizontal variations.

References

[1] P.W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev., 109 (1958), pp. 1492-1505.

[2] J.B. Keller, Stochastic Equations and Wave Propagation in Random Media, Proc. Symp. Appl. Math., 16, 145–170, Am. Math. Soc., McGraw–Hill, New York, 1964.

[3] A.Ishimaru, Wave propagation and scattering in random media, vol. I,II, Academic Press, New York, 1978.

[4] K.Aki and R.S.Wu, Multiple scattering and energy transfer of seismic waves – separation of scattering effect from intrinsic attenuation –II. Application of the theory to Hindu Kush region, in Seismic wave scattering and attenuation, vol. I (R.S.Wu and K.Aki, eds), pp. 49-80, 1988.

[5] M.C. Fehler and H. Sato, Seismic Wave Propagation and Scattering in the Heterogeneous Earth, AIP Press, 1998.

[6] P. Sheng, Introduction to wave scattering, localization, and mesoscopic phenomena, Academic Press, 1995.

[7] A. Figotin and L.A. Pastur, Spectra of Random and Almost-Periodic Operators, Springer-Verlag, 1992.

[8] R. Knapp, G. Papanicolaou and B.S. White, Transmission of waves by a nonlinear random medium, J. Stat. Physics, 63, pp. 567-584, 1991.

[9] J.C. Bronski, Nonlinear Wave Propagation in a Disordered Medium J. Stat. Phys., 92, 5/6, 1998.

[10] J.C. Bronski Nonlinear Scattering and the Analyticity Properties of Solitons, Journal of Nonlinear Science, 8, pp. 161-182, 1998.

[11] J.C. Bronski, D.W. McLaughlin and M.J. Shelley, Stability of Time-Harmonic Localized Solutions in a One Dimensional Disordered Nonlinear Medium, J. Stat. Phys., 88, 1997.

[12] J. Garnier, Transmission of solitons through random media, to appear in SIAM, J. Appl. Math.

[13] J. Garnier, PhD Thesis, Ecole Polytechnique, 1996.

[14] J.F. Claerbout, Imaging the Earth’s Interior, Blackwell, Boston, 1985.

[15] N. Bleistein, J.K. Cohen and J.W. Stockwell Jr., Mathematics of Multidimensional Seismic Inversion, pp. 1-260, Colorado School of Mines, 1998.

[16] W.W. Symes, Mathematics of Reflection Seismology, pp. 1-85, Rice University, 1995.

[17] L.Ryzhik, G.Papanicolaou and J.Keller, Transport Equations for Elastic and Other Waves in Random Media, Wave Motion 24, pp. 327-370, 1996.

[18] M.C.W. van Rossum and Th.M.Nieuwenhuizen, Multiple scattering of classical waves: from microscopy to mesoscopy and diffusion, to appear in Phys. Rev. E.

[19] P.Gérard, P.Markowich, N.Mauser and F.Poupaud, Homogenization limits and Wigner transforms, Comm. Pure Appl. Math., 50, pp. 323-380, 1996.

[20] L.Ryzhik, G.Papanicolaou and J.Keller, Transport equations in a half space, Comm. PDE’s, 22, pp. 1869-1911, 1997.

[21] P.Stott, A transport equation for the multiple scattering of electromagnetic waves by a turbulent plasma, Jour. Phys. A, 1, pp. 675-689, 1968.

[22] K.Watson and J.L.Peacher, Doppler shift in frequency in the transport of electromagnetic waves in an underdense plasma, Jour. Math. Phys., 11, pp. 1496-1504, 1970.

[23] K.Watson, Multiple scattering of electromagnetic waves in an underdense plasma, Jour. Math. Phys., 10, pp. 688-702, 1969.

[24] C.W.Law and K.Watson, Radiation transport along curved ray paths, Jour. Math. Phys., 11, pp. 3125-3137, 1970.

[25] Yu.Barabanenkov, A.Vinogradov, Yu.Kravtsov and V.Tatarskii, Application of the theory of multiple scattering of waves to the derivation of the radiative transfer equation for a statistically inhomogeneous medium, Radiofizika, 15, pp. 1852-1860, 1972. English translation pp. 1420-1425.

[26] I.M.Besieris and F.D.Tappert, Propagation of frequency modulated pulses in a randomly stratified plasma, Jour. Math. Phys., 14, pp. 704-707, 1973.

[27] M.S.Howe, On the kinetic theory of wave propagation in random media, Phil. Trans. Roy. Soc. Lond., 274, pp. 523-549, 1973.

[28] I.M. Besieris, W. Kohler and H. Freese, A transport-theoretic analysis of pulse propagation through ocean sediments, J. Acoust. Soc. Am., 72, pp. 937-946, 1982.

[29] Yu.Barabanenkov, Yu.Kravtsov, V.Ozrin and A.Saichev, Enhanced backscattering in optics, Progress in Optics, 29, pp. 67-190, 1991.

426 George Papanicolaou

[30] P. Martin and G. Emch, A rigorous model sustaining van Hove’s phenomenon, Helv. Phys. Acta, 48, p. 59, 1967.

[31] H. Spohn, Derivation of the transport equation for electrons moving through random impurities, Jour. Stat. Phys., 17, pp. 385-412, 1977.

[32] G.Dell’Antonio, Large time small coupling behavior of a quantum particle in a random potential, Annals of Inst. H. Poincare, section A, 39, p. 339, 1983.

[33] T. Ho, L. Landau and A. Wilkins, On the weak coupling limit for a Fermi gas in a random potential, Rev. Math. Phys., 5, pp. 209-298, 1993.

[34] L. Erdos and H.T. Yau, preprint, 1997.

[35] A.Fannjiang, unpublished notes, 1991.

[36] S.Chandrasekhar, Radiative transfer, Dover, New York, 1960.

[37] G.Papanicolaou, L.Ryzhik and J.Keller, On the stability of the P to S energy ratio in the diffusive regime, Bulletin of the Seismological Society of America, 86, pp. 1107-1115, 1996.

[38] J.A. DeSanto and G.S. Brown, Analytical techniques for multiple scattering from rough surfaces, in Progress in Optics XXIII (E. Wolf, ed.), pp. 1-62, Elsevier, Amsterdam, 1986.

[39] A.Voronovich, Wave Scattering from Rough Surfaces, Springer Verlag, 1991.

[40] G. Bal, J.B. Keller, G. Papanicolaou and L. Ryzhik, Transport Theory for Acoustic Waves with Reflection and Transmission at Interfaces, preprint, 1998.

[41] R.L.Weaver, On diffuse waves in solid media, Jour.Acoust.Soc.Am., 71, pp. 1608-1609, 1982.

[42] G. Bal, PhD Thesis, University of Paris VI, 1997.

[43] O. Dorn, Diffuse tomography, PhD Thesis, University of Münster, 1997.

[44] L.Brekhovskikh, Waves in Layered Media, Academic Press ,1980.

[45] R. Burridge, G. Papanicolaou, P. Sheng, and B. White, Probing a random medium with a pulse, SIAM J. Appl. Math., 49, pp. 582-607, 1989.

[46] G. Papanicolaou and S. Weinryb, A functional limit theorem for waves reflected by a random medium, Applied Mathematics and Optimization, 30, pp. 307-334, 1994.

[47] W. Kohler and G. Papanicolaou, Power statistics for wave propagation in one dimension and comparison with radiative transport theory I,II, J. Math. Phys., 14, pp. 1733-1745, 1973, and 15, pp. 2186-2197, 1974.

Documenta Mathematica · Extra Volume ICM 1998 · I · 403–427

[48] M. Asch, W. Kohler, G. Papanicolaou, M. Postel and B. White, Frequency Content of Randomly Scattered Signals, SIAM Review, 33, pp. 526-629, 1991.

[49] P. Sheng, B. White, Z.Q. Zhang and G.Papanicolaou, Wave localization and multiple scattering in randomly-layered media, in ‘Scattering and Localization of Classical Waves in Random Media,’ Ping Sheng, editor, World Scientific, pp. 563-619, 1990.

[50] B. White, P. Sheng, M. Postel, and G. Papanicolaou, Probing through cloudiness: Theory of statistical inversion for multiply scattered data, Phys. Rev. Lett., 63, pp. 2228-2231, 1989.

[51] B. White, P. Sheng, and B. Nair, Localization and backscattering spectrum of seismic waves in stratified lithology, Geophysics, 55, pp. 1158-1165, 1990.

[52] F. Herrmann, A scaling medium representation, a discussion of well-logs, fractals and waves, PhD Thesis, Delft University, 1997.

[53] M. Asch, W. Kohler, G. Papanicolaou, M. Postel and B. White, Statistical inversion from reflections of spherical waves by a randomly layered medium, Waves in Random Media, 6, pp. 293-334, 1996.

[54] J.F. Clouet and J.P. Fouque, A time-reversal method for an acoustical pulse propagating in randomly layered media, Wave Motion, 25, pp. 361-368, 1997.

[55] A. Nachbin and G. Papanicolaou, Water waves in shallow channels of rapidly varying depth, JFM, 241, pp. 311-332, 1992.

[56] W. Kohler, G. Papanicolaou and B. White, Localization and mode conversion for elastic waves in randomly layered media I,II, Wave Motion, 23, pp. 1-22 and 181-201, 1996.

[57] R.F. O’Doherty and N.A. Anstey, Reflections on amplitudes, Geophysical Prospecting, 19, pp. 430-458, 1971.

[58] J.F. Clouet and J.P. Fouque, Spreading of a pulse traveling in random media, Annals of Applied Probability, 4, pp. 1083-1097, 1994.

[59] R. Burridge, P. Lewicki and G. Papanicolaou, Pulse stabilization in a strongly heterogeneous layered medium, Wave Motion, 20, pp. 177-195, 1994.

[60] K. Solna, PhD Thesis, Stanford University, 1997.

George Papanicolaou

Department of Mathematics

Stanford University Stanford CA 94305 papanico@math.stanford.edu URL: http://georgep.stanford.edu 428

Documenta Mathematica · Extra Volume ICM 1998 · I