# The ZH ratio method for long-period seismic data: inversion for S-wave velocity structure

Tomoko Yano, T. Tanimoto, Luis Rivera

► To cite this version:

Tomoko Yano, T. Tanimoto, Luis Rivera. The ZH ratio method for long-period seismic data: inversion for S-wave velocity structure. Geophysical Journal International, 2009, 179 (1), pp.413-424. 10.1111/j.1365-246X.2009.04293.x. hal-02963281

HAL Id: hal-02963281

https://hal.science/hal-02963281

Submitted on 27 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The ZH ratio method for long-period seismic data: inversion for S-wave velocity structure

Tomoko Yano,<fnref n="1" /> T. Tanimoto<fnref n="1" /> and L. Rivera<fnref n="2" />

Accepted 2009 June 10. Received 2009 May 27; in original form 2008 April 5

## SUMMARY

The particle motion of surface waves, in addition to phase and group velocities, can provide useful information for S-wave velocity structure in the crust and upper mantle. In this study, we applied a new method to retrieve velocity structure using the ZH ratio, the ratio between vertical and horizontal surface amplitudes of Rayleigh waves. Analysing data from the GEOSCOPE network, we measured the ZH ratios for frequencies between 0.004 and 0.05 Hz (period between 20 and 250 s) and inverted them for S-wave velocity structure beneath each station. Our analysis showed that the resolving power of the ZH ratio is limited and final solutions display dependence on starting models; in particular, the depth of the Moho in the starting model is important in order to get reliable results. Thus, initial models for the inversion need to be carefully constructed. We chose PREM and CRUST2.0 in this study as a starting model for all but one station (ECH). The eigenvalue analysis of the least-squares problem that arises for each step of the iterative process shows a few dominant eigenvalues which explains the cause of the inversion's initial-model dependence. However, the ZH ratio is unique in having high sensitivity to near-surface structure and thus provides complementary information to phase and group velocities. Application of this method to GEOSCOPE data suggest that low velocity zones may exist beneath some stations near hotspots. Our tests with different starting models show that the models with low-velocity anomalies fit better to the ZH ratio data. Such low velocity zones are seen near Hawaii (station KIP), Crozet Island (CRZF) and Djibuti (ATD) but not near Reunion Island (RER). It is also found near Echery (ECH) which is in a geothermal area. However, this method has a tendency to produce spurious low velocity zones and resolution of the low velocity zones requires further careful study. We also performed simultaneous inversions for volumetric perturbation and discontinuity-depth perturbation. While its formulation and inversion were straightforward, there seemed to be a difficult trade-off problem between volumetric perturbation and discontinuity-depth perturbation.

**Key words:** Surface waves and free oscillations; Seismic tomography; Wave propagation.

## 1 INTRODUCTION

Surface waves continue to be an excellent source of information for retrieving S-wave velocity structure in the crust and upper mantle. However, the analysis of surface waves has largely been limited to the use of phase and group velocities. We have been developing a complementary approach that uses the particle motion of Rayleigh waves to retrieve shear wave velocity structure. We have published an observational technique to recover the ratios between the vertical and horizontal motions, hereafter referred to as the ZH ratio (Tanimoto & Rivera 2008). Such an approach actually has a long history (e.g. Boore & Nafi Toksöz 1969) but has made little progress

in the last 40 yr for long-period data. This lack of progress is in part due to the scarcity of good three three-component seismic data as well as the difficulties in measuring particle motions from real data (e.g. Ferreira & Woodhouse 2006).

The use of particle motion data in Rayleigh waves has been popular in shallow crustal studies (e.g. Hinzen et al. 2003; Okada 2003; Scherbaum 2003). Our approach is similar but we incorporate the effects of gravity in the analysis as it is not negligible for low frequency Rayleigh waves of about 10 mHz and below.

In this paper, we (1) invert ZH ratio data from GEOSCOPE stations for S-wave velocity structure beneath stations and (2) analyse the strengths and limitations of this approach. Since we discussed

<footnote n="1">¹Institute for Crustal Studies and Department of Earth Science, University of California, Santa Barbara, CA 93106, USA.
E-mail: tomoko@crustal.ucsb.edu</footnote>

<footnote n="2">²EOST-IPGS, 5 rue Rene Descartes, F67084 Strasbourg Cedex, France</footnote>
many aspects of this observational technique in Tanimoto & Rivera (2008), we will mainly focus on the inversion of the ZH ratios for S-wave velocity structures.

Our plan for this paper is as follows. In Section 2, before delving into the inversion analysis for S-wave velocity structure, we will discuss our views on the effects of anisotropy in the ZH ratios and why we do not pursue it in this paper. In Section 3, various aspects of inversion for S-wave velocity structure will be discussed. In Section 4, we will present our results for six GEOSCOPE stations and discuss major features in the results. Some stations near hotspots show low-velocity zones. However, since our analysis of the inversion shows that the results have a strong dependence on the initial model, we perform a test in Section 5. In Section 6, we present results of inversion including perturbations to depth of discontinuities. We will conclude in Section 7.

## 2 ANISOTROPY

The proposed approach in this paper is to collect Rayleigh wave data from many earthquakes at a seismic station and measure the ZH ratio, $η$, between vertical and horizontal amplitudes as a function of frequency

$$ \eta(f) = \frac{Z(f)}{H(f)}, \qquad (1) $$

where $f$ is the frequency, $Z(f)$ is the vertical amplitude and $H(f)$ is the horizontal amplitude. Then, this quantity $η(f)$, which was measured for a certain frequency band, is inverted for S-wave velocity.

Since the crust and the mantle have been known to be anisotropic, it is of some interest to examine the effects of anisotropy on the ZH ratio. The parameter $η(f)$ is the ratio of the two eigenfunctions of Rayleigh waves at the surface (vertical to horizontal). Since eigenfunctions vary with azimuth in anisotropic media, the ZH ratios should also show azimuthal variations.

In order to examine such azimuthal variations, we plotted $η$ as a function of azimuth. Figs 1(a)–(c) show the results for three stations, Kipapa (KIP), Crozet Island (CRZF) and Reunion Island (RER). The abscissas are the azimuths of propagation direction, measured clockwise from north, and the ordinates are the observed values of $η$. For each station, measurements at four frequencies from 10 to 40 mHz are shown. The cumulative statistics from all azimuths are shown on the right-hand side of each figure. Bars show the number of data points within 0.05 intervals of $η$. The maximum count is given at the top right-hand corner of each figure.

In these figures, each solid circle corresponds to a measurement of $η$ from an earthquake larger than magnitude 6; we analysed practically all earthquakes larger than magnitude 6 for available GEOSCOPE stations from the time of their initial installation up through 2005. Time span of available data varies from station to station, but we were able to obtain continuous data at most 17 yr at (INU) and at least 11 yr at CRZF within chosen stations, ECH, INU, KIP, RER and CRZF. One of the most striking features in these figures is the scatter of measured values of $η$. This large scatter supports the statement in Ferreira & Woodhouse (2006) which showed difficulty for measuring $η$. However, as we pointed out in Tanimoto & Rivera (2008), the predominant peaks emerge in the cumulative statistics as shown on the righthand side of each figure.

Our concern here is, however, whether we can see the effects of anisotropy from these results. Typically, the magnitude of anisotropy for surface wave velocity is from 1 to 10 per cent. We should thus

<figure><img src="image_2.png" /><figcaption>Figure 1. (a) The ZH ratios as a function of azimuth at KIP for four different frequencies, 10–40 mHz. (b) As in Fig. 1(a), for staion CRZF. (c) As in Fig. 1(a), for station RER.</figcaption></figure>

expect similar variations of $η$; this may vary within a small factor, because variations in eigenfunctions may be different from those of eigenfrequencies.

It seems clear from Figs 1(a)–(c) that the scatter in $η$ measurements is too large to clearly observe anisotropy. Though there could be azimuthal variations in $η$ on the order of 1–20 per cent that may be buried in scattered points, it is clear that there exist no clean signal for systematic azimuthal dependence. Therefore, in this paper, we will not discuss the effect of anisotropy on the ZH ratio ($η$). We will simply determine the mean and the standard error from the
<figure><img src="image_3.png" /><figcaption>Figure 1. (Continued.)</figcaption></figure>

cumulative statistics and use them in the inversion for the isotropic S-wave velocity structure.

# 3 METHOD OF INVERSION

## 3.1 Basic framework

Let us define the observed ZH ratios as $η_0(f_i)$ where the subscript i denotes observation. The subscript i in $f_i$ is introduced here because we measure $η$ at discrete frequency intervals, say for $f_1, f_2, \dots, f_N$ at $N$ different frequencies.

We denote theoretical values for $η$ as $η_T$. All the calculations are made with a spherical Earth. $η_T$ is a functional of density [$ρ(r)$], P-wave velocity [$α(r)$], and S-wave velocity [$β(r)$], where r is the radius of the earth. Thus, we write theoretical values of $η$ as $η_T(f_i; ρ, α, β)$. Hereafter, if there is no confusion, we will drop the arguments for $η_T$.

The basic framework of our approach is to seek the best structure, $ρ, α$ and $β$, which minimize the following quantity $J$:

$$J = \sum_{i=1}^{N} w_i \{\eta_0(f_i) - \eta_T(f_i; \rho, \alpha, \beta)\}^2, \quad (2)$$

where $w_i$ is a weight for the ith datum. Strictly speaking, the exact form of $J$ will have terms related to damping.

This minimization problem is non-linear because the dependence of $η_T$ on $ρ, α$ and $β$ is non-linear. We thus start from a model and iteratively improve the fit between theory and data in order to minimize $J$. For this iteration process, the most important relationship is the one between $η_T(f_i; ρ, α, β)$ and $η_T(f_i; ρ + δρ, α + δα, β + δβ)$ where $δρ, δα$ and $δβ$ are perturbations to density, P-wave velocity and S-wave velocity, respectively. This relationship can be written as

$$\eta_T(f_i; \rho + \delta\rho, \alpha + \delta\alpha, \beta + \delta\beta) = \eta_T(f_i; \rho, \alpha, \beta) + \delta\eta_T(f_i; \rho, \alpha, \beta), \quad (3)$$

where $δη_T$ is the perturbation to $η_T$ due to changes in density, P- and S-wave velocity and we can write

$$\frac{\delta\eta_T}{\eta_T} = \int_0^R \left\{ K_\rho(f_i) \frac{\delta\rho}{\rho} + K_\alpha(f_i) \frac{\delta\alpha}{\alpha} + K_\beta(f_i) \frac{\delta\beta}{\beta} \right\} dr. \quad (4)$$

In this formula, the depth resolution kernels for density, P-wave velocity and S-wave velocity are denoted as $K_ρ, K_α$ and $K_β$. We computed these kernels using numerical differentiation as we described before (e.g. Tanimoto & Alvizuri 2006; Tanimoto & Rivera 2008). Tsuboi & Saito (1983) published an analytical formula for a flat-layered medium. However, since we incorporate the effects of gravity which is not negligible in low frequency band (about 10 mHz and below), we simply used our numerical differentiation scheme.

Numerical examples of depth sensitivity kernels are shown in Fig. 2, which were computed for the PREM model (Dziewonski & Anderson 1981) with modifications to the ocean; the oceanic layer was replaced by the crustal values from the layer below. They show, for example, that the data $η$ are most sensitive to S-wave velocities in the upper 250 km at a frequency of 5 mHz (a period of 200 s) and the upper 70–80 km at 20 mHz (period 50 s). The effects of density and those of P-wave velocity are not negligible, however. The depth extent of sensitivity for $η$ is also much shallower than that for phase or group velocity of Rayleigh waves at the same frequency. Typically, the depth extent is only half of phase and group velocity kernels.

Since the data are most sensitive to S-wave velocity and it is not possible here to resolve density and P-wave velocity independently, we reduce the number of independent parameters by using

<figure><img src="image_4.png" /><figcaption>Figure 2. Depth sensitivity kernels for 0S50 (top left-hand panel), 0S100 (top right-hand panel), 0S200 (bottom left-hand panel), 0S300 (bottom right-hand panel) computed from the model PREM with modifications to the ocean layer.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Inversion results with different $c_\alpha$ at Inuyama (INU). (top panel) Observed ZH ratios (solid circles) and $1 - \sigma$ standard errors (vertical bars) as well as theoretical values by the initial model, and two final models with different $c_\alpha$ are shown. (bottom panel) The starting S-wave velocity model (red dashed line), the inversion results with $c_\alpha = 0.6$ (blue) and the inversion results with $c_\alpha = 0.8$ (green). (b) As in Fig. 3(a) except this result is for Reunion Island (RER).</figcaption></figure>

the relations

$$ \frac{\delta\rho}{\rho} = c_{\rho} \frac{\delta\beta}{\beta}, \quad \frac{\delta\alpha}{\alpha} = c_{\alpha} \frac{\delta\beta}{\beta}, \qquad (5) $$

where $c_\rho$ and $c_\alpha$ are proportionality constants we assign. This reduces the number of terms in (4) from three to one as

$$ \frac{\delta\eta_T}{\eta_T} = \int_0^R K'_\beta(f_i) \frac{\delta\beta}{\beta} dr, \qquad (6) $$

where

$$ K'_{\beta}(f_i) = c_{\rho} K_{\rho}(f_i) + c_{\alpha} K_{\alpha}(f_i) + K_{\beta}(f_i). \qquad (7) $$

The laboratory value for $1/c_\alpha$, often denoted as R, varies between 0.9 and 1.6 (Anderson & Isaac 1995; Masters et al. 2000). The anharmonic effects can raise this upper limit to 1.8 or even larger (Karato 1993; Karato et al. 1999), especially for the lower mantle. Since our study focuses on the crust- and upper-mantle region, we simply use $c_\alpha = 0.8$ ($R \approx 1.25$), about the mean of the laboratory value, and 0.4 for the value of $c_\rho$. The effects of using different proportionality constants for $c_\rho$ and $c_\alpha$ are not necessarily negligible but are not large. The comparison of results using a different ratio between P- and S-wave perturbation are shown in Figs 3(a) and (b). Fig. 3(a) is for the station INU in Japan; the observed values for $\eta$ with standard errors are shown at top by solid circles and vertical bars. The starting S-wave velocity model is shown as dashed (red) line in the bottom panel and its theoretical values for $\eta$ of the starting model are also shown as dashed (red) lines in the top panel. The inverted S-wave velocity structure with $c_\alpha = 0.6$ is shown as a solid (blue) line and that with 0.8 is shown as a solid (green) line. For density, we used $c_\rho = 0.4$ for both cases. The same analysis for another station RER is shown in Fig. 3(b). In both cases (Figs 3a and b), the final S-wave velocity models are quite similar, although there are slight deviations that can be recognized in the figures. We must regard those differences as uncertainties from the choice of the proportionality constants in the above formula. Effects from $c_\rho$ are smaller as the kernels in Fig. 2 indicate.

## 3.2 Linearization and iteration

In order to minimize J in eq. (2), we adopt a parametrization scheme that uses a stack of homogeneous layers for $\delta\beta/\beta$. This allows us to replace the integral by a summation. Then for each step of iteration, we get a least-squares problem to solve for $\delta\beta/\beta$ in each homogeneous layer. The thickness we adopted for each homogeneous layer was 5 km.

The least-squares problem that arises at each step of iteration can be written as

$$ \mathbf{A}^T \mathbf{W} \mathbf{A} \mathbf{x} = \mathbf{A}^T \mathbf{W} \mathbf{b}, \qquad (8) $$

where an element of A is an integration of depth sensitivity kernel over a homogeneous layer,

$$ a_{ij} = \int_{r_j}^{r_{j+1}} K'_{\beta}(f_i) dr, \qquad (9) $$

where i is the ith datum (frequency) and j refers to a layer number. An element in b is the difference between observation and theory $\eta_o(f_i) - \eta_T(f_i)$. The matrix W contains the weighting factor for each datum and is a diagonal matrix with $w_i(i = 1, \dots, N)$ as diagonal elements. The solution vector x contains perturbations to $\delta\beta/\beta$ and is used to update S-wave velocity; density and P-wave velocity are updated through the relations in (5).

We used simple damped least-squares with diagonal damping and first-derivative damping for S-wave velocity as differentiated with respect to depth. Typically, we selected a damping constant of about 10 per cent of the largest eigenvalue. We also tested 20 per cent damping (with more iterations) and ensured that the final solutions are the same. We stop the iterations when the variance (misfit) becomes less than 1 per cent of the original variance.

## 3.3 Dependence on initial models

One of the most important questions in non-linear optimization is the uniqueness of solutions and the dependence of solutions on their
starting models. We have found, unfortunately, that our problem
shows some dependence on the initial models.

Two examples of the inversion are shown in Figs 4(a) and (b); Fig. 4(a) is the case of Arta Tunnel in Djibouti (ATD) and Fig. 4(b) is in the case of CRZF. In both figures, the two left-hand panels show the results of inversion when PREM was the starting model. The right-hand panels show the results when the starting model was CRUST 2.0 (Bassin *et al.* 2000) in the crustal region and PREM for the mantle. Fit to the observed η are shown in top panels and variations of S-wave velocity structure are shown in bottom panels. The starting models are shown as dashed (red) lines and the final models are shown as solid (blue) lines.

In both cases, the inversion results from PREM introduce low-
velocity zones below the Moho, while the results with CRUST 2.0 do
not show such features. There are many parameters in the inversion
but it appears that the depth of the Moho plays an important role for
the final solutions. PREM places the Moho at 24.4 km while CRUST
2.0 places the Moho at 17 km for ATD and the Moho at 10.6 km
for CRZF. In our initial model, the layer of water on the surface
was replaced by the sedimental layer below. Thus, thickness of the
crust beneath some stations may be thicker than that in CRUST
2.0 because ocean depth is added to the crust. Even if CRUST 2.0
gives an oceanic layer, replacing this oceanic layer by the layer
below can be justified by the fact that some stations are located on
island.

3.4 Major information content in the ZH ratio

It is of some interest to ask why the inversion of the ZH ratio
[η(f)], shows such strong dependence on initial models. We have
done an eigenvalue–eigenvector analysis on the normal equations
in (8) which provides some perspective on this question.

Fig. 5 shows the largest six eigenvalues of the normal matrix (A^T W A) for the inversion of KIP. In this case, the starting model is PREM modified with CRUST2.0. Although eigenvalues have units, what is important here is the relative size of eigenvalues. In most cases, the largest eigenvalue is about five to six times larger than the second largest and is about 30 times larger than the third largest. This disparity clearly suggests that our solutions are dominated by a few eigenvectors which correspond to large eigenvalues. These eigenvectors have peaks in the shallow region. In this particular case, the final solutions are primarily dominated by the first eigenvector with smaller contribution from the second one and perhaps the third one. Other eigenvectors have little influence on the solutions. This dependence on a few large eigenvalues and eigenvectors explains that the information content of data is limited to near surface.

It should also be noted that the shape of eigenvector for the largest eigenvalue has a peak at a depth with large swings both on the shallower side and on the deeper side (Fig. 5). This feature has effects that velocity matching at the peak location tend to introduce large swings of opposite sign on both sides because the solution is basically dominated by the first eigenvector. Often, we see creation of spurious low-velocity zones in results, as in Figs 4(a) and (b), and such features may be due to this particular shape of eigenvector. Therefore, a claim for detection of low-velocity zones needs to be examined carefully.

Depth resolution kernels constructed for the same problem are
shown in Fig. 6. Target depths are shown by arrows in six panels in
this figure. Clearly, good resolution occurs only close to the surface,
perhaps only down to 50 km. Their shape is also constrained by the
shape of eigenvectors in Fig. 5, because a solution can be expressed

<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Effects of different initial models for ATD (Arta Tunnel in Djibuti). The result when PREM was the initial model (left-hand panel) and the result when PREM modified by CRUST2.0 (crust) was the initial model (right-hand panel). The initial models are shown as dashed (red) lines and the final models are shown as solid (blue) lines. (b) As in Fig. 4(a), for CRZF.</figcaption></figure>

as a linear combination of eigenvectors, for a single iteration of the
inversion (Parker 1977).

4 MODELS FOR VARIOUS STATIONS

In this section, we summarize our analysis for six GEOSCOPE
stations. Figs 7(a)–(f) show the fit to the ZH ratio data (top) and the
starting and final models (bottom) for each station. For both the fit
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The largest six eigenvalues of the normal matrix ($A^TWA$) for the inversion of KIP. These were computed from PREM+CRUST2.0 as the starting model.</figcaption></figure>

and S-wave velocity variations, the starting models are shown in red dash lines and the final models are shown in blue solid lines. Except for ECH, the starting model was PREM modified by CRUST 2.0. For ECH, we used the results of local study (Edel *et al.* 1975) as starting models.

One of the interesting aspects of the final S-wave velocity structure is the existence of a low-velocity zone under some stations. ATD has one at a relatively shallow depth between about 40 and 80 km. Both CRZF and KIP have a low-velocity zone at about 80–150 km. ECH also has a low-velocity zone at about 40–80 km.

Results for ATD in Fig. 7(a) is slightly different from Fig. 4(a) because it includes further iteration to improve the low frequency end of data. This station (ATD) is in Djibuti at the southern end of the Red Sea and is known to be near a hotspot. This shallow low-velocity zone may be related to the station's proximity to fractures associated with opening in the oceans in recent past. On the other hand, the depth range of low-velocity zones under CRZF and KIP, both known hotspots, are clearly deeper, perhaps due to the thicker lithosphere above thermal anomalies associated with hotspots. It may also be of some interest to note that RER, the Reunion Island, does not show an indication of a low-velocity zone. This may be due to the smaller size of mantle plumes whose signature could be wiped out by averaging over the long-wavelength scale of our data set.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Depth resolution kernels computed for KIP. Target depths are shown by arrows.</figcaption></figure>

ECH is not near a hotspot but is close to a geothermally active area in the Rhinegraben and we suspect that the low-velocity zone may be related to that activity.

## 5 TEST FOR ROBUSTNESS OF LOW-VELOCITY ZONE

While the existence of these low-velocity zones is interesting, our analyses in Section 3 indicated the limited resolving power of the ZH ratios as well as the dependence on initial model of our non-linear inversion. Therefore, we applied tests to those results with low-velocity anomalies to see the robustness of solutions. The test involved the modification of the final S-wave velocity structure with smoother velocity variations (linear gradient) which did not include low-velocity zones. Then we re-inverted data and examined how the solutions change after such modifications.

We show the results of such tests for two stations, KIP (Fig. 8a) and ECH (Fig. 8b). The solutions shown in the last section were used as the original S-wave velocity structures before modification. They are plotted as red dots and are basically the same as in Figs 8(c) and (e). For ECH, there were two local models, models 1 and 2, introduced in Fig. 8(c). In this section, the case of model 1 is used to demonstrate the robustness.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a–f) The ZH ratios (top panel) and the starting and final models (bottom panel) for each station. Starting models (red dashed line) were PREM + CRUST2.0 for all locations except ECH. Final models are shown as blue lines.</figcaption></figure>

We then modified these original structures by eliminating the low-velocity zones and replacing them with linear gradients. The models with these linear gradients are shown in both Figs 9(a) and (b) as green lines. The inversion was carried out again from these modified models as the new starting models. The final converged models after these inversions are given by blue lines. Figs 8(a) and (b) also show the fit of ZH ratios between data and theory (top) using the same colours and lines.

Both cases (KIP and ECH) show that, although the inversions were restarted from models without low-velocity zones (green lines), the re-inverted results (blue lines) clearly approached the results in the previous inversion (red dots). This conver-

gence to the original solutions does not necessarily mean that these are unique solutions and the existence of low-velocity anomalies are guaranteed in general. However it does indicate that, within the neighbourhood of the solutions in Figs 7(a)–(f), the solutions with low-velocity zones are preferred by the data since the solutions are the local minima of these initial models.

The frequency range of data that are sensitive to the difference between the linear velocity increase model (green) and the final model (blue) can be identified in the top panels of Figs 8(a) and (b). It starts from the lowest frequency in data at 0.004 Hz and extends to about 0.02 Hz.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (Continued.)</figcaption></figure>

## 6 PERTURBING DEPTHS OF DISCONTINUITY

Two referees suggested that simultaneous inversion for volumetric perturbation and depths of discontinuity may be helpful in retrieving the correct structure. In this section, the formulation

is modified to include perturbations in depths of discontinuity and the inversions were carried out for two stations, ATD and KIP. In principle, such inversions are straightforward and seem to work well. We will point out, however, that there is a difficult trade-off problem between volumetric perturbations and depths of discontinuity.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** (a) (bottom panel) S-wave velocity structure under KIP for the original solution are shown as red dots, the modified linear gradient model as green, and the solution to the second inversion (started from the green model) as blue. Low-velocity zone reappears in the second inversion. (top panel) Fit to ZH ratio data by each model. Differences between green and blue are the effects of low-velocity zone. (b) Same as Fig. 8(a) except that this is for ECH.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 9. Volumetric perturbation (left-hand panel) and boundary perturbation (right-hand panel) for two modes (frequencies) (eq. 11). In the left-hand panels, kernels for density are green, those for P-waves are red and those for S-waves are blue. Boundary perturbations are expressed in 1/km.</figcaption></figure>

## 6.1 Formulation

Inclusion of boundary perturbation terms change (6) to

$$ \frac{\delta \eta_T}{\eta_T} = \int_0^R K'_{\beta}(f_i) \frac{\delta \beta}{\beta} dr + \sum_k H_k(f_i) \delta h_k, \quad (10) $$

where $H_k(f_i)$ is the effect of boundary perturbation of the kth boundary and $\delta h_k$ is its perturbation. We compute $H_k$ by numerical differentiation in the unit of 1/km.

Fig. 9 shows examples of volumetric kernels and $H_k$ for two modes, approximately at frequencies 0.02 and 0.04 Hz. These are computed for a starting model at ATD. Volumetric perturbation terms are in the left-hand panels and boundary perturbations are in the right-hand panels. There are four discontinuities for the given depth range in the figure but the effects are dominated by three shallow discontinuities. The fourth one is at 6290 km.

As we consider an inversion for parameters with different dimensions awkward, we non-dimensionalize $\delta h_k$ by introducing a scaling parameter $L$. We rewrite (10) as

$$ \frac{\delta \eta_T}{\eta_T} = \int_0^R K'_{\beta}(f_i) \frac{\delta \beta}{\beta} dr + \sum_k [H_k(f_i)L] \left(\frac{\delta h_k}{L}\right). \quad (11) $$

We solve for $\delta\beta/\beta$ and $\delta h_k/L$. Obviously, inverting with (10) is equivalent to making a particular choice $L = 1$ (km) in (11) and there could be infinitely many options with different $L$. In essence, this will amount to using different relative damping parameters between volumetric terms and boundary perturbations terms.

## 6.2 Results for ATD and KIP

Here we show the results for ATD and KIP using two values of $L$, $L = 1$ and 10. In these inversions, we have 140 parameters for volumetric perturbation terms ($\delta\beta/\beta$). They are constant perturbation within each 5 km interval from the surface to depth 700 km. There are nine discontinuities in the starting model but the effects are dominated by the shallowest three discontinuities for both stations.

Improvement of fit between theory and data during twenty iterations are shown in Fig. 10. The cases with $L = 1$ are shown. The results with $L = 10$ are very similar. Variance for the 20th model is below 1 per cent of the original variance (by the starting model).

<figure><img src="image_13.png" /><figcaption>Figure 10. Improvement of fit during the iteration process. Circles are data and theoretical curves for the starting model to the 20th model are shown. The result for ATD at top and that for KIP at bottom. These are the results with $L = 1$ but those with $L = 10$ look very similar.</figcaption></figure>

Reduction of variance at each iteration step is shown in Fig. 11, which shows very similar behaviours between $L = 1$ and 10.

Retrieved structures are shown in Figs 12(a) and (b). Fig. 12(a) shows the results for ATD. S-wave velocity models in the left-hand panel and changes of three shallow discontinuities in the right-hand panel. Dash line (black) in the left panel is the starting model and the final model for $L = 1$ is shown by red and the final model for $L = 10$ is shown by blue. The final S-wave velocity profiles basically match (red and blue dash), at least on the scale of plot in Fig. 12(a) and both differ significantly from the starting model.

The same colour scheme to distinguish $L = 1$ and 10 is used in the right-hand panel and shows the results for depths of discontinuity. The results with $L = 1$ (red circles) basically show no changes in depths of discontinuity with iteration number. The depths do change for the inversion with $L = 10$; the first boundary and the third boundary (Moho) in this region (Afar) become more than 1 km deeper than those in the starting model.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Comparison of variance reduction between $L = 1$ (red) and $L = 10$ (blue). They are basically the same and reaches 1 per cent of the original at about the 20th iteration.</figcaption></figure>

The results for KIP are shown in Fig. 12(b). Similar to the results for ATD (Fig. 12a), S-wave velocity profiles in the left-hand panel are basically the same between $L = 1$ and $L = 10$. Both are again significantly different from the starting model (black dash). Depths of three discontinuities vary very little with $L = 1$ but seem to change almost up to a kilometre for the second layer with $L = 10$.

## 6.3 Solution dependence on L

While S-wave velocity profiles in Figs 12(a) and (b) are close between $L = 1$ and $L = 10$, it is troubling to find the differences in depths of discontinuities based on the choice of L. In fact, use of other values of L results in even different depths of discontinuity.

Some insights can be gained by examining the eigenvectors of A^TWA for this problem. In Figs 13(a) and (b), we show the eigenvectors associated with the four largest eigenvalues (for a model at ATD); parameters (horizontal axis) from 1 to 140 correspond to S-wave velocity and the rightmost nine parameters (141–149) correspond to discontinuities. The left-hand end corresponds to S-wave velocity at the surface and depth increases to the right-hand side up to 140 at 700 km.

Fig. 13(a) shows the case for $L = 1$ and the four largest eigenvalues decrease from 9.89 to 0.073. The eigenvectors are dominated by S-wave velocity at shallow depths. Small amplitudes can be recognized for boundary terms (from 141 to 149) but they are much smaller than S-wave velocity effects. As solutions are dominated by eigenvectors of large eigenvalues, it is not surprising to see very little depth changes for the discontinuities with $L = 1$.

Fig. 13(b) shows the case with $L = 100$. In this case, all four eigenvectors have significant amplitudes both for S-wave velocity and boundary perturbations. The inversion for this case (not shown in this paper) leads to larger perturbation to depths of discontinuities than the results with $L = 10$. In normal mode studies, this L is the radius of the Earth (6371 km). Such a choice would put even more emphasis on boundary perturbation terms.

It would be interesting to examine the effects from different choice of L, but getting different solutions with different L is clearly a problem as we do not know how to pick the right L.

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 12.** (a) (left-hand panel) Results for ATD. S-wave velocity results are shown. Two results for $L = 1$ (red) and $L = 10$ (blue dash) match closely. The starting model (black dash) is significantly different. (right-hand panel) Changes for the three shallow discontinuities. For $L = 1$, depths of discontinuity move very little. For $L = 10$, depth changes can be as large as 1 km in the final model at the right-hand end. (b) Same with Fig. 12(a) except that this is for KIP. Again S-wave velocity profiles match closely between $L = 1$ and $L = 10$ but depths of discontinuities differ about 1 km for the second boundary.</figcaption></figure>

## 7 CONCLUSION

Using the ZH ratios measured for some GEOSCOPE stations, we inverted for S-wave velocity structure. We found that the ZH ratio is sensitive to shallow depths and provides complementary information to phase and group velocities. While this approach is useful, the non-linear inversion often showed dependence on the initial model. The depth of Moho seems important, since different Moho depths tend to lead to different results. Careful selection of initial models seems very important for the success of this inversion.
<figure><img src="image_17.png" /><figcaption>**Figure 13.** (a) Eigenvectors of the four largest eigenvalues with $L = 1$. Parameters from 1 to 140 are related to S-wave velocity. Parameters from 141 to 149 are discontinuity-depth perturbation terms. As these eigenvectors basically dominate the solution, basically no changes are seen for depths of discontinuity. (b) Same with 13a except that this is for $L = 100$. Larger weighting by $L$ bring some amplitudes to boundary perturbation terms in the four largest eigenvalues. Depths of discontinuities move much more wildly in this case.</figcaption></figure>

We formulated simultaneous inversions for S-wave velocity profile and perturbations to depths of discontinuity. We introduced a scaling length $L$ that normalizes the boundary depth perturbation. The inversion is straightforward but the solutions for the depths of discontinuity depend on this scaling parameter. As we do not know how to pick the right scaling length, resolving the depths of discontinuity by the inversion of ZH ratios seems quite difficult.

The cause of the lack of resolving power was analysed by the eigenvalue-eigenvector analysis of the normal equation that arose for each iterative step. This analysis showed the dominance of a few sets of eigenvalues and eigenvectors. The second largest eigenvalue was about one-sixth of the first one and the third largest was about one-thirtieth, making only small contributions to the final so-

lutions. Therefore, the ZH ratios contain relatively small amount of information for deeper structure. They contain, however, important information for S-wave velocity structure at shallow depths.

Low-velocity zones were indicated in our solutions under ATD, CRZF, KIP and ECH but not under RER. The anomaly under ATD is shallower than those under CRZF and KIP. We tested these anomalies by inverting from modified structures and showed that as long as similar starting models are used, such as PREM modified with CRUST 2.0, structures with low-velocity zone fit data better. These results carry interesting implications but as the method may produce spurious low-velocity zones, further analysis is required to completely prove such claims.

## ACKNOWLEDGMENTS

We thank the GEOSCOPE group for excellent maintenance of the network and making data available. This study was partly supported by the Southern California Earthquake Center (SCEC). SCEC is funded by NSF Cooperative Agreement EAR-0529922 and USGS Cooperative Agreement 07HQAG0008. The SCEC contribution number for this paper is 1166.

## REFERENCES

* Anderson, O.L. & Isaac, D.G., 1995. Elastic constants of mantle minerals at high temperature, in *Mineral Physics and Crystallography*, ed. Ahrems, T.J., AGU Ref. Shelf 2.

* Bassin, C., Laske, G. & Masters, G., 2000. The current limits of resolution for surface wave tomography in North America, *EOS, Trans. Am. geophys. Un.*, **81**, 897.

* Boore, D. & Nafi Toksöz, M., 1969. Rayleigh wave particle motion and crustal structure, *Bull. seism. Soc. Am.*, **59**, 331–346.

* Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary Reference Earth Model, *Phys. Earth planet. Int.*, **25**, 297–356.

* Edel, J.B., Fuchs, K., Gelbke, C. & Prodehl, C., 1975. Deep structure of the southern Rhinegraben area from seismic refraction investigations, *J. Geophys.*, **41**, 333–356.

* Ferreira, A. & Woodhouse, J.H., 2006. Observations of long period Rayleigh wave ellipticity, *Geophys. J. Int.*, **169**, 161–169, doi:10.1111/j.1365-246X.2006.03276.x.

* Hinzen, K.-G., Weber, B. & Scherbaum, F., 2004. On the resolution of H/V measurements to determine sediment thickness, a case study across a normal fault in the lower Rhine Embayment, Germany, *J. Earthq. Eng.*, **8**, 909–926.

* Karato, S., 1993. Importance of anelasticity in the interpretation of seismic tomography, *Geophys. Res. Lett.*, **20**, 1623–1626.

* Karato, S., Karki, B. & Wentzcovitch, R., 1999. Origin of seismic velocity heterogeneity in the mantle, *EOS, Trans. Am. geophys. Un.*, **80**, 939.

* Masters, G., Laske, G., Bolton, H. & Dziewonski, A., 2000. The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle: implications for chemical and thermal structure, in *AGU Monograph Series, 'Seismology and Mineral Physics', Vol. 117, pp. 63–87, eds Karato, S., Forte, A.M., Liebermann R.C., Masters, G. & Stixrude, L., Am. Geophys. Union.

* Okada, H., 2003. *Geophysical Monograph Series, The Microtremor Survey Method*, Vol. 12, Society of Exploration Geophysicists.

* Parker, R.L., 1977. Linear inference and underparametrized models, *Rev. Geophys. Space Phys.*, **15**, 446–456.

* Scherbaum, F., Hinzen, K.-G. & Ohrnberger, M., 2003. Determination of shallow shear wave velocity profiles in the Cologne, Germany area using ambient vibrations, *Geophys. J. Int.*, **152**, 597–612.

* Takeuchi, H. & Saito, M., 1972. Seismic surface waves, in *Methods in Computational Physics*, Vol. 11, pp. 217–295, ed. Bolt, B.A., Academic Press, New York.
Tanimoto, T. & Alvizuri, C., 2006. Inversion of the HZ ratio of microseisms for S-wave velocity in the crust, *Geophys. J. Int.*, **165**, 323–335, doi:10.1111/j.1365-246X.2006.02905.x.

Tanimoto, T. & Rivera, L., 2005. Prograde Rayleigh wave particle motion, *Geophys. J. Int.*, **162**, 399–405.

Tanimoto, T. & Rivera, L., 2008. The ZH ratio method for long-period seismci data: sensitivity kernels and observational techniques, *Geophys. J. Int.*, **172**, 187–198, doi:10.1111/j.1365-246X.2007.03609.x.

Tsuboi, S. & Saito, M., 1983. Partial derivatives of Rayleigh wave particle motion, *J. Phys. Earth*, **31**, 103–113.