Power spectra of random heterogeneities in the solid earth

Recent seismological observations focusing on the collapse of an impulsive wavelet revealed the existence of small-scale random heterogeneities in the earth medium. The radiative transfer theory (RTT) is often used for the study of the propagation and scattering of wavelet intensities, the mean square amplitude envelopes through random media. For the statistical characterization of the power spectral density function (PSDF) of the random fractional fluctuation of velocity inhomogeneities in a 3-D space, we use an isotropic von Kármán-type function characterized by three parameters: the root mean square (RMS) fractional velocity fluctuation, the characteristic length, and the order of the modified Bessel function of the second kind, which leads to the power-law decay of the PSDF at wavenumbers higher than the corner. We compile reported statistical parameters of the lithosphere and the mantle based on various types of measurements for a wide range of wavenumbers: photo-scan data of rock samples; acoustic well-log data; and envelope analyses of cross-hole experiment seismograms, regional seismograms, and teleseismic waves based on the RTT. Reported exponents of wavenumber are distributed between −3 and −4, where many of them are close to −3. Reported RMS fractional fluctuations are on the order of 0.01–0.1 in the crust and the upper mantle. Reported characteristic lengths distribute very widely; however, each one seems to be restricted by the dimension of the measurement system or the sample length. In order to grasp the spectral characteristics, eliminating strong heterogeneity data and the lower mantle data, we have plotted all the reported PSDFs of the crust and the upper mantle against wavenumber for a wide range (10⁻³–10⁸ km⁻¹). We find that the spectral envelope of those PSDFs is well approximated by the inverse cube of wavenumber. It suggests that the earth-medium randomness has a broad spectrum. In theory, we need to re-examine the applicable range of the Born approximation in the RTT when the wavenumber of a wavelet is much higher than the corner. In observation, we will have to carefully measure the PSDF on both sides of the corner. We may consider the obtained power-law decay spectral envelope as a reference for studying the regional differences. It is interesting to study what kinds of geophysical processes created the observed power-law spectral envelope at different scales and in different geological environments in the solid earth medium.

1 Introduction

The first image of the solid earth is composed of spherical shells, for example, PREM (preliminary reference Earth model) (Dziewonski and Anderson, 1981). As seismic networks were deployed on the regional scale and worldwide, the velocity tomography based on the ray tracing method revealed 3-D heterogeneous structures at various scales; however, spatial variations in the resultant velocity structure are essentially smooth compared with seismic wavelengths. Aki and Chouet (1975) first put a focus on long-lasting coda waves of small earthquakes and interpreted them as scattered waves by small-scale random heterogeneities. They proposed to measure the scattering coefficient g, the scattering power per unit volume as a measure of medium heterogeneity. They analyzed the mean square (MS) amplitude time trace of coda waves as an incoherent sum of scattered wave power by using the Born approximation (e.g., Chernov, 1960), which is a simplified version of the radiative transfer theory (RTT). There have been many measurements of the total scattering coefficient g_iso supposing isotropic scattering (e.g., Sato, 1977a) in various seismotectonic environments. The total scattering coefficient of S waves is reported to be on the order of 10⁻² km⁻¹ for 1–20 Hz in the lithosphere, and it marks a higher value beneath active volcanoes (e.g., Sato et al., 2012; Yoshimoto and Jin, 2008). There were precise measurements of regional variations in g_iso as Carcolé and Sato (2010) in Japan and Eulenfeld and Wegler (2017) in the US. Hock et al. (2004) analyzed medium heterogeneity in Europe from the analyses of teleseismic waves using the modified energy flux model (Korn, 1993). There were also measurements of the anisotropic scattering coefficient from the analysis of S coda envelopes Jing et al. (e.g., 2014); Zeng (e.g., 2017).

Aki and Chouet (1975) derived the angular dependence of the scattering coefficient of scalar waves from the power spectral density function (PSDF) of the fractional velocity fluctuation using the Born approximation. Sato (1984) extended the envelope synthesis of scalar waves to the whole envelope synthesis of three-component seismograms from the P onset to S coda on the bases of the single scattering approximation of the RTT. His syntheses explain how seismogram envelopes in different back azimuths vary depending on the source fault mechanism. Extension to the multiple scattering case was developed by using Stokes parameters (e.g., Margerin et al., 2000; Margerin, 2005; Przybilla et al., 2009; Sanborn et al., 2017). We also note that Monte Carlo simulation was developed to stochastically solve the RTT (e.g., Hoshiba et al., 1991; Gusev and Abubakirov, 1987; Yoshimoto, 2000). For the data processing, it is more appropriate to stack MS envelopes of observed seismograms for comparison with the averaged intensity time traces stochastically synthesized by the RTT (e.g., Shearer and Earle, 2004; Rost et al., 2006; da Silva et al., 2018).

When the central wavenumber of a wavelet increases much larger than the corner wavenumber of the PSDF, the wavelet around the peak value is mostly composed of narrow-angle scattering around the forward direction. In such a case, the Born approximation becomes inappropriate; however, the phase shift modulation based on the parabolic approximation is useful, which is called the phase screen approximation. As an extension of the RTT with the phase screen approximation, the Markov approximation was also used for the analysis of envelope broadening and peak delay with increasing travel distance (e.g., Sato, 1989; Saito et al., 2002; Takahashi et al., 2009). Kubanza et al. (2007) measured regional differences in the lithospheric heterogeneity from the partitioning of seismic energy of teleseismic P waves into the vertical and transverse components based on the Markov approximation.

There have been various kinds of measurements of the PSDF of the random velocity fluctuation, where the PSDF is often supposed to be a von Kármán type. In the following section, the main objective is to compile reported PSDF measurements in various scales in different geological environments of the solid earth: photo scanning of small rock samples, acoustic well logs, array analyses of teleseismic waves, waveform analyses using finite difference (FD) simulations, and analyses of seismogram envelopes on the basis of the RTT. We enumerate their statistical parameters and plot their PSDFs against wavenumber. We will show that the envelope of all the PSDFs is well approximated by a power-law decay curve. Then, we will discuss the results obtained and a few problems in the envelope synthesis theory for such random media and the geophysical origin of such power spectra.

2 Statistical characterization of random media

We consider the propagation of scalar waves as a simple model, where the inhomogeneous velocity is given by $V\left(\mathit{x}\right)=V_{\mathrm{0}}(\mathrm{1}+\mathit{\xi }(\mathit{x}\left)\right)$. The fractional fluctuation ξ(x) is supposed to be a random function of space. We imagine an ensemble of random media {ξ(x)}, where 〈ξ(x)〉=0. Angular brackets mean the ensemble average. We suppose that random media are homogeneous and isotropic, then we statistically characterize them by using the autocorrelation function (ACF):

$$\begin{array}{}\text{(1a)}& {\displaystyle }R\left(\mathit{x}\right)=R\left(r\right)=〈\mathit{\xi }\left(\mathit{y}\right)\mathit{\xi }(\mathit{y}+\mathit{x})〉,\end{array}$$

where $r=\left|\mathit{x}\right|$. The MS fractional fluctuation as a measure of the strength of randomness is supposed to be small, and ${\mathit{\epsilon }}^{\mathrm{2}}\equiv R\left(\mathrm{0}\right)\ll \mathrm{1}$. The Fourier transform of ACF gives the PSDF:

$$\begin{array}{}\text{(1b)}& {\displaystyle }P\left(\mathit{m}\right)=P\left(m\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\iiint }}R\left(\mathit{x}\right)e^{-i\mathit{m}\mathit{x}}\mathrm{d}\mathit{x},\end{array}$$

where wavenumber $m=\left|\mathit{m}\right|$. In some literature, (2π)⁻³ is used as a prefactor in the right-hand side of Eq. (2b).

Figure 1 (a) Log–log plot of PSDF vs. wavenumber m in 3-D space (von Kármán type, κ=0.1, 0.5, and 1; Henyey–Greenstein type, HG, κ=0; Gaussian type, G). (b) Linear plot of ACF vs. lag distance r.

2.1 Several types of random media

There are several types of PSDF and ACF characterized by a few parameters.

2.1.1 von Kármán type

The ACF is written by using a modified Bessel function of the second kind of order κ and characteristic length a:

$$\begin{array}{}\text{(2a)}& {\displaystyle }R\left(r\right)={\displaystyle \frac{{\mathrm{2}}^{\mathrm{1}-\mathit{\kappa }}}{\mathrm{\Gamma }\left(\mathit{\kappa }\right)}}{\mathit{\epsilon }}^{\mathrm{2}}{\left({\displaystyle \frac{r}{a}}\right)}^{\mathit{\kappa }}{K}_{\mathit{\kappa }}\left({\displaystyle \frac{r}{a}}\right)\mathrm{for}\mathit{\kappa }>\mathrm{0},\end{array}$$

which is an exponential type $R\left(r\right)={\mathit{\epsilon }}^{\mathrm{2}}{e}^{-r/a}$ when $\mathit{\kappa }=\mathrm{1}/\mathrm{2}$. In the case of space dimension d, the PSDF is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }P\left(m\right)={\displaystyle \frac{{\mathrm{2}}^d{\mathit{\pi }}^{\frac{d}{\mathrm{2}}}\mathrm{\Gamma }(\mathit{\kappa }+\frac{d}{\mathrm{2}}){\mathit{\epsilon }}^{\mathrm{2}}{a}^d}{\mathrm{\Gamma }\left(\mathit{\kappa }\right){\left(\mathrm{1}+a^{\mathrm{2}}{m}^{\mathrm{2}}\right)}^{\mathit{\kappa }+\frac{d}{\mathrm{2}}}}}\mathrm{for}\mathit{\kappa }>\mathrm{0}\\ \text{(2b)}& {\displaystyle }& {\displaystyle }\propto m^{-\mathrm{2}\mathit{\kappa }-d}\mathrm{for}m\gg a^{-\mathrm{1}}.\end{array}$$

The PSDF obeys a power-law decay at large wavenumbers: $P\left(m\right)\propto m^{-\mathrm{2}\mathit{\kappa }-\mathrm{3}}$ for the 3-D case and $P\left(m\right)\propto m^{-\mathrm{2}\mathit{\kappa }-\mathrm{1}}$ for the 1-D case, where κ corresponds to the Hurst number. In the following, we will basically use a von Kármán-type function for characterizing the earth-medium heterogeneity.

Especially for an anisotropic case, we define the von Kármán-type PSDF in 3-D (e.g., Wu et al., 1994; Nakata and Beroza, 2015):

$$\begin{array}{}\text{(3)}& {\displaystyle }P\left(\mathit{m}\right)={\displaystyle \frac{{\mathrm{2}}^{\mathrm{3}}{\mathit{\pi }}^{\frac{\mathrm{3}}{\mathrm{2}}}\mathrm{\Gamma }(\mathit{\kappa }+\frac{\mathrm{3}}{\mathrm{2}}){\mathit{\epsilon }}^{\mathrm{2}}{a}_x{a}_y{a}_z}{\mathrm{\Gamma }\left(\mathit{\kappa }\right){\left(\mathrm{1}+a_x^{\mathrm{2}}{m}_x^{\mathrm{2}}+a_y^{\mathrm{2}}{m}_y^{\mathrm{2}}+a_z^{\mathrm{2}}{m}_z^{\mathrm{2}}\right)}^{\mathit{\kappa }+\frac{\mathrm{3}}{\mathrm{2}}}}}\mathrm{for}\mathit{\kappa }>\mathrm{0}.\end{array}$$

2.1.2 Henyey–Greenstein type

For a case formally corresponding to κ=0 of the von Kármán-type PSDF, we define the Henyey–Greenstein type ACF and PSDF in 3-D (Henyey and Greenstein, 1941):

$$\begin{array}{}\text{(4a)}& {\displaystyle }& {\displaystyle }R\left(r\right)={\mathit{\epsilon }}^{\mathrm{2}}{K}_{\mathrm{0}}\left({\displaystyle \frac{r}{a}}\right),\text{(4b)}& {\displaystyle }& {\displaystyle }P\left(m\right)={\displaystyle \frac{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}{\mathit{\epsilon }}^{\mathrm{2}}{a}^{\mathrm{3}}}{(\mathrm{1}+a^{\mathrm{2}}{m}^{\mathrm{2}}{)}^{\mathrm{3}/\mathrm{2}}}}\approx \mathrm{2}{\mathit{\pi }}^{\mathrm{2}}{\mathit{\epsilon }}^{\mathrm{2}}{m}^{-\mathrm{3}}\mathrm{for}m\gg a^{-\mathrm{1}}.\end{array}$$

Note that parameter ε² characterizes P but ε²≠R(0) since R(r) diverges as r→0.

2.1.3 Gaussian type

Gaussian-type ACF and PSDF are also used because they are mathematically tractable.

$$\begin{array}{}\text{(5a)}& {\displaystyle }& {\displaystyle }R\left(r\right)={\mathit{\epsilon }}^{\mathrm{2}}{e}^{-\frac{r^{\mathrm{2}}}{a^{\mathrm{2}}}},\text{(5b)}& {\displaystyle }& {\displaystyle }P\left(m\right)=\sqrt{{\mathit{\pi }}^{\mathrm{3}}}{\mathit{\epsilon }}^{\mathrm{2}}{a}^{\mathrm{3}}{e}^{-\frac{m^{\mathrm{2}}{a}^{\mathrm{2}}}{\mathrm{4}}}.\end{array}$$

We plot those PSDFs against wavenumber and ACFs against lag distance in Fig. 1.

3 Measurements of random heterogeneities

There are several kinds of measurements for estimating statistical parameters characterizing random media. Here we principally collect measurements supposing a von Kármán-type function for isotropic randomness; however, we include a few measurements supposing anisotropic randomness and a Gaussian-type function. On a small scale, the photo-scan method is applied to small rock samples. Acoustic well logs are available in deep wells drilled in the shallow crust. When the precise velocity tomography result is available, we can directly calculate the PSDF. In seismology, the most conventional method is to analyze seismograms of natural earthquakes or artificial explosions after traveling through the earth heterogeneity. It is better to focus on MS amplitude envelopes (intensity time traces) since phases are complex and caused by random heterogeneities. Comparing observed seismogram envelopes with envelopes synthesized in random media, we can evaluate von Kármán parameters. For the synthesis, we can use the FD simulations, the RTT with the Born approximation, and the RTT with the phase screen approximation that is equivalent to the Markov approximation. For each reported measurement, we enumerate the target region, data and the method, the measured PSDF as a function of wavenumber m, von Kármán parameters (κ, ε, a), the frequency range, the wavenumber range, and the reference in Tables 1–3. Note that measurements of heterogeneity listed in the Tables are by no means the only ones. Especially in seismological measurements, we estimate the wave number range from the frequency range by using the typical velocity of the target medium. In the Tables, the parameter values in parentheses (…) are a priori fixed in the measurement. Then, we plot obtained PSDFs against wavenumber in Figs. 2–5. When the estimated parameter value is given by a range, a value in square brackets […] is used to represent for plotting PSDFs in the figures. Measurement with a label with an asterisk ^* is insufficient for plotting the PSDF in the figures.

Sivaji et al. (2002)Sivaji et al. (2002)Sivaji et al. (2002)Fukushima et al. (2003)Fukushima et al. (2003)Shiomi et al. (1997)Shiomi et al. (1997)Wu et al. (1994)Holliger (1996)Holliger (1996)

Table 1 Reported von Kaŕmán parameters of rock samples and acoustic well logs. A value in parentheses ( ) is a priori assumed in the measurement. A label with an asterisk ^* is insufficient for plotting the PSDF.

3.1 Photo scan of the rock surface

The photo-scan method uses a scanner to take a picture of the polished flat surface of a small rock sample (e.g., Sivaji et al., 2002; Spetzler et al., 2002; Fukushima et al., 2003). For the case of a granite sample, those papers classified color images on a straight line into three types of mineral grains: quartz, plagioclase, and biotite. Assigning a typical velocity V_P or V_S to each mineral grain, they constructed a velocity profile along the line. Then, they estimated the 1-D PSDF of the velocity fractional fluctuation. They measured 1-D PSDFs of granite and gabbro samples fixing κ=0.5 as R1–R5. Figure 2a shows estimated 1-D PSDFs, where the wavenumber range is on the order of 1 mm⁻¹. We note that raw 1-D PSDFs in Figs. 4 and 5 of Fukushima et al. (2003) decay a little slower than those of R4 and R5 in Fig. 2a, especially at large wavenumbers.

3.1.1 Conversion from 1-D PSDF into 3-D PSDF

In the case of isotropic randomness, we evaluate the 1-D PSDF from the 3-D ACF along the z axis at $x=y=\mathrm{0}$ as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{P}_{\mathrm{1}-D}\left(m_z\right)\equiv \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{R}_{\text{3-D}}\left(\mathrm{0},\mathrm{0},z\right)e^{-i{m}_{z}z}\mathrm{d}z\\ {\displaystyle }& {\displaystyle }=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\left[{\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{2}\mathit{\pi }\right)}^{\mathrm{3}}}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\iiint }}{P}_{\text{3-D}}\left(m_x^{\prime },m_y^{\prime },m_z^{\prime }\right)e^{i{m}_z^{\prime }z}\mathrm{d}{\mathit{m}}^{\prime }\right]e^{-i{m}_{z}z}\mathrm{d}z\\ \text{(6a)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{2}\mathit{\pi }\right)}^{\mathrm{2}}}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\iint }}{P}_{\text{3-D}}\left(m_x^{\prime },m_y^{\prime },m_z\right)\mathrm{d}{m}_x^{\prime }\mathrm{d}{m}_y^{\prime }.\end{array}$$

Substituting Eq. (3a) into the above equation, we have

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{P}_{\mathrm{1}-D}\left(m_z\right)\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{2}\mathit{\pi }\right)}^{\mathrm{2}}}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\iint }}{\displaystyle \frac{\mathrm{8}{\mathit{\pi }}^{\mathrm{3}/\mathrm{2}}{\mathit{\epsilon }}^{\mathrm{2}}{a}^{\mathrm{3}}\mathrm{\Gamma }\left(\mathit{\kappa }+\mathrm{3}/\mathrm{2}\right)}{\mathrm{\Gamma }\left(\mathit{\kappa }\right){\left[\mathrm{1}+a^{\mathrm{2}}\left({m_x^{\prime }}^{\mathrm{2}}+{m_y^{\prime }}^{\mathrm{2}}+m_z^{\mathrm{2}}\right)\right]}^{\mathit{\kappa }+\mathrm{3}/\mathrm{2}}}}\mathrm{d}{m}_x^{\prime }\mathrm{d}{m}_y^{\prime }\\ \text{(6b)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{2}{\mathit{\pi }}^{\mathrm{1}/\mathrm{2}}\mathrm{\Gamma }\left(\mathit{\kappa }+\mathrm{1}/\mathrm{2}\right){\mathit{\epsilon }}^{\mathrm{2}}a}{\mathrm{\Gamma }\left(\mathit{\kappa }\right){\left(\mathrm{1}+a^{\mathrm{2}}{m}_z^{\mathrm{2}}\right)}^{\mathit{\kappa }+\mathrm{1}/\mathrm{2}}}}.\end{array}$$

Thus, we can evaluate the 3-D PSDF from the 1-D PSDF using parameters ε, κ, and a of 1-D PSDF.

Supposing the randomness is isotropic, we evaluate corresponding 3-D PSDFs of R1–R5 and plot them in Fig. 2b.

Figure 2 (a) 1-D PSDF vs. wavenumber for rock samples and acoustic well logs. (b) Converted 3-D PSDF vs. wavenumber, where the randomness is supposed to be isotropic. See labels in Table 1.

3.2 Acoustic well logs in boreholes

An acoustic well log is obtained from the measurement of the travel time of an ultrasonic pulse along the wall of a borehole. Measurements W1 (volcanic tuff) and W2 (tertiary to pre-tertiary) in Japan clearly show power-law decay with κ=0.225 and 0.045, respectively; however, a corner is not clearly seen in each PSDF. Measurement W4 at the deep well KTB in Germany shows κ=0.10. Measurement W3 in the same well shows that the exponent of wavenumber is −0.97, which formally corresponds to a negative κ. Measurement W5 at Cajon Pass in California shows κ=0.11. All these measurements show very small κ values close to 0. Shiomi et al. (1997) made a list of reported exponents of wavenumber, which shows that most κ values are smaller than 0.25. Measurement of a seems to be restricted by the sample length. We enumerate those measurements in Table 1 and plot their 1-D PSDFs against wavenumber in Fig. 2a. Figure 2b plots the corresponding 3-D PSDFs of W4 and W5.

We note that Wu et al. (1994) measured anisotropy of randomness from the analysis of well logs obtained from two parallel wells at KTB: the ratio of characteristic scales in horizontal to vertical directions $a_h/a_z=\mathrm{1.8}$ (see Eq. 3) as shown in W3.

Nakata and Beroza (2015)da Silva et al. (2018)Rothert and Ritter (2000)Gaebler et al. (2015)Przybilla et al. (2009)Calvet and Margerin (2013)Sens-Schönfelder et al. (2009)Sens-Schönfelder et al. (2009)Emoto et al. (2017)Takemura et al. (2009)Kobayashi et al. (2015)Takahashi et al. (2009)Saito et al. (2005)Sanborn et al. (2017)Sanborn et al. (2017)Wang and Shearer (2017)Wang and Shearer (2017)Petukhin and Gusev (2003)Hock et al. (2004)Hock et al. (2004)Yoshimoto et al. (1997b)Takahashi et al. (2009)Morioka et al. (2017)Sens-Schönfelder et al. (2009)Furumura and Kennett (2005)Aki (1973)Capon (1974)Powell and Meltzer (1984)Flatté and Wu (1988)

Table 2 Reported von Kaŕmán parameters of the lithosphere including the crust and the uppermost mantle. A value in parentheses ( ) is a priori assumed in the measurement. When the estimated value is given by a range, a value in square brackets [ ] is used for plotting the PSDF. A label with an asterisk ^* is insufficient for plotting the PSDF.

3.3 Velocity tomography

There have been measurements of velocity tomography at various scales, from which we can calculate the PSDF and then estimate von Kármán-type parameters. This method depends on the spatial resolution of the tomography result. Measurement L1 in Table 2 is calculated from the precise V_P tomography result of the shallow crust, Los Angeles, California: the exponent of wavenumber is −3.08 (κ=0.04). Anisotropic randomness is also reported: a_z=0.1 and a_h=0.5 km (see Eq. 3), we show those in Fig. 3a. Measurement M2 in Table 3 is evaluated from the 2-D PSDF of the V_S tomography result of the upper mantle in a low wavenumber range. Although there is a resolution limit of the tomography method, the exponent of wavenumber is between −2 and −3, which means $\mathrm{0}<\mathit{\kappa }<\mathrm{0.5}$. We note that Fig. 8 of Mancinelli et al. (2016a) shows that the 1-D PSDF estimated from the V_P tomography result in the upper mantle (Meschede and Romanowicz, 2015) covers that of MU2 (κ=0.05, ε=0.1, a=2000 km) for the wavenumber range from $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{4}}$ to 10⁻² km⁻¹.

Chevrot et al. (1998)Mancinelli et al. (2016a)Shearer and Earle (2004)Shearer and Earle (2004)Cormier (1999)Mancinelli et al. (2016b)Zhang et al. (2018)Margerin and Nolet (2003)Bentham et al. (2017)

Table 3 Reported von Kaŕmán parameters of the upper mantle and the lower mantle. A value in parentheses ( ) is a priori assumed in the measurement. When the estimated value is given by a range, a value in square brackets [ ] is used for plotting the PSDF. A label with an asterisk ^* is insufficient for plotting the PSDF.

Figure 3 3-D PSDF vs. wavenumber for (a) the lithosphere (the crust and uppermost mantle), (b) strong heterogeneities and array data analyses in the lithosphere. See labels in Table 2.

3.4 Array analysis of teleseismic P waves

Teleseismic P waves registered by a large aperture array were used for the evaluation of the 3-D PSDF of the lithosphere beneath the array: LA1 and LA2 of Table 2 in Montana and LA3 in southern California used amplitude and phase coherence analyses, where a Gaussian-type PSDF (Eq. 5b) was assumed because of mathematical simplicity. As shown in Fig. 3b, they drop very fast as wavenumber increases. Later Flatté and Wu (1988) developed the angular coherence analysis in addition to the above methods. Analyzing teleseismic P waves registered at NORSAR, they proposed an overlapping two-layer model LA4, which is composed of a band-limited flat spectrum from the surface to 200 km of depth and m⁻⁴ spectrum (κ=0.5, ε=0.01–0.04) for depths from 15 to 250 km. It means κ<0.5 and the decay of their PSDF is much smaller than that of Gaussian types (not shown in Fig. 3b).

Figure 4 3-D PSDF vs. wavenumber for the upper and lower mantle. See labels in Table 3.

3.5 Finite difference (FD) simulations

FD simulation is often used for the numerical simulation of waves in an inhomogeneous velocity structure. For the evaluation of average MS amplitude envelopes, we have to repeat simulations of the wave propagation through random media having the same PSDFs that are generated by using different random seeds. There are several measurements of statistical parameters using FD such as L9–L11 and ML4 in Tables 2 and 3. Measurement LS5 focused on the fact that the subducting oceanic plate is an efficient waveguide for high-frequency seismic waves: estimated anisotropic parameters are κ=0.5 and ε=0.02 with a_p=10 km and a_t=0.5 km in the parallel and transverse directions, respectively. Note that ML2 supposes a Gaussian-type function.

Figure 5 3-D PSDF vs. wavenumber for all the data.

3.6 Analyses of seismogram intensities (MS amplitude envelopes)

The RTT is essentially stochastic to directly synthesize the intensity (the average MS amplitude envelope) of a wavelet propagating through random media. There are two conventional methods on the basis of the RTT: one uses the Born approximation and the other uses the phase screen approximation based on the parabolic approximation when the wavenumber is larger than the corner. The former neglects the phase shift, but the latter correctly considers the phase shift.

3.6.1 Scalar wave scattering by a single obstacle

We here study the deterministic scattering of scalar waves by a single spherical obstacle (radius a=5 km and velocity anomaly $\mathit{\epsilon }=+\mathrm{0.05}$) embedded in a homogeneous medium (V₀=4 km s⁻¹) as a mathematical model. The Born approximation calculates spherically outgoing scattered waves putting the incident plane wave of wavenumber k_c in the interaction term of the first-order perturbation equation. From the scattering amplitude we evaluate the total scattering cross section σ₀ as a measure of scattering power of the obstacle. The resultant σ₀ monotonously increases with frequency as shown by a blue line in Fig. 6. As the wavenumber increases (ak_c≫1), the phase shift increases as the incidence plane wave penetrates into the obstacle. Putting the phase modulated wave according to the parabolic approximation (the phase screen approximation) into the interaction term of the first-order perturbation equation, we calculate the scattering amplitude and then the total scattering cross section. This is the distorted-wave Born approximation with the phase screen approximation, which is also referred to as the eikonal approximation. This approximation predicts that σ₀ (a red line in Fig. 6) saturates at high frequencies and converges to 2πa², which is twice the geometrical section area of the obstacle as predicted by shadow scattering (e.g., Landau and Lifshitz, 2003, p. 519 and 543). We recognize that the conventional Born approximation is still accurate even for ak_c>1; however, it works well only for ${\mathit{\epsilon }}^{\mathrm{2}}{a}^{\mathrm{2}}{k}_{\mathrm{c}}^{\mathrm{2}}\mathit{\lesssim }O\left(\mathrm{0.1}\right)$. We should use the distorted-wave Born approximation with the phase screen approximation for ${\mathit{\epsilon }}^{\mathrm{2}}{a}^{\mathrm{2}}{k}_{\mathrm{c}}^{\mathrm{2}}\mathit{\gtrsim }O\left(\mathrm{1}\right)$. The two approximations predict nearly the same σ₀ value in the intermediate range. We note that 2εak_c is the phase shift on the center line after passing the obstacle. Note that the phase screen approximation is not applicable for ak_c<1 since it is based on the parabolic approximation.

Interpreting ε and a as the RMS fractional fluctuation and the characteristic length of uniformly distributed random media, we may use the inequality ${\mathit{\epsilon }}^{\mathrm{2}}{a}^{\mathrm{2}}{k}_{\mathrm{c}}^{\mathrm{2}}\ll O\left(\mathrm{1}\right)$ or ${\mathit{\epsilon }}^{\mathrm{2}}{a}^{\mathrm{2}}{k}_{\mathrm{c}}^{\mathrm{2}}\mathit{\lesssim }O\left(\mathrm{0.1}\right)$ as a criterion of the Born approximation used in the RTT.

Figure 6 Deterministic scattering of scalar waves by a high-velocity sphere. (a) Log–log plot of the total scattering cross section against frequency. (b) Semi-log plot for zoomed in graph. The Born approximation and the distorted-wave Born approximation with the phase screen approximation are drawn by blue and red lines, respectively.

3.6.2 RTT with the Born approximation

For uniformly distributed random media characterized by P(m), the Born approximation leads to the scattering coefficient at wavenumber k_c into scattering angle ψ:

$$\begin{array}{}\text{(7a)}& {\displaystyle }g(k_{\mathrm{c}},\mathit{\psi })={\displaystyle \frac{k_{\mathrm{c}}^{\mathrm{4}}}{\mathit{\pi }}}P\left(\mathrm{2}{k}_{\mathrm{c}}\sin {\displaystyle \frac{\mathit{\psi }}{\mathrm{2}}}\right),\end{array}$$

which is axially symmetric. The total scattering coefficient is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{g}_{\mathrm{0}}\left(k_{\mathrm{c}}\right)\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\oint g(k_{\mathrm{c}},\mathit{\psi })\mathrm{d}\mathrm{\Omega }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}g(k_{\mathrm{c}},\mathit{\psi })\sin \mathit{\psi }\mathrm{d}\mathit{\psi }\\ \text{(7b)}& {\displaystyle }& {\displaystyle }=\underset{\mathrm{0}}{\overset{\mathrm{2}{k}_{\mathrm{c}}}{\int }}{g}_{\ker }(k_{\mathrm{c}},m)\mathrm{d}m,\end{array}$$

where $m=\mathrm{2}{k}_{\mathrm{c}}\sin \frac{\mathit{\psi }}{\mathrm{2}}$. The integral kernel in the wavenumber space is given by

$$\begin{array}{}\text{(7c)}& {\displaystyle }{g}_{\ker }(k_{\mathrm{c}},m){\displaystyle }={\displaystyle \frac{k_{\mathrm{c}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi }}}mP\left(m\right).\end{array}$$

The upper bound of the integral is twice the wavenumber. As an example, Fig. 7 shows plots of P(m) (blue) vs. m, and g_ker(m) vs. m at 0.1 Hz (red) and 1 Hz (green) for the case of κ=0.5, ε=0.05, a=1 km, and V₀=4 km s⁻¹. As shown at the upper-right corner, the scattering pattern at 1 Hz (green) has a large lobe into the forward direction; however, it becomes isotropic as the frequency decreases. Dots on each g_ker curve show corresponding scattering angles.

Figure 7 Plot of P(m) (blue, right scale) and the spectral kernel of the scattering coefficient g_ker(k_c(f_c),m) at f_c=0.1 Hz (red, left scale) and 1 Hz (green, left scale) according to the Born approximation. Scattering angles are marked by dots on each trace. For the case of the frequency band between 0.1 and 1 Hz, the phase screen approximation based on the parabolic approximation covers the wavenumber range from 0 to the upper bound (line at the bottom); however, the Born approximation covers the range from 0 to twice the upper bound (line next to the bottom). We also use these line styles in Figs. 3–5 and 9.

In the framework of the RTT, the Monte Carlo simulation is a simple method to stochastically synthesize the wavelet intensity time trace. A particle carrying unit intensity is shot randomly from a point source, and its trajectory is traced with the increment of time steps. The occurrence of scattering is stochastically tested by the inequality g₀V₀Δt> Random[0,1] at every time step of Δt, and $g(k_{\mathrm{c}},\mathit{\psi })/\left(\mathrm{4}\mathit{\pi }{g}_{\mathrm{0}}\right(k_{\mathrm{c}}\left)\right)$ is used as the probability of scattering into angle ψ. Note that Random[0,1] is a uniform random variable between 0 and 1. Since g₀V₀Δt is chosen to be small enough, scattering does not occur every time step but intermittently. As a simple example, Fig. 8a schematically illustrates the flowchart of the Monte Carlo simulation for the isotropic radiation from a point source in uniform random media. At lapse time t, dividing the number of particles n registered in a spherical shell of radius r and a thickness Δr by the total number of particles N and the shell volume 4πr²Δr, we calculate the intensity Green function G(r,t). The intensity time trace I(r,t) is calculated by the convolution of G(r,t) and the source intensity time function S(t) in the time domain. It is easy to introduce a layered structure of background velocity and intrinsic absorption into the simulation code.

Figure 8 Flowchart of the Monte Carlo simulation code according to the RTT for the scalar wavelet intensity in uniform random media. (a) RTT with the Born approximation. (b) RTT with the phase screen approximation.

The RTT for the scalar wave case can be extended to the elastic vector wave case by using Stokes parameters. There are four scattering modes: PP, PS, SS, and SP scatterings, and the S-wave scattering coefficients are not axially symmetric (see Sato et al., 2012, Fig. 4.7). Many papers (e.g., Shearer and Earle, 2004; Przybilla et al., 2009) suppose proportional relations $\mathit{\delta }{V}_{\mathrm{p}}/V_{\mathrm{P}\mathrm{0}}=\mathit{\delta }{V}_{\mathrm{S}}/V_{\mathrm{S}\mathrm{0}}=\mathit{\xi }$ and $\mathit{\delta }\mathit{\rho }/{\mathit{\rho }}_{\mathrm{0}}=\mathit{\nu }\mathit{\xi }$ based on the empirical Birch's law, which reduce three fractional fluctuations into one (e.g., Sato et al., 2012, Eqs. 4.58 and 4.59).

The RTT with the Born approximation has been often used not only for the analyses of S coda envelopes but also for the whole seismogram envelope from the P onset via P coda through S wave until S coda (see Tables 2 and 3). This method has been often used not only for the analyses of regional seismograms propagating through the lithosphere but also for the analyses of teleseismic waves propagating through the mantle. This method is not only applied to direct P phase waves but also PcP and PKP_prec phase waves and so on. In this review, we neglect intrinsic attenuation parameters a priori assumed or measured in each paper. For a given wavenumber range (k_l,k_u) (gray) in Fig. 7, each PSDF curve using this method in Figs. 3–5 and 9 is drawn by a dotted line for (0,k_l) and a solid line for (k_l,2k_u) as the line next to the bottom of Fig. 7. As indicated by dots on the g_ker curves, the wavenumber interval of the solid line reflects wide-angle scattering and that of dotted line reflects narrow-angle scattering around the forward direction.

Most measurements of S waves in the lithosphere cover the wavenumber range up to 100 km⁻¹. Measurement L2 analyzed cross-hole seismograms on the order of kHz by using 2-D RTT, of which the wavenumber range is as high as on the order of 1 m⁻¹. Measurement MU2 for teleseismic P-wave envelopes at long periods in the upper mantle shows that the characteristic scale a=2000 km is much larger than those of MU3 and ML1 at short periods. Several measurements a priori suppose κ=0.5; however, most measurements show κ<0.5 except L6. Measurements ML3 and MW1 propose the HG type (see Eq. 4b) corresponding to κ=0 for the lower or whole mantle. We note that Mancinelli et al. (2016b) proposed an alternative model of 3-D PSDF $\propto m^{-\mathrm{2.6}}$ in addition to ML3 (not shown in Fig. 4).

3.6.3 RTT with the phase screen approximation

When ak_c≫1, scattering mostly occurs within a narrow angle around the forward direction. At a large travel distance, the wavelet just after the onset is mostly composed of those waves. The phase screen approximation correctly calculates the phase shift modulation. For the incidence of a plane wave into the z direction, the mutual coherence function (MCF) of the phase shift modulated waves for an increment Δz is given by

$$\begin{array}{}\text{(8a)}& {\displaystyle }{\displaystyle }\mathrm{\Phi }(k_{\mathrm{c}},r_{\perp },\mathrm{\Delta }z)=e^{-k_{\mathrm{c}}^{\mathrm{2}}\left(A\right(\mathrm{0})-A(r_{\perp }\left)\right)\mathrm{\Delta }z}.\end{array}$$

The longitudinal integral of the ACF is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }A\left(r_{\perp }\right)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}R({\mathit{x}}_{\perp },z)\mathrm{d}z\\ \text{(8b)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{(\mathrm{2}\mathit{\pi }{)}^{\mathrm{2}}}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\iint }}P({\mathit{m}}_{\perp },m_z=\mathrm{0})e^{i{\mathit{m}}_{\perp }{x}_{\perp }}\mathrm{d}{\mathit{m}}_{\perp },\end{array}$$

where x_⊥ is the transverse coordinate vector and $r_{\perp }=\left|{\mathit{x}}_{\perp }\right|$ Sato et al. (2012, Eq. 9.60). Taking the Fourier transform of MCF Φ with respect to transverse coordinates, we have

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\stackrel{\mathrm{˘}}{\mathrm{\Phi }}(k_{\mathrm{c}},k_{\perp },\mathrm{\Delta }z)={\displaystyle \frac{\mathrm{1}}{(\mathrm{2}\mathit{\pi }{)}^{\mathrm{2}}}}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\iint }}\mathrm{\Phi }(k_{\mathrm{c}},r_{\perp },\mathrm{\Delta }z)e^{i{k}_{\perp }{x}_{\perp }}\mathrm{d}{\mathit{x}}_{\perp }\\ \text{(8c)}& {\displaystyle }& {\displaystyle }\underset{\mathrm{\Delta }z\to \mathrm{0}}{⟶}\mathit{\delta }\left({\mathit{k}}_{\perp }\right).\end{array}$$

Since ${\iint }_{-\mathrm{\infty }}^{\mathrm{\infty }}\stackrel{\mathrm{˘}}{\mathrm{\Phi }}(k_{\mathrm{c}},k_{\perp },\mathrm{\Delta }z)\mathrm{d}{k}_{\perp }=\mathrm{1}$, interpreting $\stackrel{\mathrm{˘}}{\mathrm{\Phi }}(k_{\mathrm{c}},k_{\perp },\mathrm{\Delta }z)$ as the probability of ray-bending angle $\mathit{\psi }={\tan }^{-\mathrm{1}}\frac{k_{\perp }}{k_{\mathrm{c}}}$ per increment Δz=V₀Δt, we can stochastically synthesize the intensity by using the Monte Carlo simulation (e.g., Williamson, 1972; Takahashi et al., 2008; Saito et al., 2008). As a simple example, Fig. 8b schematically illustrates the flowchart of the RTT with the phase screen approximation for the isotropic radiation from a point source in uniform random media. Different from the Born approximation, narrow-angle ray bending occurs at every time step. The intensity Green function can be obtained in the same manner as the RTT with the Born approximation. This approximation synthesizes the intensity time trace having a delayed peak from the onset and a decaying tail of early coda at large travel distances. This approximation can not synthesize the late coda intensity since wide-angle scattering is neglected. The Markov approximation is known as a stochastic extension of the phase screen method for the two-frequency mutual coherence function (e.g., Saito et al., 2002). If we focus on the intensity time trace around the peak arrival, the Markov approximation and the RTT with the phase screen approximation show good coincidence (see Sato and Emoto, 2018, Fig. 8).

When this approximation is used, $k_{\mathrm{c}}\gg a^{-\mathrm{1}}$ is a priori supposed. Most of this type of measurement reads the peak delay and the envelope width of each seismogram envelope. There is some merit to the fact that the peak delay measurement is rather insensitive to intrinsic absorption. In NE Japan, a κ value beneath a volcano LS2 is smaller than those in the fore-arc side L12 and L13. Note that narrow-angle scattering around the forward direction dominates in teleseismic wavelets even if the Born approximation is used for the analysis. Narrow-angle scattering is mostly produced by the PSDF in low wavenumbers compared with k_c. For a given wavenumber range (k_l,k_u) (gray) in Fig. 7, each PSDF curve using this method in Fig. 3 is drawn by a dotted line for (0,k_l) and a solid line for (k_l,k_u) as the bottom line of Fig. 7.

Figure 9 3-D PSDF vs. wavenumber for the crust and the upper mantle. Data of Gaussian-type, anisotropy type, strong heterogeneity, the lower mantle, and the whole mantle are excluded. The light-gray straight line visually fits to most spectral envelopes.

3.7 Characteristics of reported PSDFs

3.7.1 All the data

Some measurements a priori assumed κ=0.5; however, most of measurements report κ<0.5. In the mantle, κ is very small and close to zero, and a HG-type function is also proposed. The RMS fractional fluctuation ε is on the order of 0.1 for rock samples and well-log data, and in the range from 0.01 to 0.1 in the lithosphere and the upper mantle. Large values are reported for the shallow crust L16 and beneath a volcano LS3; however, smaller values are reported for the lower mantle. The characteristic scale a varies a lot depending on measurements. The corner wavenumber a⁻¹ is not clearly seen in PSDFs of acoustic well logs. Some measurements report anisotropy: W3 of well-logs, L1 of velocity tomography in the shallow crust, and LS5 in the subducting oceanic slab. The characteristic length in the vertical direction is smaller than the horizontal direction in the shallow crust, and that in the transverse direction is smaller than that in the direction parallel to the subducting slab.

Plotting PSDFs against wavenumber is more informative for understanding the random heterogeneity compared with enumerating statistical parameter values. Figure 5 shows the plot of 3-D PSDF vs. wavenumber for all the data covering a wide wavenumber range (10⁻⁴–10⁸ km⁻¹). We recognize that the Gaussian type PSDFs show a very different behavior from others, which suggests the Gaussian-type assumption is inappropriate. PSDFs in the lower mantle take smaller values, and those for volcanoes and for the shallow crust take larger values than others.

3.7.2 Lithosphere and the upper mantle except strong heterogeneity, Gaussian, and anisotropy types

Eliminating data supposing a Gaussian-type data LA1–LA3, strong heterogeneity data LS1–LS4, anisotropy-type data L1 and LS5, and the lower mantle and the whole mantle data ML1–ML4 and MW1–MW2 from Fig. 5, we plot the rest of the data for the crust and the upper mantle in Fig. 9. Most ε values are in the range of 0.01–0.1, most κ values are less than or equal to 0.5 and many of them are close to 0, and the high wavenumber end of the power-law decay branch of each PSDF is not far from each corner wavenumber.

We draw a power-law decay line PSDF$\left(m\right)=\mathrm{0.01}{m}^{-\mathrm{3}}$ km³ (gray) visually fitting to most PSDF envelopes for a very wide range of wavenumbers (10⁻³–10⁸ km⁻¹). This line is not the average of PSDFs. This line looks like an extension of MU2 in the upper mantle into higher wavenumbers of the shallow crust.

4 Discussions

4.1 Measurements

It will be necessary for us to measure the small-scale randomness of sedimentary rock samples. More measurements are necessary in the wavenumber range 10³–10⁵ km⁻¹ since there are few measurements.

In most PSDF measurements, each power-law decay branch is short since the Born approximation senses the spectrum up to twice the wavenumber. It will be necessary to measure how each power-law decay branch varies with wavenumber increasing. It will be necessary to estimate the corner a⁻¹ in each measurement with a wide wavenumber range covering, sufficiently large enough, both sides of the corner. The flat part, the low-wavenumber side of each PSDF drawn by a dotted line in figures is also important as the cause of narrow-angle scattering.

Although most measurements used in this review analyzed intrinsic attenuation, we did not enumerate them in this review since different assumptions were used in different measurements. It will be necessary for us to systematically measure the PSDF of random heterogeneity in conjunction with intrinsic attenuation.

We should note that there are large variations in δln V_S∕δln V_P and $\mathit{\nu }\equiv \mathit{\delta }\ln \mathit{\rho }/\mathit{\delta }\ln V_{\mathrm{S}}$ in the earth. Koper et al. (1999) estimated δln V_S∕δln V_P to be in the range 1.1–1.5 in the Tonga Slab. Romanowicz (2001) estimated δln V_S∕δln V_P to be larger than 2.5 in the lower mantle at larger length scales. Many measurements reported use ν=0.8 for the synthesis, which is appropriate for the shallow lithosphere. Parameter ν takes smaller values of 0.17 for volcanic-tuff (Shiomi et al., 1997) and 0.31–0.33 for sandstone and shale (Kenter et al., 2007). In the mantle, Karato (2008) estimated ν=0.23–0.42 for the S-wave velocity predicted from the temperature derivatives of seismic wave velocities and thermal expansion, and ν=0.15–0.23 including the influence of anelasticity. It will be necessary to introduce realistic δln V_S∕δln V_P and δln ρ∕δln V_S in the synthesis.

Figure 9 summarizes reported PSDF measurements supposing isotropic randomness; however, there are measurements clearly showing anisotropic randomness such as W3 and L1 for the shallow crust and LS5 for the oceanic slab. Those may reflect the effect of gravity for the creation of anisotropy. It will be necessary for us to study how a wavelet propagates through anisotropic random heterogeneity of the earth medium (e.g., Margerin, 2006).

4.2 Mathematical theory

In Sect. 3.6.1, we mentioned that the conventional Born approximation is inapplicable and the phase screen approximation is useful when the phase shift becomes large as the wavenumber increases. In order to avoid the difficulty, taking the center wavenumber of a wavelet as a reference, Sato and Emoto (2018) proposed to divided the PSDF into two components (see also Sato, 2016; Sato and Emoto, 2017). They use the Born and phase-screen approximations to the short-scale (high wavenumber) and long-scale (low wavenumber) components, P_S and P_L, respectively, in the RTT in order to simultaneously explain the envelope broadening just after the onset and the excitation of late coda waves. Figure 10 illustrates the flowchart of their Monte Carlo simulation. Their spectrum division method looks like an implementation of the distorted-wave Born approximation in the RTT since it describes wide-angle scattering for the incidence of the phase-shift modulated wave. They successfully synthesized intensity time traces that explain FD simulation results for the case of ak_c=23.6 and ${\mathit{\epsilon }}^{\mathrm{2}}{a}^{\mathrm{2}}{k}_{\mathrm{c}}^{\mathrm{2}}=\mathrm{1.39}$. It would be interesting to see how this method may be extended to polarized elastic waves.

We note that some papers numerically show that the RTT with the Born approximation works well in some cases over the above limitation. Przybilla et al. (2006) synthesizes vector-wave intensity that fits to that of the FD simulation in 2-D, even for S waves of ak_c=58 and ${\mathit{\epsilon }}^{\mathrm{2}}{a}^{\mathrm{2}}{k}_{\mathrm{c}}^{\mathrm{2}}=\mathrm{8.4}$ (see their Table 1) if the wandering effect is convolved as a result of the travel time fluctuation. Emoto and Sato (2018) show that the synthesized scalar intensity by the RTT with the Born approximation fits to that of the FD simulation in 3-D, even for the case of ak_c=23.6 and ${\mathit{\epsilon }}^{\mathrm{2}}{a}^{\mathrm{2}}{k}_{\mathrm{c}}^{\mathrm{2}}=\mathrm{1.39}$ when the wandering effect is convolved. If the earth heterogeneity is represented by a power-law decay power spectrum for such a wide wavenumber range, which means that the corner wavenumber is very low, we should carefully examine the applicability of the Born approximation in the RTT.

Figure 10 Spectrum division method. (a) Division of P into two components, P_S and P_L, with respect to the center wavenumber of a wavelet k_c as a reference, where ζ is a tuning parameter. (b) Flowchart of the Monte Carlo simulation according to the RTT with the spectrum division method. Modified from Sato and Emoto (2018).

Acoustic well-log and photo-scan methods faithfully measure the inhomogeneous elastic coefficients. The RTT used here supposes the scattering contribution of random inhomogeneity of elastic coefficients only; however, observed seismograms do not only reflect those but also the scattering contribution of pores and cracks distributed over the earth medium. It will be necessary for us to study their contribution in the intensity synthesis.

4.3 Power-law decay spectral envelope

In observation, we may take the power-law spectral envelope as a reference curve for studying the regional differences, especially in the power-law decay part of the PSDF. The characteristic length a seems to increase as the wavenumber of a wavelet decreases or as the dimension of measurement system becomes large. It reminds us that the characteristic scale of the slip distribution increases with increasing source dimension as Mai and Beroza (2002) analyzed finite-fault source inversion results (see their Fig. 12).

The power-law decay spectral envelope reminds us of the observed fractal nature of various kinds of surface topographies: Sayles and Thomas (1978a, b) show a 1-D PSDF $\propto m^{-\mathrm{2}}$ for wavelengths 10⁻⁶–10³ km although the power exponent varies from −1.07 to −3.03 for small segments; Brown and Scholz (1985) show a 1-D PSDF $\propto m^{-\mathrm{1.64}\mathrm{to}-\mathrm{3.36}}$ for the wavenumber range 10⁻⁶ to 0.1 µm⁻¹, especially for the topography of natural rock surfaces and faults. We also note that the PSDF of the refractive index fluctuation of air obeys the Kolmogorov spectrum: 3-D PSDF $\propto m^{-\mathrm{11}/\mathrm{3}}$, where $\mathit{\kappa }=\mathrm{1}/\mathrm{3}$. This spectrum is physically produced by the cascade in the turbulent flow of low viscosity air: the large eddies breaks up originating smaller eddies dissipating energy by viscosity. However, it may be difficult to apply this cascade model to the mantle since the viscosity of mantle fluid is thought to be high.

For igneous rocks such as granite, there are variations in composition of minerals and grain sizes, which depend on a variety of slow crystallization differentiations of basaltic magma. Random variations in acoustic well-log profiles reflect the complex sedimentation process during the geological history. Volcanism produces more heterogeneous structures composed of pyroclastic material and lava. For random heterogeneities in the mantle, we imagine complex mantle flow. Mancinelli et al. (2016a) referred to a marble cake mantle model (Allègre and Turcotte, 1986) containing heterogeneity mostly composed of basalt and harzburgite at many scales in the upper mantle in order to explain the power-law spectrum. Stixrude and Lithgow-Bertelloni (2007) proposed the velocity variation due to chemical and phase stability at different depths, which is a possible candidate especially for heterogeneity in the vertical direction. If we accept the power-law decay spectrum, we will have to study what kinds of geophysical mechanisms created such random medium spectra at different scales and in different geological environments in the solid earth.

4.4 Isotropic scattering coefficient

In advance of the measurements based on the RTT for anisotropic scattering presented here, there have been many measurements of the isotropic scattering coefficient g_iso in the world on the basis of the RTT with the isotropic scattering assumption (e.g., Sato et al., 2012; Yoshimoto and Jin, 2008). The isotropic scattering model is mathematically tractable, and the multiple lapse-time window analysis (Fehler et al., 1992; Hoshiba, 1993) has often been used for practical analyses. This method essentially estimates g_iso from the ratio of late coda excitation to the radiated energy irrespective of the envelope broadening. Recent measurements show that g_iso decreases with depth (e.g., Rachman and Chung, 2016; Badi et al., 2009). It will be interesting to plot the frequency dependence of reported g_iso for a wide range of frequencies and to study the relation with the obtained power spectral envelope shown in Fig. 9.

5 Conclusions

Recent seismological observations focusing on the collapse of an impulsive wavelet revealed the existence of small-scale random heterogeneities in the earth medium. The RTT has often been used for the study of the propagation of wavelet intensities, the MS amplitude envelopes. For the statistical characterization of the PSDF of random velocity inhomogeneities in a 3-D space, we have used von Kármán type functions with three parameters: the RMS fractional velocity fluctuation ε, the characteristic length a, and the order κ of the modified Bessel function of the second kind. This model leads to the power-law decay of PSDF $\propto m^{-\mathrm{2}\mathit{\kappa }-\mathrm{3}}$ at wavenumber m higher than the corner at a⁻¹. We have compiled reported statistical parameters of the lithosphere and the mantle based on various types of measurements for a wide range of wavenumbers: photo-scan data of rock samples, acoustic well-log data, and envelope analyses of cross-hole experiment seismograms, regional seismograms, and teleseismic waves based on the RTT. Reported κ values are distributed between 0 and 0.5 (PSDF $\propto m^{-\mathrm{3}\mathrm{to}-\mathrm{4}}$), where many of them are close to 0 (PSDF $\propto m^{-\mathrm{3}}$). Reported ε values are on the order of 0.01–0.1 in the crust and the upper mantle, where smaller values in the lower mantle and higher values beneath volcanoes. Reported a values distribute very widely; however, each one seems to be restricted by the dimension of the measurement system or the sample length. In order to grasp the spectral characteristics, eliminating strong heterogeneity data and the lower mantle data, we have plotted all the reported PSDFs in the crust and the upper mantle against wavenumber m for a wide range (10⁻³–10⁸ km⁻¹). We find that the envelope of those PSDFs is well approximated by a power-law decay curve 0.01 m⁻³ km³. Multiple plots of PSDFs and the power-law decay spectral envelope obtained require us to do the following. In theory, it will be necessary to examine whether the Born approximation is reliable or not if the wavenumber increases much larger than the corner; in observation, we will have to more carefully examine the behavior of each PSDF on both sides of the corner. If we accept the power-law decay spectral envelope, it suggests that the earth-medium randomness has a broad spectrum. We may consider the obtained power-law decay spectral envelope as a reference for studying the regional differences. It is interesting to study what kinds of geophysical processes created the power-law spectral envelope at different scales and in different geological environments of the solid earth medium.