# Hierarchical model for the scale-dependent velocity of seismic waves

Elastic waves of short wavelength propagating through the upper layer of the Earth appear to move faster at large separations of source and receiver than at short separations. This scale dependent velocity is a manifestation of Fermat’s principle of least time in a medium with random velocity fluctuations. Existing perturbation theories predict a linear increase of the velocity shift with increasing separation, and cannot describe the saturation of the velocity shift at large separations that is seen in computer simulations. Here we show that this long-standing problem in seismology can be solved using a model developed originally in the context of polymer physics. We find that the saturation velocity scales with the four-third power of the root-mean-square amplitude of the velocity fluctuations, in good agreement with the computer simulations.

Seismologists probe the internal structure of the Earth by recording the arrival times of waves created by a distant earthquake or explosion. Systematic differences between studies based on long and short wavelengths $`\lambda `$ have been explained in terms of a scale dependence of the velocity at short wavelengths. The velocity obtained by dividing the separation $`L`$ of source and receiver by the travel time $`T`$ increases with increasing $`L`$, because — following Fermat’s principle — the wave seeks out the fastest path through the medium (see Fig. 1). This search for an optimal path is more effective for large separations, hence the apparent increase in velocity on long length scales. It is a short-wavelength effect, since Fermat’s principle breaks down if the width $`\sqrt{L\lambda }`$ of the first Fresnel zone becomes comparable to the size $`a`$ of the heterogeneities. The scale-dependent velocity of seismic waves was noted by Wielandt more than a decade ago, and has been studied extensively by geophysicists.

A rather complete solution of the problem for small $`L`$ was given by Roth, Müller, and Snieder, by means of a perturbation expansion around the straight path. The velocity shift $`\delta v=v_0(1v_0T/L)`$ (with $`v_0`$ the velocity along the straight path) was averaged over spatially fluctuating velocity perturbations with a Gaussian correlation function (having correlation length $`a`$ and variance $`\epsilon ^2v_0^2`$, with $`\epsilon 1`$). It was found that $`\delta vv_0\epsilon ^2L/a`$ increases linearly with $`L`$. Clearly, this increase in velocity can not continue indefinitely. The perturbation theory should break down when the root-mean-square deviation $`\delta x\epsilon a(L/a)^{3/2}`$ of the fastest path from the straight path becomes comparable to $`a`$. Numerical simulations show that the velocity shift saturates on length scales greater than the critical length $`L_\mathrm{c}a\epsilon ^{2/3}`$ for the validity of perturbation theory. A theory for this saturation does not exist.

It is the purpose of this article to present a non-perturbative theory for this seismological problem, by making the analogy with a problem from polymer physics. The problem of the velocity shift in a random medium belongs to the class of optimal path problems that has a formal equivalence to the directed polymer problem. The mapping between these two problems relates a wave propagating through a medium with velocity fluctuations to a polymer moving in a medium with fluctuations in pinning energy. The travel time of the wave between source and receiver corresponds to the energy of the polymer with fixed end points. At zero temperature the configuration of the polymer corresponds to the path selected by Fermat’s principle. (The restriction to directed polymers, those which do not turn backwards, becomes important for higher temperatures.) There exists a simple solvable model for directed polymers, due to Derrida and Griffiths, that has remained unnoticed in the seismological context. Using that model we can go beyond the breakdown of perturbation theory and describe the saturation of the velocity shift on large length scales.

Hierarchical model   
We follow a recursive procedure, by which the probability distribution of travel times is constructed at larger and larger distances, starting from the perturbative result at short distances. At each iteration we compare travel times from source to receiver along two branches, choosing the smallest time. A branch consists of two bonds, each bond representing the length scale of the previous step. This recursive procedure produces the lattice of Fig. 2, called a hierarchical lattice. The lattice in this example represents a two-dimensional system, since at each step the length is doubled while the number of bonds is increased by a factor of four. For the three-dimensional version one would compare four branches at each step (each branch containing two bonds), so that the total number of bonds would grow as the third power of the length. Since most of the simulations have been done for two-dimensional systems, we will consider that case in what follows.

To cast this procedure in the form of a recursion relation, we denote by $`p_k(T)`$ the distribution of travel times at step $`k`$. One branch, consisting of two bonds in series, has travel time distribution

$$q_k(T)=_0^{\mathrm{}}𝑑T^{}p_k(T^{})p_k(TT^{}),$$
(1)   
assuming that different bonds have uncorrelated distributions. To get the probability distribution at step $`k+1`$ we compare travel times of two branches,

$`p_{k+1}(T)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑T^{}{\displaystyle _0^{\mathrm{}}}𝑑T^{\prime \prime }\delta \mathbf{\left(}T\mathrm{min}(T^{},T^{\prime \prime })\mathbf{\right)}q_k(T^{})q_k(T^{\prime \prime })`$ (2)   
$`=`$ $`2q_k(T){\displaystyle _T^{\mathrm{}}}𝑑T^{}q_k(T^{}).`$ (3)   
We start the recursion relation at step $`0`$ with the distribution $`p_0(T)`$ calculated from perturbation theory at length $`L_c`$. Iteration of eq. (3) then produces the travel time distribution $`p_k(T)`$ at length $`L=2^kL_c`$.

Equation (3) is a rather complicated non-linear integral equation. Fortunately, it has several simplifying properties. One can separate out the mean $`T_0`$ and standard deviation $`\sigma _00`$ of the starting probability distribution, by means of the $`k`$-dependent rescaling $`\tau =(T2^kT_0)/\sigma _0`$. The recursion relation (3) is invariant under this rescaling, which means that we can restrict ourselves to starting distributions $`\stackrel{~}{p}_0(\tau )=\sigma _0p_0(\sigma _0\tau +T_0)`$ having zero mean and unit variance. This is the first simplification. After $`k`$ iterations the mean $`\stackrel{~}{m}_k`$ and standard deviation $`\stackrel{~}{\sigma }_k`$ of the rescaled distribution $`\stackrel{~}{p}_k(\tau )`$ yield the mean $`T_k`$ and standard deviation $`\sigma _k`$ of the original $`p_k(T)`$ by means of

$$T_k=\sigma _0\stackrel{~}{m}_k+2^kT_0,\sigma _k=\sigma _0\stackrel{~}{\sigma }_k.$$
(4)   
The second simplification is that for large $`k`$, the recursion relation for $`\stackrel{~}{p}_k(\tau )`$ reduces to

$$\stackrel{~}{p}_{k+1}(\tau )=\frac{1}{2}\alpha \stackrel{~}{p}_k\left(\frac{1}{2}\alpha \tau +\beta \stackrel{~}{\sigma }_k+(1\alpha )\stackrel{~}{m}_k\right),$$
(5)   
with universal constants $`\alpha =1.627`$ and $`\beta =0.647`$. Under the mapping (5), the mean and standard deviation evolve according to

$$\stackrel{~}{m}_{k+1}=2\stackrel{~}{m}_k2\beta \stackrel{~}{\sigma }_k/\alpha ,\stackrel{~}{\sigma }_{k+1}=2\stackrel{~}{\sigma }_k/\alpha .$$
(6)   
The solution of this simplified recursion relation is

$$\stackrel{~}{m}_k=\frac{2^k\beta }{\alpha 1}(A\alpha ^kB),\stackrel{~}{\sigma }_k=2^kA\alpha ^k.$$
(7)   
The coefficients $`A`$ and $`B`$ are non-universal, depending on the shape of the starting distribution $`\stackrel{~}{p}_0`$. For a Gaussian $`\stackrel{~}{p}_0`$ we find $`A=0.90`$, $`B=0.95`$, close to the values $`A=1`$, $`B=1`$ that would apply if eq. (6) holds down to $`k=0`$. For a highly distorted bimodal $`\stackrel{~}{p}_0`$ we find $`A=0.71`$, $`B=0.88`$. We conclude that $`A`$ and $`B`$ depend only weakly on the shape of the starting distribution.

Scaling laws   
Given the result (7) we return to the mean and standard deviation of $`p_k(T)`$ using eq. (4). Substituting $`k=\mathrm{log}_2(L/L_c)`$ one finds the large-$`L`$ scaling laws

$`{\displaystyle \frac{T}{L}}`$ $`=`$ $`{\displaystyle \frac{T_0}{L_c}}{\displaystyle \frac{\beta }{\alpha 1}}{\displaystyle \frac{\sigma _0}{L_c}}\left[BA\left({\displaystyle \frac{L_c}{L}}\right)^p\right],`$ (8)   
$`{\displaystyle \frac{\sigma }{L}}`$ $`=`$ $`{\displaystyle \frac{\sigma _0}{L_c}}A\left({\displaystyle \frac{L_c}{L}}\right)^p.`$ (9)   
The mean travel time $`T`$ and standard deviation $`\sigma `$ scale with $`L`$ with an exponent $`p=\mathrm{log}_2\alpha =0.702`$. This scaling exponent has been studied intensively for the directed polymer problem.

For the seismic problem the primary interest is not the scaling with $`L`$, but the scaling with the strength $`\epsilon `$ of the fluctuations. Perturbation theory gives the $`\epsilon `$-dependence at length $`L_c`$,

$$1v_0T_0/L_c\epsilon ^2L_c/a,v_0\sigma _0/L_c\epsilon \sqrt{a/L_c},$$
(10)   
where $``$ indicates that coefficients of order unity have been omitted. (We will fill these in later.) Since $`L_ca\epsilon ^{2/3}`$ (as mentioned in the introduction), we find upon substitution into eq. (8) the scaling of the mean velocity shift at length $`LL_c`$:

$$\delta v/v_01v_0T/L\epsilon ^{4/3}\left[1+𝒪(L_c/L)^p\right].$$
(11)   
The mean velocity shift saturates at a value of order $`v_0\epsilon ^{4/3}`$. The exponent $`\frac{4}{3}`$ was anticipated in ref.  and is close to the value $`1.33\pm 0.01`$ resulting from simulations.

Comparison with simulations   
For a more quantitative description we need to know the coefficients omitted in eq. (10). These are model dependent. To make contact with the simulations we consider the case of an incoming plane wave instead of a point source. The perturbation theory for the mean velocity shift at length $`L_c`$ gives

$$\delta v_0=v_0\epsilon ^2\frac{L_c}{a}\frac{\sqrt{\pi }}{2}\left(1\frac{2}{\sqrt{\pi }}\frac{a}{L_c}\right).$$
(12)   
The variance at length $`L_c`$ is

$$\delta v^2_0=v_0^2\epsilon ^2\frac{a}{L_c}\sqrt{\pi }\left(1\sqrt{\pi }\epsilon ^2\frac{L_c^3}{a^3}\right).$$
(13)   
We quantify the criterion for the breakdown of perturbation theory by $`L_c=\kappa a\epsilon ^{2/3}`$, with $`\kappa =0.765`$ our single fit parameter. For the non-universal constants $`A`$ and $`B`$ we can use in good approximation $`A=1`$, $`B=1`$. The mean velocity shift in the non-perturbative regime ($`L>L_c`$) is then expressed as:

$$\delta v=\frac{\beta }{\alpha 1}\sqrt{\delta v^2_0}\left[1\left(\frac{L_c}{L}\right)^p\right]+\delta v_0.$$
(14)   
For $`L<L_c`$ we use the perturbative result (12) (with $`L_c`$ replaced by $`L`$). As shown in Fig. 3, the agreement with the computer simulations is quite satisfactory, in particular in view of the fact that there is a single fit parameter $`\kappa `$ for all curves.

Conclusion   
We have presented a non-perturbative theory of the scale-dependent seismic velocity in heterogeneous media. The saturation of the velocity shift at large length scales, observed in computer simulations, is well described by the hierarchical model — including the $`\epsilon ^{4/3}`$ scaling of the saturation velocity. We have concentrated on the case of two-dimensional propagation (for comparison with the simulations), but the $`\epsilon ^{4/3}`$ scaling holds in three dimensions as well. (The coefficients $`\alpha =1.74`$, $`\beta =1.30`$ are different in 3D.)

Our solution of the seismic problem relies on the mapping onto the problem of directed polymers. This mapping holds in the short-wavelength limit, $`\lambda a^2/L`$. To observe the saturation at $`LL_c`$ thus requires $`\lambda a\epsilon ^{2/3}`$. For $`La^2/\lambda `$ the velocity shift will decrease because the velocity fluctuations are averaged out over a Fresnel zone. There exists a perturbation theory for the velocity shift that includes the effects of a finite wavelength. It is a challenging problem to see if these effects can be included into the non-perturbative hierarchical model as well.

Acknowledgements. We are indebted to R. Snieder for suggesting this problem to us, and to X. Leyronas for valuable discussions throughout this work. We acknowledge the support of the Dutch Science Foundation NWO/FOM.