# Upper mantle stratification by P and S receiver functions

Véronique Farra, Lev Vinnik

► To cite this version:

Véronique Farra, Lev Vinnik. Upper mantle stratification by P and S receiver functions. Geophysical Journal International, 2000, 141, pp.699-712. 10.1046/j.1365-246X.2000.00118.x. insu-03596942

HAL Id: insu-03596942

https://insu.hal.science/insu-03596942

Submitted on 4 Mar 2022

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
# Upper mantle stratification by *P* and *S* receiver functions

Véronique Farra and Lev Vinnik*

Département de Sismologie, Institut de Physique du Globe de Paris, 4 Place Jussieu, BP89, F75252 Paris Cedex 05, France.

E-mail: farra@ipgp.jussieu.fr

Accepted 2000 January 7. Received 2000 January 6; in original form 1999 July 29

## SUMMARY

Seismic stratification of the upper mantle is investigated by applying two complementary techniques to the records of the Graefenberg array in southern Germany. The anisotropic *P* receiver function technique (Kosarev *et al.* 1984; Vinnik & Montagner 1996) is modified by using summary seismic events instead of individual events and different weighting functions instead of the same function for the harmonic angular analysis of the *SV* and *T* components of the *Pds* phases. The summary events provide better separation of the second azimuthal harmonic than the individual events. The parameters of the second harmonics of *SV* and *T* thus evaluated should be similar if they reflect the effects of azimuthal anisotropy. This can be used as a criterion to identify the anisotropy. To detect the *Sdp* phases and their azimuthal variations caused by azimuthal anisotropy we have developed a stacking technique, which can be termed the *S* receiver function technique. It includes axis rotation to separate interfering *P* and *S* arrivals, determination of the principal (*M*) component of the *S*-wave motion, deconvolution of the *P* components of many recordings by their respective *M* components and stacking of the deconvolved *P* components with weights depending on the level of noise and the angle between the *M* direction and the backazimuth of the event. Both techniques yield consistent results for the Graefenberg array. As indicated by the *P* receiver functions, the upper layer of the mantle between the Moho and 80 km depth is anisotropic with $dV_s/V_s$ around 0.03 and the fast direction close to 20° clockwise from north. The fast direction of anisotropy below this layer is around 110°. The boundary between the upper and the lower anisotropic layers is manifested by the detectable *Pds* and *Sdp* converted phases. Shear wave splitting in *SKS* is strongly dominated by azimuthal anisotropy in the lower layer (asthenosphere).

**Key words:** converted waves, mantle discontinuities, receiver functions, seismic anisotropy.

## INTRODUCTION

Seismic phases converted from *P* to *S* underneath the receiver (*Pds* in Fig. 1) are very useful in studies of the crust and mantle. These phases can be detected by deconvolution of the radial (*R*) or *SV* component by the vertical (*Z*) or *P* component and stacking of the deconvolved *R* or *SV* component of many events with moveout time corrections. The *R* component deconvolved by the *Z* component is termed the receiver function (Langston 1979). The *SV* component deconvolved by the *P* component is another kind of receiver function, and the general method for detection of the mantle *Pds* phases (Vinnik 1977) is referred to as stacking *P* receiver functions with moveout corrections.

<figure><img src="image_2.png" /><figcaption>Figure 1. Ray paths of *Pds* and *Sdp*.</figcaption></figure>

<footnote>* Also at: Institute of Physics of the Earth, Moscow, Russia.</footnote>
The classic receiver function technique assumes that the earth medium is isotropic. Signals in the *T* (transverse) component of the receiver function are explained by lateral heterogeneity of the isotropic medium, in particular, by dipping layers (Langston 1979; Savage 1998). Azimuthal anisotropy can be recognized if it can be assumed that the properties of the wavefield vary as a function of source azimuth (backazimuth) with a period *π*. This periodicity is expected if the axis of symmetry of the medium is horizontal (Kosarev *et al.* 1984; Levin & Park 1997), and it is still dominant if the axis is tilted by 10–20°. The anisotropic component of the wavefield can then be isolated by harmonic angular analysis of the *SV* and *T* components of the *Pds* converted phases (Kosarev *et al.* 1984; Vinnik & Montagner 1996; Girardin & Farra 1998; Vinnik *et al.* 1998). This analysis is performed by stacking *SV* and *T* components with weights depending on the azimuth.

The crustal *Pds* phases are followed in the receiver functions by crustal multiples of larger amplitude than the expected amplitudes of the mantle *Pds* phases. The slowness of phases converted at shallow depths in the mantle (100–250 km) is close to that of the crustal multiples that arrive in the same time interval. Due to their relative weakness and similar slowness, the *Pds* phases related to the shallow mantle usually cannot be easily separated from crustal multiples by stacking with move-out corrections, except in favourable circumstances (Bostock 1997). Consequently, the important region of the transition from the lithosphere to the asthenosphere still remains poorly sampled by seismic data.

A possible solution of this problem is to use the phases converted from *S* to *P* (*Sdp* in Fig. 1) rather than from *P* to *S*. In that case, the mantle converted phases *Sdp* arrive earlier than the crustal multiples. Observations of *S* phases converted to *P* from the Moho and discontinuities in the upper mantle and the mantle transition zone are reported in the literature (e.g. Bath & Stefansson 1966; Faber & Mueller 1980, 1984; Bock 1994; Bock & Kind 1991). Langston (1979) used *Sdp* converted phases in his evaluation of dipping interfaces. The *Sdp* phases related to the transition zone were observed only at long periods with a resolution much lower than that provided by the records of the *Pds* phases from the same discontinuities. Sacks *et al.* (1979) reported observations of *Sdp* phases converted at a depth of 250 km beneath the NORSAR array in Norway. From the polarity of the assumed *Sdp* phase they concluded that it was related to the bottom of the lithosphere. However, as noted by Farra *et al.* (1991), the *Sdp* phase described by Sacks *et al.* (1979), if identified correctly, corresponds to the *S* wave polarized as *SH*, which is impossible in the isotropic upper mantle.

It is usually difficult in individual records to distinguish between the true *Sdp* phases and other phases arriving in front of teleseismic *S*. Phases such as *SmP* and *PmP* (reflected from the Earth's surface *m* times) may arrive in the time window of interest and be mistaken for *Sdp*. Vinnik & Romanowicz (1991) interpreted some of the precursors to *S* as the *S* wave converted to *P* by scattering in the lithosphere between the source and the receiver. Another variety of 'false' *Sdp* can be produced by the superposition of *P* and *S* multiples between the Earth's surface and discontinuities in the mantle transition zone (Bock 1994).

What makes the 'false' *Sdp* phases and the true conversions potentially separable is their different dependence on the source depth, source function, focal mechanism and epicentral

distance. To separate the mantle *Sdp* phases from noise of various origins, the records of many events can be stacked, and the appropriate technique is the central subject of this study. To test the techniques, the records in a region with a reasonably well-understood deep structure are required; we found such records from the Graefenberg array in Germany (Fig. 2). This is the first broad-band digital array in the world (Harjes & Seidl 1978), and in more than 20 years of operation it has accumulated a large number of recordings.

There are strong indications of azimuthal anisotropy in the upper mantle beneath the array. Azimuthal anisotropy in the uppermost mantle in southern Germany was discovered from the observation of strong azimuthal velocity variations of the *Pn* phase propagating in the uppermost mantle (Bamford 1977). The fast direction of this anisotropy is around 20° clockwise from north. The polarization direction of the fast split wave in *SKS* at the Graefenberg (GRF) array deviates from this direction by up to almost 90° (Vinnik *et al.* 1994; Brechner *et al.* 1998), whereas for a one-layer model both directions should be the same. A mantle model consistent with both kinds of data includes an upper layer with the fast direction around 20°, and a lower layer of a much larger thickness with the fast direction around 100° (Vinnik *et al.* 1994). This model was confirmed by the harmonic angular analysis of the *T* component of the *Pds* phases by Vinnik & Montagner (1996).

Because the work of Vinnik & Montagner (1996) was published as a short paper, we present here the same data on the *Pds* phases in extended form, with important modifications.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Stations of the Graefenberg array: A1, B1 and C1 are three-component stations; the others are vertical-component stations.</figcaption></figure>
We then introduce our technique for the detection and analysis of the Sdp phases and apply it to teleseismic recordings of more than 100 teleseismic events. In the final section we summarize and discuss our findings.

## ANISOTROPY BENEATH THE GRF ARRAY FROM THE Pds PHASES

The components termed *P*, *SV* and *T* (transverse) are obtained by axis rotation of the recorded *Z*, *N-S* and *E-W* components of ground motion. The *P*-axis corresponds to the principal direction of *P*-wave particle motion in the vertical plane containing the source and the receiver. This direction is obtained by the analysis of the related covariance matrix. Our *SV*-axis is in the same plane and perpendicular to the *P*-axis. Note that our *P* direction is close to but not exactly the same as the direction of the *P*-wave ray. Our *SV* direction is close to but not exactly the same as the true *SV* direction, which is perpendicular to the *S*-wave ray. Our *SV* component is weakly contaminated by the stronger *P*-wave motion. The *T*-axis is perpendicular to the plane containing the *P*- and *SV*-axes. The adopted sign convention is as follows. The *SV*-axis has a positive projection in the direction from the source to the receiver, the *T*-axis is obtained by counterclockwise rotation of this projection by $\pi/2$, and the projection of the *P*-axis on the *Z*-axis is positive upwards.

In an isotropic laterally homogeneous earth, the *P* wave is coupled only with the *SV* wave. Consequently, *Pds* phases are polarized strictly in the vertical plane containing the source and the receiver, and their amplitudes are independent of the azimuth. In the presence of azimuthal anisotropy, an appreciable amount of energy of the *Pds* phases is contained in the *T* component, and the amplitude of the *SV* component depends on the azimuth. There are two different mechanisms by which the *T* component of the *Pds* phases is formed. If the *Pds* phase is converted from the boundary between two isotropic media but its wave path to the Earth's surface crosses an anisotropic medium, the *T* component can arise from the splitting of *SV*. This mechanism is exploited in the studies of mantle anisotropy using the *SKS* techniques (Vinnik et al. 1984; Savage 1999). If the *Pds* phase is converted from the discontinuity between anisotropic media with different anisotropies or from the discontinuity between isotropic and anisotropic media, the *T* component, like *SV*, is generated directly by conversion from *P* (Keith & Crampin 1977). In the process of further propagation, both the *SV* and the *T* components are modified by shear wave splitting. For hexagonal anisotropy with a horizontal axis of symmetry, *T* and the anisotropic component of *SV* vary with the azimuth with a period of $\pi$, and they can be extracted by harmonic angular analysis (Kosarev et al. 1984). To facilitate this analysis, the *SV*, *T* and *P* components are deconvolved by the *P* component. The deconvolution transforms the *P* waveform into a standard 'bump' and eliminates those differences between the records that depend on the different source functions and amplitudes of the *P* waves. The components with a period of $\pi$ can be extracted from *SV* and *T* by summation of the deconvolved components with suitable weights depending on the azimuth and moveout corrections. The moveout corrections account for the dependence of the traveltimes of the *Pds* phases on the ray parameter.

The data set of our study consists of a few tens of three-component records of teleseismic events (Table 1) obtained at stations A1, B1 and C1 (Fig. 2). The best data were obtained at A1, and they are described here in detail. To suppress signal-generated noise (randomly scattered waves) and microseisms, the records were low-pass filtered with a corner at around 7 s. The distribution of the events as a function of backazimuth is highly inhomogeneous (Fig. 3). To obtain a more homogeneous distribution and to facilitate the optimum separation of the second azimuthal harmonic, we form 'summary' events by stacking the records in narrow azimuthal sectors. Stacking here means summation and division by the number of records. In the sectors with a large number of records available, the number of summed records could reach 6, whereas in the sectors with weak seismicity the number of records was 1 or 2. The *P*, *SV* and *T* components of the summary events are shown in Fig. 4.

We follow Girardin & Farra (1998) and for each summary event define two weighting coefficients $W_i^T$ and $W_i^{SV}$ as

<table><caption>**Table 1.** List of events used for the analysis of the *Pds* phases.</caption><thead><tr><th>Date (m/d/y)</th><th>Time (hr:min:s)</th><th>Lat. deg</th><th>Lon. (deg)</th><th>Depth (km)</th><th>Baz. (deg)</th><th>Dist. (deg)</th></tr></thead><tbody><tr><td>03/09/77</td><td>14:27:56.20</td><td>41.7</td><td>131.1</td><td>556.4</td><td>42.3</td><td>74.8</td></tr><tr><td>09/04/77</td><td>15:40:57.29</td><td>51.2</td><td>178.4</td><td>34.0</td><td>8.2</td><td>78.9</td></tr><tr><td>08/16/79</td><td>21:31:26.29</td><td>41.8</td><td>130.8</td><td>588.0</td><td>42.4</td><td>74.6</td></tr><tr><td>08/24/79</td><td>04:26:54.20</td><td>9.0</td><td>-81.8</td><td>40.0</td><td>279.3</td><td>86.3</td></tr><tr><td>05/27/80</td><td>14:50:57.09</td><td>37.5</td><td>-116.1</td><td>13.0</td><td>322.1</td><td>82.6</td></tr><tr><td>07/05/80</td><td>20:25:25.70</td><td>41.9</td><td>77.5</td><td>28.0</td><td>73.5</td><td>45.5</td></tr><tr><td>01/22/81</td><td>19:34:40.09</td><td>38.2</td><td>142.7</td><td>20.0</td><td>36.5</td><td>82.5</td></tr><tr><td>09/04/81</td><td>11:15:13.60</td><td>10.0</td><td>124.0</td><td>644.5</td><td>66.1</td><td>96.7</td></tr><tr><td>10/25/81</td><td>03:22:16.48</td><td>18.2</td><td>-100.5</td><td>28.1</td><td>299.1</td><td>90.4</td></tr><tr><td>07/01/82</td><td>07:41:53.25</td><td>51.4</td><td>-177.5</td><td>47.7</td><td>7.1</td><td>78.8</td></tr><tr><td>07/04/82</td><td>01:20:06.80</td><td>27.9</td><td>137.0</td><td>536.0</td><td>45.9</td><td>88.9</td></tr><tr><td>07/31/82</td><td>06:29:15.52</td><td>51.8</td><td>176.1</td><td>37.0</td><td>9.5</td><td>78.1</td></tr><tr><td>01/24/83</td><td>23:09:21.43</td><td>12.9</td><td>93.6</td><td>78.0</td><td>87.3</td><td>75.3</td></tr><tr><td>02/14/83</td><td>08:10:03.60</td><td>55.0</td><td>-156.5</td><td>33.0</td><td>354.3</td><td>75.4</td></tr><tr><td>05/01/83</td><td>18:10:40.38</td><td>46.4</td><td>153.5</td><td>24.0</td><td>25.6</td><td>78.9</td></tr><tr><td>05/02/83</td><td>23:42:37.75</td><td>36.2</td><td>-118.6</td><td>10.2</td><td>322.5</td><td>84.3</td></tr><tr><td>06/09/83</td><td>12:49:03.85</td><td>40.2</td><td>139.0</td><td>30.6</td><td>38.0</td><td>79.4</td></tr><tr><td>06/09/83</td><td>18:46:00.90</td><td>51.4</td><td>-171.7</td><td>20.8</td><td>3.4</td><td>79.2</td></tr><tr><td>03/21/88</td><td>23:31:21.68</td><td>77.6</td><td>125.5</td><td>10.0</td><td>15.7</td><td>46.8</td></tr><tr><td>04/30/83</td><td>14:03:49.22</td><td>41.5</td><td>143.8</td><td>30.4</td><td>34.2</td><td>80.1</td></tr><tr><td>06/21/83</td><td>14:48:07.34</td><td>24.2</td><td>122.4</td><td>36.9</td><td>58.8</td><td>84.5</td></tr><tr><td>06/24/83</td><td>07:18:22.13</td><td>21.7</td><td>103.3</td><td>18.0</td><td>74.1</td><td>75.0</td></tr><tr><td>06/24/83</td><td>09:06:45.79</td><td>24.2</td><td>122.4</td><td>44.1</td><td>58.8</td><td>84.5</td></tr><tr><td>06/28/83</td><td>03:25:17.06</td><td>60.2</td><td>-138.6</td><td>18.5</td><td>345.6</td><td>68.2</td></tr><tr><td>07/07/83</td><td>20:35:37.36</td><td>-5.8</td><td>27.9</td><td>10.0</td><td>160.6</td><td>58.7</td></tr><tr><td>12/22/83</td><td>04:11:29.22</td><td>11.9</td><td>-13.5</td><td>11.3</td><td>217.9</td><td>42.9</td></tr><tr><td>03/06/84</td><td>02:17:21.26</td><td>29.4</td><td>138.9</td><td>457.4</td><td>43.7</td><td>88.6</td></tr><tr><td>04/20/84</td><td>06:31:10.63</td><td>50.1</td><td>148.7</td><td>582.0</td><td>26.9</td><td>74.1</td></tr><tr><td>04/24/84</td><td>04:11:29.09</td><td>30.9</td><td>138.4</td><td>403.1</td><td>43.3</td><td>87.1</td></tr><tr><td>10/26/84</td><td>20:22:21.82</td><td>39.2</td><td>71.3</td><td>33.0</td><td>80.4</td><td>43.1</td></tr><tr><td>05/01/85</td><td>13:27:56.13</td><td>-9.2</td><td>-71.2</td><td>599.9</td><td>258.4</td><td>92.1</td></tr><tr><td>05/14/85</td><td>13:24:57.83</td><td>-7.8</td><td>41.4</td><td>10.0</td><td>147.1</td><td>65.6</td></tr><tr><td>05/14/85</td><td>18:11:08.94</td><td>-7.8</td><td>41.4</td><td>10.0</td><td>147.1</td><td>65.6</td></tr><tr><td>05/16/85</td><td>14:20:25.15</td><td>-26.6</td><td>77.7</td><td>10.0</td><td>125.8</td><td>98.1</td></tr><tr><td>06/06/85</td><td>02:40:12.94</td><td>0.9</td><td>-28.4</td><td>10.0</td><td>228.0</td><td>59.2</td></tr><tr><td>10/05/85</td><td>15:24:02.26</td><td>62.2</td><td>-122.8</td><td>10.0</td><td>338.4</td><td>62.9</td></tr><tr><td>04/30/86</td><td>07:07:18.12</td><td>18.4</td><td>-100.6</td><td>26.5</td><td>300.0</td><td>90.8</td></tr><tr><td>03/21/88</td><td>23:31:21.68</td><td>77.6</td><td>125.5</td><td>10.0</td><td>15.7</td><td>46.8</td></tr><tr><td>08/21/89</td><td>01:09:06.63</td><td>11.9</td><td>41.9</td><td>15.8</td><td>135.5</td><td>45.4</td></tr><tr><td>05/12/90</td><td>04:50:09.06</td><td>49.0</td><td>141.9</td><td>610.8</td><td>31.5</td><td>72.9</td></tr><tr><td>05/20/90</td><td>02:21:59.77</td><td>5.0</td><td>32.1</td><td>7.2</td><td>151.4</td><td>47.9</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Locations of the events used for the analysis of the *Pds* phases. A few events at epicentral distances exceeding 90° are not shown. The epicentral distance increment is 15°. For the list of events see Table 1.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** *P*, *SV* and *T* components of the summary events used for the analysis of the *Pds* phases. Numbers on the left indicate backazimuth in degrees. Amplitudes of the *P* component are 1.0; amplitude scales for *SV* and *T* are shown in the lower right corners of the panels. Origin of the timescale corresponds to the maximum displacement of *P*.</figcaption></figure>
follows:

$$W_i^T(k, \psi) = \frac{\sin k(\psi - \phi_i)}{\sum_{j=1}^{n} \sin^2 k(\psi - \phi_j)}, \quad (1)$$

$$W_i^{SV}(k, \psi) = -\frac{\cos k(\psi - \phi_i)}{\sum_{j=1}^{n} \cos^2 k(\psi - \phi_j)}, \quad (2)$$

where $\phi_i$ is the backazimuth of the ith event and $\psi$ is a variable. Azimuthal anisotropy corresponds to $k=2$. The SV and T stacks are obtained as follows:

$$S_{SV}(k, t, \psi, d) = \sum_{i=1}^{n} W_i^{SV}(k, \psi) S V_i(t, d), \quad (3)$$

$$S_T(k, t, \psi, d) = \sum_{i=1}^{n} W_i^T(k, \psi) T_i(t, d), \quad (4)$$

where $S V_i(t, d)$ and $T_i(t, d)$ are the summary components in Fig. 4, with the time-shifts corresponding to the depth of conversion $d$. $S_{SV}$ and $S_T$ are calculated for $\psi$ in the range 0–360° (for $k=1$) and 0–180° (for $k=2$).

Our technique of stacking is broadly similar to that applied by Vinnik & Montagner (1986) to the same records, but it contains two new elements. First, we use summary events instead of individual events, and second, we use different expressions for $W_i^{SV}$ and $W_i^T$, instead of the same expression for both (SV and T) components. As shown by numerical modelling of the phases converted from discontinuities between the layers with hexagonal anisotropy and horizontal symmetry axes (Girardin & Farra 1998), the SV and T components vary with the azimuth as $-\cos 2(\psi_a - \phi)$ and $\sin 2(\psi_a - \phi)$, respectively, where the value of $\psi_a$ depends on the parameters of anisotropy. This value of $\psi$ provides maximum amplitudes of both $S_{SV}$ and $S_T$ for $k=2$. Similar waveforms of $S_{SV}$ and $S_T$ and their similar dependence on $\psi$ for $k=2$ are criteria for distinguishing between the effects of anisotropy and those of lateral heterogeneity.

Plots of $S_{SV}$ and $S_T$ are shown in Fig. 5 for $k=1$ and $k=2$. The assumed depth $d$ is 0 km. A remarkable feature of the plots for station A1 for $k=2$ is a signal that is recorded in the interval between 0 and 15 s. The waveforms in both components are very similar within the first 15 s, and their largest amplitudes in both components are observed in the same intervals of $\psi$ (0°–40° and 80°–120°). This similarity suggests that the signal is indeed formed mainly from the Pds phases in the anisotropic medium. The waveform of the signal is very different from the standard ‘bump’ in the P component (Fig. 4), which implies that this signal cannot be interpreted as the converted phase from a single discontinuity. The other wave group is seen in the T component between 20 and 40 s, but no counterpart of it is present in the SV component. Comparable amplitudes are observed in the first azimuthal harmonic ($k=1$), which means that azimuthal anisotropy with a horizontal axis of symmetry is not the only source of energy in the T component and of azimuthal variations in the SV component. A likely reason for the signal in the first azimuthal harmonic is lateral heterogeneity of the crust and mantle beneath and in the vicinity of the array. Relative to the signal in the second harmonic, the signal in the first harmonic should be regarded as noise. In the raw data (Fig. 4) the effects of anisotropy are combined with this noise, and their presence is not

k=1

k=2

<figure><img src="image_6.png" /><figcaption>Figure 5. Results of the harmonic angular analysis for the records of A1 shown in Fig. 4. Scale is shown in the lower right corner. Numbers on the left indicate the value of the angle $\psi$ in degrees.</figcaption></figure>

evident. Since both azimuthal anisotropy and lateral heterogeneity are weak, their effects are additive in the first-order approximation, and a significant fraction of noise is suppressed by azimuthal filtering with $k=2$.

Features similar to those found in the records of station A1 for $k=2$ are present on the seismograms of two other three-component stations of the array (Fig. 6), but the amplitudes of the second harmonics at B1 and C1 are weaker. The maximum amplitudes of the second harmonic are observed at B1 and C1 for $\psi=30°$ and $\psi=120°$. These values are slightly larger than at A1.

The wavefield observed at A1 was interpreted by comparing it with the synthetic seismograms. For the recorded P component and the assumed earth model, theoretical SV and T components can be computed with the aid of the Thomson-Haskell program (Kosarev et al. 1979). We denote the waveform of the incoming plane P wave in the isotropic half-space beneath the stack of anisotropic horizontal layers by $P_0(t)$. The components of motion recorded at the free surface can be expressed in the frequency domain as

$$\begin{bmatrix} SV(\omega) \\ T(\omega) \\ P(\omega) \end{bmatrix} = \mathbf{F}(\omega) P_0(\omega), \quad (5)$$

where $SV(\omega)$, $T(\omega)$ and $P(\omega)$ are the spectra of the corresponding components. Vector $\mathbf{F}(\omega)$ is composed of the transfer functions $SV_p(\omega)$, $T_p(\omega)$ and $P_p(\omega)$. From (5), $P_0(\omega)$ can be found as a function of the observed $P(\omega)$. The synthetic $SV(\omega)$ and $T(\omega)$ components, denoted $SV^*(\omega)$ and $T^*(\omega)$, can then
<figure><img src="image_7.png" /><figcaption>Figure 6. The same as Fig. 5 but for (a) B1 and (b) C1.</figcaption></figure>

be expressed as

$$SV^*(\omega) = \frac{SV_p(\omega)}{P_p(\omega)} P(\omega), \qquad (6)$$

$$T^*(\omega) = \frac{T_p(\omega)}{P_p(\omega)} P(\omega). \qquad (7)$$

The related functions of time are obtained by inverse Fourier transformation.

Our preferred model (Table 2) is only slightly different from that obtained by Vinnik & Montagner (1996), and includes two anisotropic layers. We assume that the anisotropy is hexagonal, with a horizontal axis of symmetry. Such anisotropy is described by the minimum number of parameters. The isotropic component of the model is defined by velocities $V_p$ and $V_s$ and density $\rho$. The anisotropic component is defined by $dV_s/V_s$, $dV_p/V_p$, $\eta$ and $\alpha$, where $dV_s$ is the difference between the fast and slow S velocities, $dV_p$ is the same as $dV_s$ but for the

$P$ velocity, $\eta$ is related to the anisotropic elastic constants by $\eta=F/(A-2L)$, and $\alpha$ is the fast direction. The amplitudes of the Pds phases are strongly sensitive only to $\delta V_s$ and $\alpha$. As shown by numerical modelling (Girardin & Farra 1998), $\psi_a$ for the Pds phase converted from the boundary between the isotropic crust and anisotropic mantle is equal to $\alpha+\pi$ in the mantle, and the value of $\alpha$ in the upper anisotropic layer of the model (20°) is well constrained by the data in Fig. 5. The value of $\delta V_s/V_s$ (0.03) in this layer is well constrained by the amplitude of the converted phase. The value of $dV_p/V_p$ (0.05) in the upper layer is close to that reported by Bamford (1977). The ratio between the adopted values of $dV_p/V_p$ and $\delta V_s/V_s$ is in the range of values usually reported for ultramafic rocks. The adopted $\alpha$ in the lower layer (110°) is in agreement with data on SKS splitting (Vinnik et al. 1994), and the adopted value of $\delta V_s/V_s$ together with the depth of the boundary between the two layers (80 km) yields a good fit between the synthetics and the observations. The adopted value of $dV_p/V_p$ in the lower anisotropic layer is taken to be similar to that in

<table><caption>Table 2. Model used to generate synthetics in Fig. 7.</caption><thead><tr><th>Depth (km)</th><th>V<sub>p</sub> (kms<sup>-1</sup>)</th><th>V<sub>s</sub> (kms<sup>-1</sup>)</th><th>ρ (gcm<sup>-3</sup>)</th><th>dV<sub>p</sub>/V<sub>p</sub></th><th>dV<sub>s</sub>/V<sub>s</sub></th><th>η</th><th>α (deg)</th></tr></thead><tbody><tr><td>0.0</td><td>5.80</td><td>3.40</td><td>2.60</td><td>0</td><td>0</td><td>1</td><td></td></tr><tr><td>27.0</td><td>6.90</td><td>3.80</td><td>2.80</td><td>0</td><td>0</td><td>1</td><td></td></tr><tr><td>31.0</td><td>7.96</td><td>4.39</td><td>3.31</td><td>0</td><td>0</td><td>1</td><td></td></tr><tr><td>31.0</td><td>8.0</td><td>4.49</td><td>3.38</td><td>0.05</td><td>0.03</td><td>1.10</td><td>20</td></tr><tr><td>80.0</td><td>8.0</td><td>4.49</td><td>3.38</td><td>0.05</td><td>0.03</td><td>1.10</td><td>20</td></tr><tr><td>80.0</td><td>8.0</td><td>4.49</td><td>3.38</td><td>0.05</td><td>0.03</td><td>1.10</td><td>110</td></tr><tr><td>250.0</td><td>8.0</td><td>4.49</td><td>3.38</td><td>0.05</td><td>0.03</td><td>1.10</td><td>110</td></tr><tr><td>250.0</td><td>8.56</td><td>4.67</td><td>3.38</td><td>0</td><td>0</td><td>1</td><td></td></tr></tbody></table>
the upper layer. The synthetic seismograms were calculated for the slowness values of the *P* waves of summary events for the *P* waveforms shown in Fig. 4. The synthetic *SV* and *T* components were stacked like the real data.

The synthetic stacks obtained (shown by solid lines in Fig. 7) for *k*=2 are remarkably similar to the stacks in Fig. 5 within the first 15–20 s. The complicated waveform of the signal in the *T* component (negative and positive extrema at 4 and 8 s for $\psi$=20°) is interpreted with our model as mainly an effect of the two superimposed *Pds* phases converted from the top and bottom of the upper anisotropic layer. The wave train in the *SV* component is longer than in the *T* component due to a stronger contribution from the crustal multiples. The multiples are affected by anisotropy when they are reflected from the upper boundary of the anisotropic upper mantle. The *Pds* phase related to the bottom of the lower anisotropic layer (at about 30 s in Fig. 7) cannot be identified in the real data, and the maximum depth of anisotropy in the real mantle remains unconstrained. The amplitudes of the synthetic stacks for *k*=1 in Fig. 7 are very small, which indicates that the azimuthal distribution of our summary events allows us to discriminate between the first and second azimuthal harmonics.

Below we will test our model with *Sdp* observations.

## S-wave receiver function technique

In the *P* receiver function technique the incoming wave is *P*, and the response contains *Pds* phases. Now we consider the incoming *S* and the response containing *Sdp* phases (Fig. 1).

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Results of the harmonic angular analysis for the synthetics (solid lines). Superimposed with dashed lines are the stacks obtained for Al shown in Fig. 5. Scale is shown in the lower right corner. Numbers on the left indicate the value of the angle $\psi$ in degrees.</figcaption></figure>

To detect the *Sdp* phases, the recorded three-component seismogram is decomposed into *P*, *SV* and *T* components. *SV* corresponds to the principal particle motion direction of the *S* wave in the vertical plane containing the source and the receiver. The *P*-axis is perpendicular to the *SV* axis in the same plane, and the *T* axis is perpendicular to this plane. Our definition of the *SV* and *P* components here is different from that used in the presentation of the *P* receiver function technique, but for practical purposes the difference is insignificant. In both cases the two axes are orthogonal. Orthogonality allows us to observe the secondary phases (either *Pds* or *Sdp*) without interference with the stronger parent phases (either *P* or *S*, respectively). The principal direction *M* of the *S* particle motion in the horizontal plane and the direction *O*, which is perpendicular to *M* in the same plane, are determined by the analysis of the covariance matrix related to this plane. $O(t)=0$ if the *S* wave is generated by a point shear dislocation and propagates in a laterally homogeneous isotropic medium. In the real Earth, $O(t)$ contains contributions from seismic source complexity, anisotropy in the source and receiver regions and the effects of lateral heterogeneity. At some distances *S* may consist of two or more differently polarized phases (*S*, *SKS* and *ScS*), and this again can contribute to $O(t)$. The azimuth $\theta$ of the direction *M* is controlled mainly by the backazimuth of the event and the focal mechanism of the seismic event. The $P(t)$, $M(t)$ and $O(t)$ components are deconvolved by $M(t)$. The deconvolution eliminates those differences between the records of different events, which depend on different waveforms and amplitudes of *S*.

For an anisotropic medium, the spectra of the *P*, *SV* and *T* components can be expressed as

$$ \begin{bmatrix} P(\omega) \\ SV(\omega) \\ T(\omega) \end{bmatrix} = \mathbf{F}(\omega) \begin{bmatrix} SV_0(\omega) \\ SH_0(\omega) \end{bmatrix}, \qquad (8) $$

where $SV_0(t)$ and $SH_0(t)$ are components of the incoming plane *S* wave in the isotropic homogeneous half-space beneath the stack of plane layers.

Matrix $\mathbf{F}(\omega)$ contains transfer functions and can be written as

$$ \mathbf{F}(\omega) = \begin{bmatrix} P_{SV}(\omega) & P_{SH}(\omega) \\ SV_{SV}(\omega) & SV_{SH}(\omega) \\ T_{SV}(\omega) & T_{SH}(\omega) \end{bmatrix}. \qquad (9) $$

The first equation of system (8) can be written as

$$ P(\omega) = P_{SV}(\omega)SV_0(\omega) + P_{SH}(\omega)SH_0(\omega). \qquad (10) $$

From the other two equations of system (8), $SV_0$ and $SH_0$ to be used in eq. (10) are expressed as

$$ \begin{bmatrix} SV_0(\omega) \\ SH_0(\omega) \end{bmatrix} = \frac{1}{SV_{SV}(\omega)T_{SH}(\omega) - T_{SV}(\omega)SV_{SH}(\omega)} \times \begin{bmatrix} T_{SH}(\omega) & -SV_{SH}(\omega) \\ -T_{SV}(\omega) & SV_{SV}(\omega) \end{bmatrix} \begin{bmatrix} SV(\omega) \\ T(\omega) \end{bmatrix}. \qquad (11) $$

In the isotropic medium, the transfer functions $SV_{SH}(\omega)$, $T_{SV}(\omega)$ and $P_{SH}(\omega)$ are zero. Thus, for the weakly anisotropic medium these transfer functions are first-order terms with
respect to anisotropy. If the second-order terms are neglected,

$$P(\omega) = P_R(\omega)SV(\omega) + P_T(\omega)T(\omega), \qquad (12)$$

where the first-order approximations of $P_R$ and $P_T$ can be written as

$$P_R(\omega) = \frac{P_{SV}(\omega)}{SV_{SV}(\omega)}, \quad P_T(\omega) = \frac{P_{SH}(\omega)}{T_{SH}(\omega)} - \frac{P_{SV}(\omega)SV_{SH}(\omega)}{SV_{SV}(\omega)T_{SH}(\omega)}. \qquad (13)$$

$SV(\omega)$ and $T(\omega)$ are related to $M(\omega)$ and $O(\omega)$ by

$$\begin{bmatrix} SV(\omega) \\ T(\omega) \end{bmatrix} = M(\omega) \begin{bmatrix} -\cos(\phi-\theta) \\ -\sin(\phi-\theta) \end{bmatrix} + O(\omega) \begin{bmatrix} \sin(\phi-\theta) \\ -\cos(\phi-\theta) \end{bmatrix}, \qquad (14)$$

where $\theta$ is the azimuth of the direction $M$ and $\phi$ is the backazimuth of the seismic event.

Assuming that $O(\omega)$ is much smaller than $M(\omega)$ and using eqs (12), (13) and (14), the first-order approximation of the ratio $P(\omega)/M(\omega)$ can be expressed as

$$\frac{P(\omega)}{M(\omega)} = -P_c(\omega)\cos(\phi-\theta) - P_s(\omega)\sin(\phi-\theta), \qquad (15)$$

where subscripts c and s stand for cos and sin, and

$$P_c(\omega) = P_R(\omega), \quad P_s(\omega) = P_T(\omega) - P_R(\omega) \frac{O(\omega)}{M(\omega)}. \qquad (16)$$

As will be demonstrated in the next section, seismic records suitable for processing are obtained from groups of closely spaced seismic events. Backazimuths of the events of one group are close to each other, whereas their focal mechanisms and the related values of $\theta$ are different. Estimates of $P_c(t)$ and $P_s(t)$, which provide the optimum approximation of the deconvolved $P$ components in a least-squares sense, can be obtained as follows. A misfit function $E$ is defined as

$$E(P_c(t), P_s(t)) = \sum_{i=1}^{N} w_i^2 [P_i(t) - P_c(t) \cos \Delta\theta_i - P_s(t) \sin \Delta\theta_i]^2, \qquad (17)$$

where $N$ is the number of events, $P_i(t)$ is the deconvolved $P$ component corresponding to event $i$, $\Delta\theta_i = \phi_i + \pi - \theta_i$ is simply related to the difference between the backazimuth of event $i$ and the azimuth of $M_i$, $w_i^2 = 1/\sigma_i^2$ and $\sigma_i^2$ is the variance of noise in the $i$th record. To evaluate $\sigma_i^2$, we assume that in the time interval between -20 and -60 s our records contain only noise. Let us denote

$$\mathbf{P}(t) = \begin{bmatrix} w_1 P_1(t) \\ \vdots \\ w_N P_N(t) \end{bmatrix}, \quad \mathbf{G} = \begin{bmatrix} w_1 \cos \Delta\theta_1 & w_1 \sin \Delta\theta_1 \\ \vdots & \vdots \\ w_N \cos \Delta\theta_N & w_N \sin \Delta\theta_N \end{bmatrix}. \qquad (18)$$

Vector $\mathbf{P}(t)$ is composed of the deconvolved $P$ components at time $t$ multiplied by the weight $w_i$. The solution $(P_c(t), P_s(t))$

<figure><img src="image_9.png" /><figcaption>Figure 8. Locations of the events used for the analysis of the Sdp phases. Epicentral distance increment is 15°.</figcaption></figure>
that minimizes the misfit function (17) is given by

$$ \begin{bmatrix} P_c(t) \\ P_s(t) \end{bmatrix} = (\mathbf{G}^\mathrm{T}\mathbf{G})^{-1}\mathbf{G}^\mathrm{T}\mathbf{P}(t), \qquad (19) $$

where $\mathbf{G}^\mathrm{T}$ is the transpose matrix of $\mathbf{G}$. It is clear that expression (19) describes the stacking of $P_i(t)$ with weights. Matrix $\mathbf{G}$ is independent of time $t$, so that the weights defined by $(\mathbf{G}^\mathrm{T}\mathbf{G})^{-1}\mathbf{G}^\mathrm{T}$ can be computed once and used for any $t$. The same weights can be used if the records are stacked with moveout corrections, which account for the difference in slowness between $S$ and $Sdp$.

$P_c(t)$ is easy to interpret. According to eqs (13) and (16), $P_c(\omega) = P_{SV}(\omega)/SV_{SV}(\omega)$. Therefore, in the time domain $P_c(t)$ is just the deconvolution of $P_{SV}(t)$ by $SV_{SV}(t)$. Physically it is the $P$ component (measured at the free surface) corresponding to the incident $SV$ in the isotropic half-space and deconvolved by the $SV$ observed at the free surface. By comparison, the expression for $P_s(t)$ is more complicated. It includes as the second term $P_R(\omega)$ multiplied by $O(\omega)/M(\omega)$. If the medium contains strong discontinuities in isotropic parameters (such as the Moho boundary), $P_R(\omega)$ and the related product can be large. $O(\omega)$ depends not only on the structure underneath the receiver, but also includes other, practically unknown, components. As a result of this, $P_s(t)$ does not depend only on the subsurface structure in the vicinity of the receiver, which is what we want to elucidate. We hope to find ways of interpreting $P_s(t)$ in the future, but the present report is focused on the analysis of $P_c(t)$.

## Observations of Sdp at the GRF array

The technique described in the preceding section was applied to the records of the GRF array. The array consists of sub-arrays A, B and C (Fig. 2). Each subarray includes one three-component station (A1, B1, C1) and several stations with only vertical seismographs. To suppress microseisms and scattered phases, all records were low-pass filtered with a corner period of 7–9 s. Beams of the vertical channels of each subarray were formed by assuming a theoretical slowness of $S$. Records of events at epicentral distances between 65° and 90° (Fig. 8) with a low noise and a large $M/O$ amplitude ratio were selected for final analysis. They were combined into groups with backazimuths around 38°, 59°, 80°, 280° and 355°. The numbers of selected events are 38, 27, 10, 31 and 15, respectively. The best data were obtained for subarray A, and they will be discussed in detail.

Fig. 9 demonstrates preliminary processing of the records of one event at subarray A. Among the first three traces (E–W, N–S and Z), Z is the stack of all vertical channels of the subarray, and the other traces are recorded at A1. Traces E–W and N–S are used to obtain traces R and T. Traces R and Z are used to obtain the $SV$ and $P$ components, where the $P$ component is practically free from the $S$-wave motion. The E–W and N–S components are used to obtain $M(t)$ and $O(t)$, the principal and minor components of $S$-wave motion, respectively. Finally, the $P$, $M$ and $O$ components are deconvolved by $M$, and the resulting traces are labelled $P_d$, $M_d$ and $O_d$. The signal in the $P_d$ component is much weaker than in $M_d$, and to extract it from noise reliably, many records are required. The $O_d$ component is much weaker than $M_d$, as required by the technique.

<figure><img src="image_10.png" /><figcaption>Figure 9. Example of processing the record of the $S$ phase of one event at subarray A. $P_d$, $M_d$ and $O_d$ are deconvolved $P$, $M$ and $O$ components. Scale for amplitudes of $P_d$, $M_d$ and $O_d$ is different from that for the other traces. E–W and N–S are the components of the record at A1. Z is the stack of the vertical components of A1–A5.</figcaption></figure>

Fig. 10 demonstrates the whole set of $P_d(t)$ with backazimuths around 280° for subarray A. The corresponding list of events is given in Table 3. A relatively strong and coherent signal in the time interval between -10 and 0 s is the $Sdp$ phase related to the Moho. The numbers on the left of Fig. 10 are values of $\Delta\theta_i$, and the signal related to the Moho is clearly seen for $\Delta\theta_i$ around 0, 180 and 360°.

The results of processing the records of subarray A with backazimuths around 280° are presented in Fig. 11. No moveout corrections were applied, which is acceptable for the $Sdp$ phases arriving with a lead time of not more than about 20 s. The $Sdp$ phase related to the Moho is expressed as a strong downward motion in $P_c$ with an amplitude of 0.13 at a time of -3.5 s. To see if anything besides this phase is present in the first 10–20 s in $P_c$, this strong phase should be suppressed. We assume that its waveform is similar to the waveform of $M(t)$. In practice, this is equivalent to the assumption that the Moho boundary is sharp. This phase can then be eliminated by subtracting from the $P_c$ trace the $M$ trace with appropriate scaling and time-shift. The residual trace labelled $P_c - m$ contains a phase with downward motion at a time of -13 s and an amplitude of -0.034. Standard errors of the estimates of $P_c(t)$ can be obtained with the method used by Tarantola & Valette (1982) by calculating the a posteriori covariance matrix $C_m$.
<table><caption>**Table 3.** $P_d$ traces of which are shown in Fig. 10.</caption><thead><tr><th>Date (m/d/y)</th><th>Time (hr:min:s)</th><th>Lat. (deg)</th><th>Lon. (deg)</th><th>Depth (km)</th><th>Dist. (deg)</th><th>Baz. (deg)</th></tr></thead><tbody><tr><td>02/22/78</td><td>06:07:37.0</td><td>14.249</td><td>-91.378</td><td>100</td><td>87.21</td><td>288.7</td></tr><tr><td>01/10/79</td><td>13:24:14.3</td><td>16.944</td><td>-93.543</td><td>156</td><td>86.45</td><td>292.0</td></tr><tr><td>06/04/79</td><td>06:26:42.7</td><td>15.684</td><td>-93.596</td><td>80</td><td>87.46</td><td>291.2</td></tr><tr><td>06/22/79</td><td>06:30:54.3</td><td>17.000</td><td>-94.609</td><td>107</td><td>87.05</td><td>292.8</td></tr><tr><td>11/05/79</td><td>01:51:12.9</td><td>17.828</td><td>-68.617</td><td>104</td><td>70.09</td><td>274.3</td></tr><tr><td>11/23/79</td><td>23:40:29.8</td><td>4.805</td><td>-76.217</td><td>108</td><td>84.71</td><td>271.2</td></tr><tr><td>06/25/80</td><td>12:04:56.9</td><td>4.437</td><td>-75.779</td><td>162</td><td>84.70</td><td>270.6</td></tr><tr><td>10/24/80</td><td>14:53:35.1</td><td>18.211</td><td>-98.240</td><td>72</td><td>88.25</td><td>296.3</td></tr><tr><td>09/14/81</td><td>12:44:29.8</td><td>18.320</td><td>-68.891</td><td>170</td><td>69.91</td><td>274.8</td></tr><tr><td>08/19/82</td><td>15:59:01.5</td><td>6.718</td><td>-82.680</td><td>10</td><td>87.45</td><td>277.3</td></tr><tr><td>03/08/83</td><td>17:06:36.5</td><td>11.007</td><td>-62.364</td><td>82</td><td>71.07</td><td>264.7</td></tr><tr><td>06/24/84</td><td>11:17:11.9</td><td>17.984</td><td>-69.338</td><td>23</td><td>70.44</td><td>274.9</td></tr><tr><td>07/13/86</td><td>09:12:10.7</td><td>16.061</td><td>-93.901</td><td>80</td><td>87.35</td><td>291.7</td></tr><tr><td>11/17/87</td><td>03:40:09.0</td><td>12.534</td><td>-87.030</td><td>76</td><td>85.82</td><td>284.3</td></tr><tr><td>03/10/88</td><td>06:17:23.3</td><td>10.402</td><td>-60.587</td><td>56</td><td>70.38</td><td>262.9</td></tr><tr><td>03/11/88</td><td>03:44:57.1</td><td>9.105</td><td>-82.985</td><td>22</td><td>85.84</td><td>279.1</td></tr><tr><td>05/06/88</td><td>14:46:17.0</td><td>11.493</td><td>-85.911</td><td>86</td><td>85.90</td><td>282.8</td></tr><tr><td>06/18/89</td><td>14:06:28.8</td><td>17.761</td><td>-68.811</td><td>62</td><td>70.27</td><td>274.4</td></tr><tr><td>09/16/89</td><td>23:20:53.2</td><td>16.497</td><td>-93.671</td><td>108</td><td>86.88</td><td>291.8</td></tr><tr><td>05/08/90</td><td>00:01:40.0</td><td>6.905</td><td>-82.622</td><td>9</td><td>87.27</td><td>277.4</td></tr><tr><td>03/01/91</td><td>17:30:26.1</td><td>10.939</td><td>-84.637</td><td>196</td><td>85.51</td><td>281.5</td></tr><tr><td>04/04/91</td><td>03:22:57.9</td><td>7.017</td><td>-78.153</td><td>32</td><td>84.30</td><td>274.1</td></tr><tr><td>05/04/91</td><td>03:42:54.5</td><td>9.542</td><td>-82.418</td><td>10</td><td>85.15</td><td>278.9</td></tr><tr><td>11/19/91</td><td>22:28:51.1</td><td>4.554</td><td>-77.442</td><td>21</td><td>85.69</td><td>271.9</td></tr><tr><td>03/07/92</td><td>01:53:37.8</td><td>10.210</td><td>-84.323</td><td>78</td><td>85.86</td><td>280.8</td></tr><tr><td>09/19/93</td><td>14:10:56.1</td><td>14.362</td><td>-93.325</td><td>18</td><td>88.32</td><td>290.2</td></tr><tr><td>02/08/95</td><td>18:40:25.3</td><td>4.104</td><td>-76.622</td><td>74</td><td>85.50</td><td>271.0</td></tr><tr><td>06/14/95</td><td>11:15:04.2</td><td>12.117</td><td>-88.339</td><td>33</td><td>86.96</td><td>285.0</td></tr><tr><td>10/21/95</td><td>02:38:57.1</td><td>16.840</td><td>-93.469</td><td>159</td><td>86.49</td><td>291.9</td></tr><tr><td>12/11/97</td><td>07:56:28.8</td><td>3.929</td><td>-75.787</td><td>178</td><td>85.09</td><td>270.3</td></tr><tr><td>02/03/98</td><td>03:02:00.2</td><td>15.883</td><td>-96.298</td><td>33</td><td>88.93</td><td>293.4</td></tr></tbody></table>

Following Tarantola (1987), this matrix can be written as

$$C_m = (G^T G)^{-1}. \qquad (20)$$

The diagonal terms of matrix $C_m$ are variances of the estimates of $P_c(t)$ and $P_s(t)$. The standard error of $P_c(t)$ thus determined for the data in Fig. 11 is 0.006. The actual errors can be larger than 0.006, because our analysis is based upon simplified assumptions on the properties of noise. It is very likely, however, that not only the signal from the Moho but also the weaker signal in $P_c - m$ is statistically significant.

The records of $P_c$ for all azimuths at subarray A are displayed in Fig. 12. Standard errors are 0.013, 0.010, 0.011 and 0.014 for 38°, 59°, 80° and 355°. The amplitude of the $Sp$ phase converted from the Moho depends strongly on the azimuth: the amplitude at 280° is twice as low as than at 355°. This difference is much larger than the standard error. In the records stripped of this phase the signal with a lead time of 13 s is seen not only at 280° but also at 355° and 38°, where its polarity is reversed relative to 280°. It is not seen in the traces corresponding to 59° and 80°. In Fig. 13 the synthetic traces are superimposed on the traces taken from Fig. 12. The synthetics are calculated for a model that differs from that in Table 2 only by the depth of the boundary between the anisotropic layers (90 km instead of 80 km in Table 2). The synthetic $P_c$ traces fit the observed traces very well: the amplitude variations of the signal converted from the Moho are very similar in the synthetics and the actual data. The signal detected in the $P_c - m$ traces is interpreted with our model as the $Sdp$ phase converted from

<figure><img src="image_11.png" /><figcaption>**Figure 10.** $P_d$ traces of the events with backazimuth around 280° at subarray A. For the list of events see Table 3. Origin of the timescale corresponds to the maximum displacement of $M_d$. The numbers on the left are values of $\Delta\theta_i$ in degrees. The numbers on the right are backazimuths in degrees and theoretical slowness in $s/°$.</figcaption></figure>

the boundary between the two anisotropic layers. The model predicts correctly the polarities of the observed signal and its disappearance from the traces for 59° and 80°. The signal with a 13 s lead time could not be detected when all available records were stacked, irrespective of the backazimuth. This

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Results of processing the records of subarray A with backazimuths around 280°. Amplitude of $M_d$ is 1.0. Amplitude scale for the other traces is shown on the right. Arrows indicate times of -3.5 s of the $Sdp$ phase in $P_c$ and -13 s in $P_c - m$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Results of processing the records of subarray A with backazimuths around 280°, 355°, 38°, 59° and 80°. Amplitude scale is shown in the lower left corner. Arrows indicate time of −3.5 s for *P*_c and −13 s for *P*<sub>c-m</sub>.</figcaption></figure>

implies that the related discontinuity has no isotropic com-
ponent. Data fairly similar to those in Fig. 12 are obtained for
subarray B (Fig. 14), while the data for subarray C (not shown)
are of lower quality.

The estimates of P_c shown in Figs 11–14 are obtained on the
assumption that the converted phases have the same slowness
as the S wave. At periods of 8–10 s this assumption is justified
for depths of conversion of not more than about 150 km.
Fig. 15 demonstrates the results of processing the records of

subarray A with moveout corrections for 280°, irrespective
of the backazimuth. The first-order moveout corrections are
calculated as a product of the differential slowness and the
differential distance. The reference distance is taken at 85°, and
the reference slowness is that of S. The theoretical differential
slowness for the *Sdp* phases is positive. At a time of −30 s in
both stacks there is a signal with a negative polarity and an
amplitude of 0.015. Considering the noise level, this signal is at
the limit of resolution. The required S-velocity contrast at the
respective discontinuity at a depth of 230 km is between 0.1
and 0.2 km s⁻¹. This observation demonstrates the potential
possibility of using our technique to detect signals from the
depth range usually out of reach for *P* receiver functions in
the continents. According to Bock (1994), the signal with a
similar lead time can be generated by multiples between the
Earth’s surface and the mantle transition zone. However, the
signal in our data is detected at a positive differential slowness
of around 0.5 s/°, which is characteristic of the *Sdp* phases
rather than of multiple reflections. The other signal with an
amplitude of 0.025 is detected at an azimuth of 280° at a time of
−52 s. This signal can be interpreted at the *Sdp* phase from the
410 km discontinuity. The points of conversion for this phase
are at a distance of about 1000 km from the station and at
different azimuths the data may sample very different regions.
This might explain why this signal is not so clearly detected in
the stack of all available records.

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Theoretical *P*<sub>c</sub> and *P*<sub>c-m</sub> components (dashed lines) superimposed on the real components taken from Fig. 12.</figcaption></figure>

Discussion

We have modified the technique of harmonic angular analysis
for the *Pds* phases and introduced a novel technique (the *S*
receiver function technique) for the analysis of the *Sdp* phases.
Our technique of harmonic angular analysis suggests a criterion
to identify the effects of azimuthal anisotropy and to distinguish
them from those of lateral heterogeneity. Graefenberg data
demonstrate that the recorded wavefield contains contributions
<figure><img src="image_15.png" /><figcaption>**Figure 14.** The same as Fig. 12 but for subarray B.</figcaption></figure>

from both azimuthal anisotropy and lateral heterogeneity, but these contributions can be separated with our *P* receiver function technique. Our *S* receiver function technique is complementary to the *P* receiver function technique. With this technique we have probably detected, among other signals of interest, the phase converted at a depth of 230 km. The data for this discontinuity may have important implications (Karato 1992).

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Results of processing the records of subarray A with moveout corrections, with backazimuths around 280° (left column) and all backazimuths (right column). Numbers on the left are values of the differential slowness in s/°. Times of the detected phases are shown by arrows. Amplitude scale is in the lower right corner.</figcaption></figure>

When applied to teleseismic records of the GRF array, both techniques yielded consistent results. Since the points of conversion of *Sdp* and *Pds* are different for the same *d* (Fig. 1), this consistency implies some degree of lateral homogeneity in the study region. A spectacular agreement between various kinds of data is observed at subarray A. The results of the harmonic angular analysis for *k*=2 are similar for *SV* and *T* (Fig. 5), which is characteristic of azimuthal anisotropy but is unlikely if the azimuthal variations are caused by scattered waves. Even minor features of the plots in Fig. 5 are correctly reproduced by the synthetics for the model with two anisotropic layers in Table 2 (Fig. 7). Azimuthal variations of the *Sdp* phase converted from the Moho are very similar in the actual data and the synthetics corresponding to the model (Fig. 13), and there is a signal that can be interpreted as the *Sdp* phase from the boundary between the two anisotropic layers. The synthetics and the actual data for this signal are broadly consistent, but the agreement is not so striking as for the signal converted from the Moho. To improve the agreement, the boundary between the two anisotropic layers was shifted to a depth of 90 km, instead of 80 km in the initial model, but the signal in the traces corresponding to 280° and 355° still arrived noticeably earlier than the theoretical signal. The reason can be seen in the different locations of the conversion points of the *Pds* and *Sdp* phases related to the same discontinuity.

There are some differences between the data for subarray A and those for the other subarrays that are not well understood. Anisotropic signals at B and C are weaker than at A (Fig. 6) and the value of $\psi_a$ is slightly different. The weakness of the signal could be explained by a weaker anisotropy in the upper layer, but this is inconsistent with the strong azimuthal dependence of the amplitude of the *Sdp* phase converted from the Moho at B (Fig. 14). The anisotropic signal in Fig. 5 is interpreted as a result of constructive interference between the *Pds* phases converted from the top and bottom of the upper anisotropic layer. The amplitude of the signal could then be reduced at B and C by a weakness of the phase converted from the bottom
of the layer. This explanation seems to be in a better agreement with the data in Fig. 14. We note that a strong dependence of the parameters of shear wave splitting in SKS on the azimuth, which is indicative of a two-layer anisotropic structure, was documented by Vinnik *et al.* (1994) only for station A1.

The most useful method for mapping azimuthal anisotropy in the continental mantle is based on observations of shear wave splitting in SKS and similar phases (Vinnik *et al.* 1984; Savage 1999). Unfortunately, this method combines an excellent lateral resolution with a negligible depth resolution. A lack of tight constraints on the depth of anisotropy often makes interpretations ambiguous and speculative. In this respect, the region of our study is remarkable in the amount of mutually consistent data on the depth of anisotropy, and from this point of view this region is unique. The fast direction of anisotropy in the uppermost layer of the mantle immediately below the Moho (around 20°) is now established by four independent techniques: (1) $P_n$ traveltimes (Bamford 1977); (2) azimuthal variations of the parameters of shear wave splitting in SKS (Vinnik *et al.* 1994; Brechner *et al.* 1998); (3) Pds converted phases (Vinnik & Montagner 1996; present study); (4) Sdp converted phases (present study). The thickness of this layer is around 50 km, and at a greater depth the fast direction of anisotropy becomes 100°–110°. This again is established by several independent teleseismic studies (Vinnik *et al.* 1994; Vinnik & Montagner 1996; Brechner *et al.* 1998; present study).

The origin of anisotropy in the upper layer is unclear. Fuchs (1977) proposed that parallelism between the direction of 20° and the strike of the Rhine Graben implies a relationship between the anisotropy and the present-day dynamics. Anisotropy is caused by deformation, but anisotropy with a fast direction of 20° is found at a depth of only 30 km, where the temperature is low and deformation is possible only under very high deviatoric stresses. Such stresses are possible in narrow active zones, but not in a large area without clear evidence of strong recent deformation. Alternatively, this anisotropy could be preserved in a frozen state from older orogenies. Considering the strength of this anisotropy it would be reasonable to expect that the direction of 20° is somehow imprinted in the crustal fabric, but no evidence of this has been found so far. It is possible that we are dealing with a remnant of the oceanic lithosphere, where anisotropy of similar strength is very common. As for the direction of 110°, it is found in the depth range of the asthenosphere and was most probably caused by recent mantle flow. The contribution of this anisotropy to the SKS splitting is a few times larger than that of lithospheric anisotropy, and this is very likely to be true of many other stable continental regions. As demonstrated by Dricker *et al.* (1999), the direction of 110° is dominant in SKS in central and eastern Europe, which implies flow on the scale of the continent, rather than a local phenomenon.

## ACKNOWLEDGMENTS

711We thank Uta Mundl for great help with collecting seismograms of the Graefenberg array. Helpful reviews by Michael Bostock and an anonymous referee are appreciated. Most of the data processing was performed with the program Seismic Handler written by Klaus Stammler. Partial support of this research by RFBR grant 98-05-64894 is acknowledged. Contribution IPGP no. 1654.

## REFERENCES

Bamford, D., 1977. Pn velocity anisotropy in a continental upper mantle, *Geophys. J. R. astr. Soc.*, **49**, 29–48.

Bath, M. & Stefansson, R., 1966. S-P conversion at the base of the crust, *Ann. Geofis.*, **19**, 119–130.

Bock, G., 1994. Multiples as precursors to S, SKS and ScS, *Geophys. J. Int.*, **119**, 421–427.

Bock, G. & Kind, R., 1991. A global study of S-to-P and P-to-S conversions from the upper mantle transition zone, *Geophys. J. Int.*, **107**, 117–129.

Bostock, M.G., 1997. Anisotropic upper-mantle stratigraphy and architecture of the Slave craton, *Nature*, **390**, 392–395.

Brechner, S., Klinge, K., Kruger, F. & Plenefish, T., 1998. Backazimuthal variations of splitting parameters of teleseismic SKS phases observed at the broadband stations in Germany, *PAGEOPH*, **151**, 305–331.

Dricker, I., Vinnik, L., Roecker, S. & Makeyeva, L., 1999. Upper-mantle flow in Eastern Europe, *Geophys. Res. Lett.*, **26**, 1219–1222.

Faber, S. & Mueller, G., 1980. Sp phases from the transition zone between the upper and the lower mantle, *Bull. seism. Soc. Am.*, **70**, 487–508.

Faber, S. & Mueller, G., 1984. Converted phases from the mantle transition zone observed at European stations, *J. Geophys.*, **54**, 183–194.

Farra, V., Vinnik, L.P., Romanowicz, B., Kosarev, G.L. & Kind, R., 1991. Inversion of teleseismic S particle motion for azimuthal anisotropy in the mantle: a feasibility study, *Geophys. J. Int.*, **106**, 421–431.

Fuchs, K., 1977. Seismic anisotropy of the subcrustal lithosphere as evidence for dynamical processes in the upper mantle, *Geophys. J. R. astr. Soc.*, **49**, 167–169.

Girardin, N. & Farra, V., 1998. Azimuthal anisotropy in the upper mantle from observations of P-to-S converted phases: application to southeast Australia, *Geophys. J. Int.*, **133**, 615–629.

Harjes, H.-P. & Seidl, D., 1978. Digital recording and analysis of broad-band seismic data at the Graefenberg (GRF)-array, *J. Geophys.*, **44**, 511–523.

Karato, S., 1992. On the Lehmann discontinuity, *Geophys. Res. Lett.*, **19**, 2255–2258.

Keith, C.M. & Crampin, S., 1977. Seismic body waves in anisotropic media: reflection and refraction at a plane interface, *Geophys. J. R. astr. Soc.*, **49**, 181–208.

Kosarev, G.L., Makeyeva, L.I., Savarensky, E.F. & Chesnokov, E.M., 1979. Influence of anisotropy beneath seismograph station on body waves, *Izv. Acad. Nauk, Fizika Zemli*, **2**, 26–37 (in Russian).

Kosarev, G.L., Makeyeva, L.I. & Vinnik, L.P., 1984. Anisotropy of the mantle inferred from observations of P to S converted waves, *Geophys. J. R. astr. Soc.*, **76**, 209–220.

Langston, C.A., 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves, *J. geophys. Res.*, **84**, 4749–4762.

Levin, V. & Park, J., 1997. Crustal anisotropy in the Ural mountains foredeep from teleseismic receiver functions, *Geophys. Res. Lett.*, **24**, 1283–1286.

Sacks, I.S., Snoke, J.A. & Husebye, E.S., 1979. Lithosphere thickness beneath the Baltic shield, *Tectonophysics*, **56**, 101–110.

Savage, M.K., 1998. Lower crustal anisotropy or dipping boundaries? Effects on receiver functions and a case study in New Zealand, *J. geophys. Res.*, **103**, 15 069–15 087.

Savage, M.K., 1999. Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting?, *Rev. Geophys.*, **37**, 65–106.

Tarantola, A., 1987. *Inverse Problem Theory: Methods for Data Fitting and Parameter Estimation*, Elsevier, Amsterdam.
Tarantola, A. & Valette, B., 1982. Generalized non-linear inverse problems solved using the least squares criterion, *Geophys. Space Phys.*, **20**, 219-232.

Vinnik, L.P., 1977. Detection of waves converted from P to SV in the mantle, *Phys. Earth planet. Inter.*, **15**, 39-45.

Vinnik, L.P. & Montagner, J.-P., 1996. Shear wave splitting in the mantle Ps phases, *Geophys. Res. Lett.*, **23**, 2449-2452.

Vinnik, L.P. & Romanowicz, B.A., 1991. Origin of pre-cursors to teleseismic S waves, *Bull. seism. Soc. Am.*, **814**, 1216-1230.

Vinnik, L.P., Kosarev, G.L. & Makeyeva, L.I., 1984. Lithospheric anisotropy as indicated by SKS and SKKS waves, *Doklady—Earth Sci. Sections*, **278**, 39-43 (published in English by Scripta Technica, New York).

Vinnik, L.P., Krishna, V.G., Kind, R., Bormann, P. & Stammler, K., 1994. Shear wave splitting in the records of the German Regional Seismic network, *Geophys. Res. Lett.*, **21**, 457-460.

Vinnik, L., Chevrot, S. & Montagner, J.-P., 1998. Seismic evidence of flow at the base of the upper mantle, *Geophys. Res. Lett.*, **25**, 1995-1998.