Scale lengths of heterogeneities under Tibet

1.   Introduction

Mapping seismic velocity variations in the interior of the Earth provides a powerful way to interrogate dynamics of the Earth (Dziewonski, 2005). Subducted slabs are cold and usually associated with fast seis-mic velocities; and hot upwellings like the deep mantle plumes and partial melts under the mid ocean ridges are commonly associated with slow velocities. These features can be traced deep into the mantle by seismic tomography (Grand et al., 1997). The spatial reso-lution of global tomography is on the order of thou-sands of kilometers. However many geological units such as the oceanic crust thickness, seamounts, and slabs are on the order of tens of kilometers. In the de-terministic tomography, gross volumetric averages of such structures largely miss the geological complex-ity, giving very reductionist visions of the structure and processes occurring in the dynamic mantle system. Although progresses have been made in terms of new inversion methods and more seismometers deployed, we are still far from producing tomograms with a spa-tial resolution on the scale of these interesting tectonic units. In the traveltime tomography, there is signifi-cant power in the traveltime residuals that are not ac-counted for by the final velocity model (Maceira et al., 2011). The residuals can be explained by random het-erogeneities in the model, which are missed by the de-terministic tomography.

Therefore, it is natural to view the seismic velocity variations in terms of a deterministic component plus a stochastic component. Small-scale statistical veloc-ity variations also carry important information about the Earth's dynamics. The perturbation strengths and length scales of these small-scale heterogeneities provide a way to study convective mixing of the mantle, which may have far-reaching consequences in petrology, min-eralogy and geochemistry (Bercovici, 2007).

The first statistical study of the Earth interior dated back to Aki (1973), which predated the mod-ern deterministic seismic tomography. Following Aki's initial work to characterize heterogeneities under the Large Aperture Seismic Array (LASA), Montana, sev-eral other studies were carried out for Norwegian Seis-mic Array (NORSAR) and southern India (Capon, 1974; Capon and Berteussen, 1974; Berteussen et al., 1977; Flatté and Wu, 1988). For inverting the hetero-geneity spectrum, several types of coherence functions are proposed. Aki (1973) used the transverse coher-ence function (TCF) or the Chernov theory (Chernov, 1960). Flatté and Wu (1988) extended the method to the angular coherence function (ACF). And then Wu and Flattxe (1990) and Chen and Aki (1991) indepen-dently proposed the joint transverse and angular coher-ence function (JTACF). However, these important theo-retical progresses are all for a homogeneous background velocity. Recent work by Zheng and Wu (2008) has ex-tended the theories to having a 1D variable background velocity, suitable for global studies.

Establishing statistical characteristics of random heterogeneities under various global localities and dif-ferent tectonic units is useful to understand mantle con-vection and mixing. In this paper, we document results under Tibet using the Hi-CLIMB (Himalayan–Tibetan Continental Lithosphere during Mountain Building) data (Nábelěk et al., 2009).

2.   Data and method

In this paper, we perform the transverse coher-ence function (TCF) analysis; the assumption is that the incident wave is a plane wave. We used seismic data from the Hi-CLIMB seismic experiment (Figure 1). The Hi-CLIMB seismic line extends from the in-terior of the Tibet southwards to Nepal. The south-ernmost stations are in the Ganga sedimentary basin where thick (≥6 km) sediments are documented. For those stations, unusual amplifications of seismic waves are observed due to soft sediments. Therefore, in our analysis, we exclude data from those stations and only use stations whose elevation is greater than 4 300 m. We only use direct P waves from distant earthquakes whose depths are greater than 40 km and epicentral distances greater than 40 degrees to avoid upper mantle triplica-tions. The earthquakes are mainly in subduction zones around the west Pacific and Indonesia (Figure 1). We loosely divide the earthquakes into three groups. The nominal Japan group includes events in Alaska, Kuriles, Japan and Marianas, whose backazimuths are less than 100 degrees. The Tonga group consists of earthquakes with back azimuths larger than 100 degrees and epicen-tral distances greater than 60 degrees. The Indonesia group has earthquakes with backazimuths larger than 100 degrees and epicentral distances greater than 40 but less than 60 degrees. The TCF measures the spatial correlation of the natural logarithmic amplitude (logA) or phase between two stations at a particular frequency. For each event, we measure the logA (or phase) at 0.5 Hz and remove the linear trend to obtain the logA (or phase) fluctuations across the seismic array. These fluc-tuations are then used to form the coherence functions that depend on the distance between the stations. Co-herence functions from multiple earthquakes are then averaged and resampled using the bootstrap technique (Efron and Tibshirani, 1993) to compute the statisti-cal errors and to test for statistical significance. It is remarkable that TCFs for both logA and phase for the three separate event groups are similar (Figure 2), which strongly argues that the origin of heterogeneities giving rise to these fluctuations is below the Hi-CLIMB seismic stations. Using these logA and phase TCF data, we can invert for the heterogeneity spectrum under stations.

Figure  1.  (a) A map showing Hi-CLIMB seismic stations with elevation greater than 4 300 m; (b) Earthquakes (red stars) used in this study. White dashed circles are epicentral distances (every 30°) from a nominal point (longitude 84.9°E, latitude 31.5°N) in Tibet.

DownLoad: Full-Size Img PowerPoint

Figure  2.  Comparison between observed TCF data (dots with error bars of 90% confidence) and modeled TCFs (thick lines) using heterogeneity spectra in Figure 3 for Japan, Tonga and Indonesia earthquakes.

DownLoad: Full-Size Img PowerPoint

The unknown depth-dependent heterogeneity spectrum P(z, κ) (where z is the depth variable and κ is the horizontal wavenumber) and the observed co-herence functions can be related mathematically (Zheng and Wu, 2008). In the frequency (ω) domain, logA (u) and phase (φ) coherence functions between two plane waves identified by their incidence angles, φ₁ and φ₂ (in 3D geometry, φ includes both the back azimuth and the incident angle and it points in the same direction as the slowness vector; but in the TCF case, φ₁=φ₂), can be written as, respectively

and

where u(x) is the measurement of the logA fluctuation at station x for the plane wave. Therefore, 〈u(x₁)u(x₂)〉is the observed logA coherence function which depends on the spatial lag ρ=x₂-x₁ between station x₂ and x₁. However, in the TCF case and if the random medium is statistically isotropic (e.g., independent of direction), the coherence function only depends on the station dis-tance ρ=|ρ|=R(z). The symbol 〈·〉 is the ensemble av-erage operator over random media and J0(·) the Bessel function. L is the thickness of the heterogeneous region. The WKBJ amplitude term reads

where k(z) = ω/c(z) is the background wavenumber at depth z and η(z) is the vertical wavenumber. The phase function ϑ is defined as

where

is the familiar τ function (e.g., Buland and Chapman, 1983) which contains the kinematic information for the wave propagation; c(z) is the background velocity at depth z and p is the slowness of the plane wave. The validity of equations (1) and (2) has been verified by full-wave finite difference numerical simulation for wave propagation in random media in Zheng and Wu (2008). In equation (1) (or (2)), the left hand side is the ob-served coherence function using logA (or phase) fluctu-ation data and the right hand side is a weighted sum-mation of the spatially-filtered depth-dependent hetero-geneity spectra P(z, κ) over the depth z.

The inverse problem consists of solving for the spectrum P(z, κ) given the observed coherence func-tions. We use a vector m to denote the unknown spec-trum P(κ_i)≥0 and A for the matrix of kernels. In com-puting A, one has to pay attention to the oscillatory in-tegral. Assuming the heterogeneous medium is just one stationary layer, we can use the standard least squares method to minimize

in which d contains the observables (i.e., logA and phase coherence functions on the Earth's surface). One can also control the smoothness of the model m by mini-mizing the following objective function with positivity constraint for m (Constable et al., 1987):

in which

is the smoothing regularization operator. This is equiva-lent to solving the following problem in the least squares sense:

For a multi-layered random model, the similar inverse procedure follows.

3.   Inversion results

For each event group (Japan, Tonga and Indone-sia), we simultaneously use both the logA and the phase TCFs for that group to invert for the heterogeneity spectrum under the stations. We assume that hetero-geneities exist in a single layer and are statistically sta-tionary in space. In the inversion, this unknown hetero-geneous layer is superimposed on a background velocity model, which is the IASPEI91 model (Kennett and En-gdahl, 1991). The TCF inversion has weak constraint on the layer thickness (Wu and Xie, 1991). So we can only invert for the average statistical property of hetero-geneities for the layer. Because the phase TCF has large sensitivity at the shallow depths and the logA TCF is sensitive to deeper structures, the layer thickness can-not be arbitrary. If the layer thickness is too small, only the phase TCF can be fitted well but not the logA. On the other hand, if the layer thickness is too large, logA data can be fitted but not the phase. Earlier studies as-sumed that the heterogeneities were in the lithosphere. So we can adopt a layer thickness of 120 km. The layer thickness has a trade-off with the random heterogeneity perturbation strength. Using different layer thicknesses can only affect the root-mean-square (RMS) velocity of the heterogeneities, not the shape of the inverted spectrum. Therefore, the scale lengths of the hetero-geneities can still be inferred. In our parameterization of the spectrum, the largest length scale is 400 km and the smallest is 20 km. This is consistent with the array aperture and the average station spacing. Within this range, we invert for 61 discrete spectral values. Linear interpolation between discrete spectra is adopted in the inversion. We tried different regularization parameter μ and the optimal value is determined by the well-known L-curve, which corresponds to a smooth model that can fit the data (Figure 2). We also inspected other μ values around the optimal value. The results are stable with respect to the choice of μ. The inverted models show that the spectral amplitude decreases from the largest scale to about 80 km (Figure 3). This part of the spec-trum explains the long-wavelength features for both the logA and phase TCFs. However, in order to fit the short-wavelength rapid fluctuations in the logA TCF, small-scale heterogeneities must be used. Heterogeneities on the scales of ∼50 km and ∼25 km are also inverted (Fig-ure 3). We have to clarify here that the heterogeneity scale is meant to be the spatial wavelength of the het-erogeneities. Over the spatial distance of a wavelength, there are a positive anomaly and a negative anomaly. Random realizations show that these inverted spectra (Figure 3), if isotropic, correspond to RMS velocity fluc-tuations of ∼1.9%, 1.8% and 1.6%, for the Japan, Tonga and Indonesia groups, respectively, with respect to the background velocity model.

Figure  3.  Inverted heterogeneity spectra at 0.5 Hz. For each event group, we simultaneously model both the phase and the logA TCF data.

DownLoad: Full-Size Img PowerPoint

4.   Discussion and conclusions

Regarding inversion of small-scale statistical het-erogeneities, early studies assumed that the hetero-geneities are Gaussian in statistics (e.g., Aki, 1973) or power-law (Bataille and Flatté, 1988; Flatté and Wu, 1988; Sato and Fehler, 1998) by a parametric approach. However, we directly invert for the average heterogene-ity spectrum using both logA and phase data with-out assuming a specific form for the spectrum a pri-ori. Our results showed that the average heterogeneity spectrum under the Tibet is not a single-scale Gaussian type spectrum but includes several length scales. This is a robust observation despite that the heterogeneity layer thickness is only approximately constrained. It is also plausible that under Tibet several layers of random heterogeneities exist, with each having a distinct spec-trum. However, given the depth resolving power of the TCF, inverting multi-layer spectra has not been possi-ble. However with JTACF data this is feasible (Zheng and Wu, 2008).

Global deterministic seismology revealed signifi-cant velocity variations in the mantle on the length scale of thousands of km, which are frequently inter-preted as cold downwellings and hot upwellings in term of plates and plumes. Due to limitation in the method-ology and data coverage, such deterministic methods cannot produce information about small-scale hetero-geneities. However information about length scales on the order of tens of km is very useful because they cor-respond to spatial dimensions of familiar and pervasive tectonic units such as subducted oceanic plates, sub-ducted oceanic crusts, depleted harzburgite layers at the mid ocean ridges, etc. This represents a gap in scale. Our method in principle can provide such information even for the deep lower mantle. Small-scale heterogene-ity spectra may yield information how the mantle is convecting. For large-scale heterogeneities (> 500 km), global and regional tomography showed that the spec-trum is power-law and essentially red (Chevrot et al., 1998), meaning the mantle convection primarily destroy large scale heterogeneities. However, in our results, we observe that the spectral amplitude decreases at the large-scale end. But the same decaying behavior stops at length scales of tens of km (Figure 3). This may in-dicate that the mantle convection is not efficient in de-stroying small-scale heterogeneities.