Computation of Large Anisotropic Seismic
Heterogeneities (CLASH)

Éric Beucler, Jean-Paul Montagner

► To cite this version:

Éric Beucler, Jean-Paul Montagner. Computation of Large Anisotropic Seismic Heterogeneities (CLASH). Geophysical Journal International, 2005, 165, pp.447-468. 10.1111/j.1365-246X.2005.02813.x. insu-03603142

HAL Id: insu-03603142

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

Submitted on 9 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
# Computation of Large Anisotropic Seismic Heterogeneities
(CLASH)

Éric Beucler<fnref n="1" />,* and Jean-Paul Montagner<fnref n="2" />

Accepted 2005 September 9. Received 2005 September 7; in original form 2003 December 12

## SUMMARY

A general tomographic technique is designed in order (i) to operate in anisotropic media; (ii) to account for the uneven seismic sampling and (iii) to handle massive data sets in a reasonable computing time. One *modus operandi* to compute a 3-D body wave velocity model relies on surface wave phase velocity measurements. An intermediate step, shared by other approaches, consists in translating, for each period of a given mode branch, the phase velocities integrated along ray paths into local velocity perturbations. To this end, we develop a method, which accounts for the azimuthal anisotropy in its comprehensive form. The weakly non-linear forward problem allows to use a conjugate gradient optimization. The Earth's surface is regularly discretized and the partial derivatives are assigned to the individual grid points. Possible lack of lateral resolution, due to the inescapable uneven ray path coverage, is taken into account through the *a priori* covariances on parameters with laterally variable correlation lengths. This method allows to efficiently separate the 2ψ and the 4ψ anisotropic effects from the isotropic perturbations. Fundamental mode and overtone phase velocity maps, derived with real Rayleigh wave data sets, are presented and compared with previous maps. The isotropic models concur well with the results of Trampert & Woodhouse. Large 4ψ heterogeneities are located in the tectonically active regions and over the continental lithospheres such as North America, Antarctica or Australia. At various periods, a significant 4ψ signature is correlated with the Hawaii hotspot track. Finally, concurring with the conclusions of Trampert & Woodhouse, our phase velocity maps show that Rayleigh wave data sets do need both 2ψ and 4ψ anisotropic terms.

**Key words:** azimuthal anisotropy, global phase velocity maps, Rayleigh surface waves, seismic tomography.

## 1 INTRODUCTION

Since the early 1980's, seismologists have made rapid progress in determining the global 3-D structure of the mantle (Dziewonski 1984; Woodhouse & Dziewonski 1984, 1986; Nataf et al. 1986; Tanimoto 1990; Montagner & Tanimoto 1991; Su et al. 1994; Li & Romanowicz 1995; Trampert & Woodhouse 1996; Grand et al. 1997; van der Hilst et al. 1997; Ekström & Dziewonski 1998; Laske & Masters 1998; Ritsema et al. 1999; Mégnin & Romanowicz 2000; Ritsema & van Heijst 2000; Debayle et al. 2005, among many others). Such improvements were made possible thanks to the impressive increase of high-quality data sets and to the incorporation of sophisticated wave propagation theories taking into account anisotropy and anelasticity or finite frequency kernels. To derive 3-D global body wave velocity models different types of data can be considered. Traveltimes (and/or waveforms) of body waves and surface waveforms are the most widely used. When focusing on approaches relying on surface wave velocity measurements—well suited for global and regional studies—the inversion procedures are obviously different from one to each other. Yet, a common feature shared by several of them is a regionalization step, which consists in translating the integrated velocities between two points, into local perturbations (e.g. Woodhouse & Dziewonski 1984; Montagner 1986; Trampert & Woodhouse 1995; Laske & Masters 1996; Ekström et al. 1997; van Heijst & Woodhouse 1999; Barmin et al. 2001;

<footnote n="1">¹Laboratoire de Planétologie et de Géodynamique, Université de Nantes, BP 92205, 2 chemin de la Houssinière, 44322 Nantes CEDEX 3, France.
E-mail: eric.beucler@univ-nantes.fr</footnote>

<footnote n="2">²Laboratoire de Sismologie Globale, I.P.G.P., 4 place Jussieu, 75252 Paris CEDEX 05, France</footnote>

<footnote>*Author to whom correspondence should be addressed.</footnote>
Trampert & van Heijst 2002; Trampert & Woodhouse 2003; Debayle & Sambridge 2004). The long-wavelength structure of the resulting global maps to date seems to be well determined, up to the harmonic degree 6. However, as pointed out by Laske (1995), the shorter-wavelength structures vary significantly according to the different parametrizations and the different groups. And, it turns out that a large part of the discrepancy can be attributed to the inversion techniques. It is well known, for instance, that heterogeneity amplitudes are not satisfactorily constrained and depend significantly upon the *a priori* constraints or *ad hoc* damping factors. Another major difference relies on the introduction or not of the azimuthal anisotropy. Addressing these two points is of great importance to further assess reliability of tomographic models.

In this paper, we present our method, entitled Computation of Large Anisotropic Seismic Heterogeneities (hereafter referred to as CLASH), to produce global isotropic and anisotropic phase velocity maps. This method is not fundamentally new since it relies on a classical azimuthal decomposition of waves horizontally propagating, but the idea of using a regular grid and assigning the partial derivatives to the individual grid points and the introduction of laterally variable *a priori* operators depending on the ray density and the azimuthal coverage make it different from others.

A brief description of the forward problem and of the major distinctive features of the CLASH are firstly given. Resolution tests using synthetic data sets, enable to assess the reliability of our approach. Finally, applications to real Rayleigh wave data sets are presented and our models are compared with previous results. The issue of the neglect of one anisotropy kind is addressed, by using synthetic data sets, and is detailed in the Appendix B.

## 2 THEORY AND INVERSE METHOD

### 2.1 Forward problem

Using the Backus' (1965) harmonic tensor decomposition and the Rayleigh's principle, Smith & Dahlen (1973, 1975) investigated the first-order effects of a small anisotropy on the propagation of Love and Rayleigh surface waves. The azimuthal dependence, due to the presence of a slight anisotropic elastic stratified half-space, is a homogeneous trigonometric polynomial of degree 4 in $\psi$, where $\psi$ denotes the azimuth of the horizontal wave vector, measured clockwise from the north. For a given point (located at $\theta, \phi$) and for a given angular frequency $\omega$, C the local phase velocity is defined as a perturbation of an isotropic reference one, denoted as $C_0$,

$$C(\theta, \phi, \omega, \psi) = C_0(\theta, \phi, \omega) + \delta C(\theta, \phi, \omega, \psi). \quad (1a)$$

The velocity perturbation $\delta C$ can be expanded as Fourier series in $\psi$,

$$\delta C(\theta, \phi, \omega, \psi) = \frac{1}{2C_0(\theta, \phi, \omega)} [A_1(\theta, \phi, \omega) + A_2(\theta, \phi, \omega) \cos 2\psi + A_3(\theta, \phi, \omega) \sin 2\psi + A_4(\theta, \phi, \omega) \cos 4\psi + A_5(\theta, \phi, \omega) \sin 4\psi]. \quad (1b)$$

The five $A_n$ coefficients depend on 21 anisotropic elastic parameters describing the medium through where the surface waves are propagating (Smith & Dahlen 1973). They are depth integral functions involving canonical harmonic elastic components, as derived in Backus (1970). These 21 harmonic components are an alternative way to describe exactly the same set as the well known 81 elastic properties, since the latter are constrained by 60 symmetry conditions. Although the $A_n$ coefficients are obviously different when considering Love or Rayleigh wave, the azimuthal dependence described in eq. (1b) holds for both cases.

In 1986, Montagner & Nataf performed the same first-order development but using Cartesian elastic coefficients. They derived the same expression as in eq. (1b) and related the five $A_n$ coefficients to 13 depth functions—which are linear combinations of elastic coefficients— involving Love or Rayleigh displacement eigenfunctions in an unperturbed medium.

For a given point, the contribution of the total perturbation referred to as $A_1$ eq. (1b) does not depend of $\psi$ and, therefore, describes the isotropic effects. The fast-axis direction of wave propagation is defined by the linear combination of the four other coefficients.

For a given angular frequency $\omega$, we consider that the jth traveltime of a surface wave train, propagating over the Earth's surface, can be expressed, by dividing the corresponding ray path into m segments, as

$$\frac{\Delta_j}{\langle C_j(\omega) \rangle} = \sum_{k=1}^{m} \frac{l_{jk}}{C_{jk}(\theta_{jk}, \phi_{jk}, \omega, \psi_{jk})}, \quad (2)$$

where $\langle C_j(\omega) \rangle$ is the integrated (or averaged) velocity along the ray path epicentral distance (in km), noted as $\Delta_j$. The way $\langle C_j(\omega) \rangle$ is determined, is beyond the scope of this article; it is only assumed this datum is correctly measured. Following a discrete approach, $C_{jk}$ represents the local velocity for each of the m segments, which length is referred to as $l_{jk}$. For the sake of clarity, geographical variables, $\theta_{jk}$ and $\phi_{jk}$, are omitted in the following. Since the anisotropic effects are taken into account, the local velocity depends on the azimuth, referred to as $\psi_{jk}$, of the local unit tangent vector pointing towards the kth segment direction.

Introducing a global reference isotropic velocity, noted as $C_0(\omega)$, eq. (2) becomes

$$\frac{\Delta_j}{\langle C_j(\omega) \rangle} = \sum_{k=1}^{m} \frac{l_{jk}}{C_0(\omega)} \left[ 1 + \frac{\delta C_{jk}(\omega, \psi_{jk})}{C_0(\omega)} \right]^{-1}, \quad (3)$$
and by using a Taylor–McLaurin expansion, we obtain

$$ \frac{C_0(\omega)}{\langle C_j(\omega) \rangle} - 1 = \frac{1}{\Delta_j} \sum_{k=1}^{m} l_{jk} \sum_{n=1}^{\infty} (-1)^n \left[ \frac{\delta C_{jk}(\omega, \psi_{jk})}{C_0(\omega)} \right]^n . \quad (4) $$

Since $\langle C_j(\omega) \rangle$ refers to an independent datum, the jth component of the data vector, hereafter referred to as **d**, is equal to the left member of eq. (4). The data space dimension is then directly given by the amount of integrated velocity measurements. In practice, *n* is set to 3 in order to prevent numerical truncating effects (see Beucler *et al.* 2003, Appendix A).

Following eq. (1b),

$$ \frac{\delta C_{jk}(\omega, \psi_{jk})}{C_0(\omega)} = \frac{1}{2C_0^2(\omega)} [A_{1jk}(\omega) + A_{2jk}(\omega) \cos 2\psi_{jk} + A_{3jk}(\omega) \sin 2\psi_{jk} + A_{4jk}(\omega) \cos 4\psi_{jk} + A_{5jk}(\omega) \sin 4\psi_{jk}], \quad (5) $$

represents a linear combination of five unknowns and hence is referred to as parameters or model, indifferently. By defining a new set of points (located as well on the Earth's surface), such as $\forall j, \forall k, \exists i : \theta_i = \theta_{jk}$ and $\phi_i = \phi_{jk}$, the forward problem can be finally written as

$\mathbf{d} = g(\mathbf{p}),$

$$ \frac{C_0(\omega)}{\langle C_j(\omega) \rangle} - 1 = \frac{1}{2\Delta_j C_0^2(\omega)} \sum_{k=1}^{m} l_{jk} \sum_{n=1}^{\infty} (-1)^n [A_{1i}(\theta_i, \phi_i, \omega) + A_{2i}(\theta_i, \phi_i, \omega) \cos 2\psi_{jk} \\ + A_{3i}(\theta_i, \phi_i, \omega) \sin 2\psi_{jk} + A_{4i}(\theta_i, \phi_i, \omega) \cos 4\psi_{jk} + A_{5i}(\theta_i, \phi_i, \omega) \sin 4\psi_{jk}]^n. \quad (6) $$

The index *i* now refers to the component of the parameter vector **p**. Thus, **p** is composed of five subvectors (**p**1... **p**5), each of them corresponding to the isotropic or the anisotropic perturbations. The parameter space dimension is entirely conditioned by the new basis gridding the Earth's surface and by the resolution degree of the derived model.

## 2.2 Inverse method

When considering eq. (6), a classical gradient least-squares optimization (Tarantola & Valette 1982b) is well suited to solve our problem. Since the latter is widely employed in many geophysical studies, only a brief presentation is given while the main particularities of the CLASH are detailed in the following.

In order to generate an improved parameter estimate ($\mathbf{p}_k$), starting from an initial parameter distribution ($\mathbf{p}_0$), we use an iterative algorithm as expressed in eq. (23) in Tarantola & Valette (1982b),

$$ \mathbf{p}_k = \mathbf{p}_0 + \mathbf{C}_p \mathbf{G}_{k-1}^T \mathbf{M}^{-1} [\mathbf{d} - g(\mathbf{p}_{k-1}) + \mathbf{G}_{k-1}(\mathbf{p}_{k-1} - \mathbf{p}_0)], \quad (7) $$

with

$$ \mathbf{M} = [\mathbf{C}_d + \mathbf{G}_{k-1} \mathbf{C}_p \mathbf{G}_{k-1}^T], \quad (8) $$

where **G** is the partial derivative matrix, $\mathbf{G}_{k-1} = \frac{\partial g(\mathbf{p}_{k-1})}{\partial \mathbf{p}_{k-1}}$. $\mathbf{C}_d$ and $\mathbf{C}_p$ are the *a priori* covariance operators on data and parameters, respectively. The weak non-linearity of the forward problem allows to choose any model as starting vector. For each iteration *k*, the compromise between modifying the model to better fit the data while keeping a good balance with the *a priori* constraints imposed on the parameters, is measured by the misfit function, whose image is

$$ S(\mathbf{p}_k) = \frac{1}{2} \left[ (\mathbf{g}(\mathbf{p}_k) - \mathbf{d})^T \mathbf{C}_d^{-1} (\mathbf{g}(\mathbf{p}_k) - \mathbf{d}) + (\mathbf{p}_k - \mathbf{p}_0)^T \mathbf{C}_p^{-1} (\mathbf{p}_k - \mathbf{p}_0) \right]. \quad (9) $$

When considering this inversion scheme as an intermediate step of a whole procedure deriving 3-D seismic velocity model, the values now treated as data sets result from previous epicentre-to-station velocity measurements. Based on smallest a priori assumptions, the *a posteriori* variance on a given integrated velocity is the only value set as *a priori* variance on the corresponding datum.

### 2.2.1 Discretization of the parameter space

In most existing tomographic methods producing global phase velocity maps, the parameter space is expanded into spherical harmonics. In practice, all phase velocity maps are truncated at a given arbitrary degree and Trampert & Snieder (1996) as Spetzler & Trampert (2003) showed that a wrong estimate in the poorly sampled regions can bias the amplitudes in the best-covered parts of the model, even at low degrees. In order to prevent truncating effects, other approaches rely on cubic B-spline basis on a triangular grid of knots (Wang & Dahlen 1995; Wang *et al.* 1998) or on the Delaunay triangulation on a sphere (Barmin *et al.* 2001; Levshin *et al.* 2005).

The CLASH relies on another type of parametrization which is, in some respect, similar to the two latter approaches. The gridding method used here offers advantages not to present some strong metric disparities between two neighbouring points located in the equatorial or in the polar regions, and to simplify the regularization of the $\mathbf{C}_p$ operator. The main idea is to evenly discretize the model space in the metric sense: the Earth's surface is divided into a number of small circles, parallel to the equator plane. They are evenly spaced out in latitude following a given input distance value. Starting from one geographical pole, grid points on each small circle are determined such that two successive knots are equally distant of the same distance value. The whole grid is rearranged and finally optimized, in the least-squares sense, such as the distance between two points on the surface, corresponding to two successive parameter vector components, is the closest to the input distance value.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Laterally variable *a priori* covariances on parameters. The Earth is regularly discretized in 2036 grid points (+). On the left, the density and the azimuthal distribution of ray paths surrounding the considered grid point (white ⊕) provide naturally a well-determined problem. Hence, the resulting covariances for this point match a Gaussian distribution over the spherical surface. On the right, since the regional ray path coverage around the considered grid point (white ⊕), using initial correlation lengths of 600 km, is not sufficient to ensure a well-determined problem, they are modified and the resulting distribution exhibits an eastward stretching. See Section 2.2.3 for details.</figcaption></figure>

All the models derived by using synthetic or real data sets and presented in this paper are built with an input distance of 500 km. As shown in Fig. 1, the Earth's surface is consequently discretized in 2036 grid points.

### 2.2.2 Relationships between data and parameters

Most tomographic methods rely so far on the great-circle approximation. Even though it was shown a long time ago that velocity heterogeneities induce focusing and defocusing effects of a ray bundle, the agreement between results of studies accounting for off-great-circle propagation (Wong 1989; Laske & Masters 1996, 1998; Wang *et al.* 1998) and others is, however, very good, provided that seismic wavelength is smaller than the heterogeneity length scale. Considering that off-great-circle effects for large length-scale tomography (Spetzler *et al.* 2001), and especially for minor arcs (Woodhouse & Wong 1986) are of second order, it is also assumed here that the seismic waves are propagating along great-circle paths joining epicentres to stations.

One main point of the CLASH relies on the relationships between the geographical distribution of ray paths (composing a given data set) and the location of the parameter components. Since these geometrical relationships determine the length and the azimuth of each segment, referred in eq. (6) to as *l*<sub>jk</sub> and *ψ*<sub>jk</sub>, respectively, they are of great importance. Assuming that the main issue is to know whether a given ray path is intersected by others, our whole approach primarily depends on the geographical ray distribution and on the amount of ray path intersections.

In order to illustrate our approach few rays are sketched in Fig. 2. The first step consists in computing all crossing points between all the ray paths. This operation is equivalent to divide a given epicentral distance in as many segments (minus one) as many times it is intersected by other rays. For each intersection, the nearest grid point (represented by +) is determined. Each grid point is secondly orthogonally projected on the corresponding ray part, defined by two successive crossing points. As shown in Fig. 2, this operation provides an exact determination of the length and the azimuth of the part of the ray attached to the nearest grid point.

When looking at Fig. 2, it becomes obvious that some ray segments are not always related to the nearest grid points, which is only due to the small amount of ray paths shown for the sake of clarity. If a few more rays, propagating in this area and crossing the considered ray path, are added, this bias rapidly vanishes. We may add that the minimum condition requirements and the computation of the $C_p$ operator ensure that such artefacts are correctly accounted for and can degrade the output model resolution but cannot affect the reliability of the results.

As a concluding remark, the major differences with other approaches are

(i) the effective length of a segment is determined only through the path intersections, which avoids smearing effects due to a preferentially oriented ray bundle, usually met in the actual geographical distributions of events and stations and

(ii) each segment is associated with only one grid point instead of several, using barycentres or linear interpolations, which greatly lightens the computations.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Sketch illustrating the relationships between the ray paths and the grid points (+). For a given ray path (bold line joining S to R) all intersections with others are computed. For each crossing point (●), the nearest grid point is determined and the part of the epicentral distance related to it, referred to as *l*<sub>jk</sub> (eq. 6), is computed. It is symbolized by the base length of the corresponding hatched triangle.</figcaption></figure>

### 2.2.3 *A priori* covariance operator on parameters

When considering the azimuthal anisotropy in its comprehensive form, the parameter vector **p** is composed of five subvectors (eq. 6). Consequently, the *a priori* covariance operator on parameters **C**<sub>p</sub> is modelled by five submatrices. Since the different subvectors describing the velocity perturbation are independent, there is no covariance between elements belonging to two different submatrices.

**C**<sub>p</sub> is usually employed to control the variation range of the model vector components when data variances are too large. In addition we impose in the CLASH, as a first priority, that **C**<sub>p</sub> ensures a well-determined problem over the whole parameter space. This implies that the correlation lengths between two components of one model subvector, can vary according to the spatial and/or the azimuthal ray density. Consequently, in a poorly sampled region, the output model resolution is degraded but the amount of linear equations is always sufficient to well solve the problem.

This atypical feature is realized in three steps: for a given initial correlation length and by using the ray path intersections, the first stage determines how many ray segments are related to a given grid point and with which azimuth. For a given grid point, when the local azimuthal distribution criterion is not matched—that is, not enough rays and/or not evenly distributed in the angular sense—the correlation length is extended, in order to increase the amount of linear equations involving the corresponding model vector component. For the two 2ψ terms for example, it means that, within the area defined by a given correlation length, at least four local azimuth values have to be evenly spaced over the [0, π] angular domain (mod π). In practice, the local azimuthal distribution criterion is the same for both the 0ψ and the 2ψ parameters and stricter for the 4ψ. The second step consists in a regularization of the correlation length distribution to avoid strong disparities between neighbouring grid points and then to generate a positive definite operator (Tarantola & Valette 1982a). The covariances between two components, belonging to the same submatrix, are finally computed by using a spherical probability density function (Fisher 1925, 1953). Two examples of covariance distributions, corresponding to two different coverages, are shown in Fig. 1.

On the left (Fig. 1a), the area surrounding the considered grid point (white ⊕) is sufficiently covered, with a good azimuthal distribution. Hence the initial correlation length does not need to be increased and the resulting covariances for this point match a Gaussian distribution over the spherical surface. On the other hand, in Fig. 1(b), since the region located in the mid-Atlantic ocean (east of Brazil) is not satisfactorily covered, correlation lengths of concerned components are largely increased. Due to the symmetry conditions, the covariances corresponding to the grid point located in Brazil (white ⊕) are modified and the resulting distribution exhibits an eastward stretching.

In comparison with approaches relying on fixed correlation lengths over the parameter space, our choice prevents from any horizontal smearing phenomenon. It ensures as well an efficient separation between 0ψ, 2ψ and 4ψ perturbations, and avoids the propagation of wrong measurements from the poorly sampled regions to the best-covered parts of the model.

Finally, it might also be noted that the regular mapping with individual grid knot assignments and with laterally variable covariances take into account the irregular seismic sampling of the Earth's surface, and, to some extent, present similarities with sophisticated parametrization schemes (Sambridge 1999; Nolet & Montelli 2005).
# 3 RESOLUTION TESTS USING SYNTHETIC DATA

## 3.1 Preliminaries

In order to assess the reliability and the robustness of the CLASH, many synthetic tests were performed but only the most relevant are detailed in this section. For the sake of clarity, the results of synthetic tests showing the effects of the neglect of anisotropic terms are detailed in Appendix B.

### 3.1.1 Methodology

Descriptions and differences between all synthetic tests are summarized in Table 1. For each experiment, the same methodology is followed:

(i) Starting from a given synthetic set of parameters (referred to as input model), velocities integrated along spherical paths, denoted as $\langle C_j(\omega) \rangle$ in eq. (6), are computed by using the forward problem equation.

(ii) The synthetic data sets are composed of 9288 epicentre-to-station integrated velocities. The coverage is identical to the ray path distribution used for the inversion of the real 100 s Rayleigh fundamental mode measurements, and is shown in Fig. 3(a). Epicentre and station locations come from the Harvard CMT catalogue and from the Federation of Digital Broad-Band Seismograph Networks (Fig. 4), respectively.

(iii) For each synthetic test, the same *a priori* variances on parameters are set. They correspond to standard deviation values of 5 per cent and 1 per cent for isotropic and anisotropic parameters, respectively.

(iv) Both isotropic and anisotropic output models are plotted using equal area Hammer projection. They are centred on the longitude $\phi = 180^{\circ}$.

### 3.1.2 Graphical representation of the output models

For both the first two synthetic test results and the models obtained by the inversions of real data sets, the models are displayed in terms of isotropic and maximum anisotropic perturbations. The isotropic perturbations are expressed in per cent with respect to the reference velocity, and are mapped with colour scales (Figs 6 and 9). The peak-to-peak amplitudes of anisotropy are represented by the segments, which correspond to the fast propagation axes. The length of segments reflects the amplitude perturbations (in per cent, with respect to the reference velocity) and the azimuths of fast axes are represented by the segment directions. When considering $2\psi$ or $4\psi$ anisotropic parameters separately (Figs 10 and 11), the amplitudes and the azimuths can be readily found, which is not feasible for the comprehensive form of anisotropy, say $2\psi+4\psi$ perturbations (Figs 12 and 13). Consequently, the amplitude perturbations and the fast-axis azimuths are numerically computed, using eq. (1b), between 0 and $180^{\circ}$, with a step of $0.25^{\circ}$, as shown in Fig. 5.

### 3.1.3 Graphical representation of the *a posteriori* variances

From a theoretical point of view, values of the *a posteriori* standard deviations of parameters are given by the square root of the diagonal elements of the *a posteriori* covariance matrix, referred to as $C_{pp}$ (Tarantola 1987, p. 189),

$$C_{pp} = C_p - C_p G_{k-1}^T M^{-1} G_{k-1} C_p, \qquad (10)$$

with $M$ as defined as in eq. (8). For all parameters but isotropic, such estimates are however meaningless, therefore, the *a posteriori* variances on $2\psi$ and $4\psi$ anisotropic parameters are turned into uncertainties in the amplitude perturbations and in the fast-axis azimuths. To this end, for a given grid point, the *a posteriori* anisotropic standard deviations are added to and subtracted from the corresponding parameter values, making 16 new anisotropic perturbation curves, as a function of $\psi$, as shown in Fig. 5. For each curve, the maximum of the perturbation

<table><caption>**Table 1.** Description and values of $\chi^2$ for the resolution tests using synthetic data sets. $p_0$ and $p_f$ refer to the starting and the final models, respectively.</caption><thead><tr><th></th><th>Input model</th><th>Inverted parameters</th><th>χ<sup>2</sup>(p<sub>0</sub>)</th><th>χ<sup>2</sup>(p<sub>f</sub>)</th><th>Variance red. (per cent)</th><th>Figure(s) n°</th></tr></thead><tbody><tr><td>Test #1</td><td>0ψ, 2ψ, 4ψ</td><td>0ψ, 2ψ, 4ψ</td><td>1.5067</td><td>0.0021</td><td>99.8</td><td>6(a), 7</td></tr><tr><td>Test #2</td><td>0ψ, 2ψ, 4ψ</td><td>0ψ, 2ψ, 4ψ</td><td>4.5269</td><td>0.0033</td><td>99.9</td><td>6(b), B1</td></tr><tr><td>Test #3</td><td>0ψ, 2ψ, 4ψ</td><td>0ψ, 2ψ</td><td>4.5269</td><td>0.0697</td><td>98.4</td><td>B2(a)</td></tr><tr><td>Test #4</td><td>0ψ, 2ψ, 4ψ</td><td>0ψ</td><td>4.5269</td><td>0.2489</td><td>94.5</td><td>B2(b)</td></tr><tr><td>Test #5</td><td>0ψ, 2ψ</td><td>0ψ, 2ψ, 4ψ</td><td>4.3117</td><td>0.0016</td><td>99.9</td><td>B3(a)</td></tr><tr><td>Test #6</td><td>0ψ, 2ψ</td><td>0ψ</td><td>4.3117</td><td>0.0986</td><td>97.7</td><td>B3(b)</td></tr><tr><td>Test #7</td><td>0ψ</td><td>0ψ, 2ψ, 4ψ</td><td>4.2026</td><td>0.0005</td><td>99.9</td><td>B4(a)</td></tr><tr><td>Test #8</td><td>0ψ</td><td>0ψ, 2ψ</td><td>4.2026</td><td>0.0005</td><td>99.9</td><td>B4(b)</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Ray path coverages used for both the computation of synthetic data sets and the inversions of real data sets. They correspond to (a) the fundamental mode at 100 s and (b) the second overtone at 150 s.</figcaption></figure>

and the corresponding azimuth are computed numerically. The so called *a posteriori* uncertainties in the amplitude perturbation and in the azimuth are finally given by the difference between the largest and the lowest maxima and by the widest angular domain between all azimuths, respectively.

### 3.1.4 Quantitative comparisons in the data space

In order to compare two different inversion results, it is possible to use a misfit estimate in the data space. For a given iteration *k*, we may define $\chi^2$, as a function of the model $\mathbf{p}_k$,
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Locations of the 141 FDSN stations (IRIS and GEOSCOPE) processed for the construction of our Rayleigh wave phase velocity maps. The same station set is used for the resolution tests using synthetic data.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** 2ψ + 4ψ anisotropic perturbations are computed numerically, using eq. (1b), between 0° and 180° (bold line), thereby providing the amplitude and the azimuth of the fast axis. Adding to and subtracting from the *a posteriori* standard deviations make 16 new curves, similarly computed (dashed lines). For a given point, when the *a posteriori* variances on both 2ψ and 4ψ parameters are (a) very small, the 16 curves almost perfectly fit the bold one. Consequently, uncertainties in the amplitude and in the azimuth are lower than 0.1 per cent and 1°, respectively. On the other hand, when the *a posteriori* variances are (b) much larger, they provide 16 various curves (dashed lines). The resulting uncertainty in the azimuth is large as well as in the amplitude.</figcaption></figure>

$$ \chi^2(\mathbf{p}_k) = \frac{1}{N} (g(\mathbf{p}_k) - \mathbf{d})^T \mathbf{C}_d^{-1} (g(\mathbf{p}_k) - \mathbf{d}), \qquad (11) $$

where N is the number of data. When iterations start from the null model, the $\chi^2(\mathbf{p}_0)$ value yields a good estimate of how the reference velocity explains the data set. In Tables 1 and 2, the values of $\chi^2$, for both the synthetic tests and the real data inversions, are reported for the starting and for the final model, noted as $\mathbf{p}_f$.

## 3.2 Synthetic test results

In the first experiment, the input model, in which the data are synthetized, is strongly unrealistic (Table 1). As shown in Fig. 6(a), it is composed of four square anomalies, two of them being purely isotropic, located on the west of Central America and eastward from Australia. The two others contain isotropic and both 2ψ and 4ψ anisotropic perturbations (from Manchuria to Australia and on the east of Japan). They are realized using step functions and, therefore, are not outwardly well suited to be recovered by an inverse method providing smooth varying parameters. After inversion of the synthetic data set, both the shapes and the amplitudes of isotropic and anisotropic heterogeneities are very correctly retrieved.

This can be seen qualitatively in Fig. 6(a) and also quantitatively when considering the variance reduction in Table 1. Only two small and very weak secondary oscillations are hardly visible. The fast-axis directions are almost identical to the input ones in the two anisotropic regions and, as a consequence, no anisotropic perturbation is observed elsewhere. As shown in Fig. 7, areas where anisotropy is reliably
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Resolution tests using synthetic data sets. The blue segments represent the fast-axis directions which are scaled to the anisotropic perturbation amplitude, expressed in per cent. For the Test #1 (a), the input model is strongly unrealistic while for the Test #2 (b), the input model is realistic but highly heterogeneous.</figcaption></figure>

detected are clearly highlighted through the *a posteriori* variances, especially the uncertainties in the fast-axis azimuths, shown in Fig. 7(c). As an example, few components of the parameter vector, corresponding to the grid points located on the east of Honshu (Japan), are not satisfactorily retrieved. It might be noted that this is the unique region where the fast-axis directions are not well recovered. The bad retrieving is due to the rapid variation of the input model anisotropic terms over a too small distance, compared to the 600 km correlation lengths. Considering the *a priori* covariance operator on parameters, such rotations of fast-axes are impossible to retrieve but, as seen in Fig. 7(c), this wrong measurement is clearly outlined. Thus, computation of the *a posteriori* variances on parameters, and consequently, the transformation into uncertainties in the anisotropic amplitudes and in the fast-axis azimuths is an efficient tool to assess whether a region is reliably constrained.

The second resolution test is realized by using a much more realistic input model. As shown in Fig. 6(b), the input model encloses smooth isotropic and both 2ψ and 4ψ anisotropic heterogeneities. In order to be consistent with the Earth's surface discretization and the initial value of the correlation lengths, the smallest length scale of heterogeneities to be recovered is set to 600 km. An interesting feature of the input model is to present much larger perturbation amplitudes, than the *a priori* variances set in C<sub>p</sub>, up to 12 per cent for the 0ψ terms and up to 3 per cent for the anisotropic terms. When comparing the value of χ<sup>2</sup>(p<sub>0</sub>) for this test with respect to the Test #1, reported in Table 1, the input model is undoubtedly highly heterogeneous. It might be noted that such a discrepancy in the data set is greater than those measured for the real Rayleigh wave data set, as reported in Table 2.

The output model is presented in Fig. 6(b) and it turns out that almost all anomalies, isotropic and anisotropic, are very well recovered. Even though such large perturbation amplitudes are not *a priori* expected, the derived model converges to the right amplitudes, after three iterations. The very good agreement between input and output models can be seen as well in the data space through the χ²(*p*<sub>f</sub>) value. As shown in the quantitative comparison maps, presented in Fig. B1, the four output anisotropic subvectors are almost identical to the corresponding input ones and only very small differences are noticeable for the 0ψ parameters. The largest distances, between isotropic input and output models, are located in a well-known poorly sampled area of the Earth's surface (around the Kerguelen archipelago) and reach values of 1.1, which corresponds approximately to a velocity difference of 1 per cent. Concerning the anisotropic parameters, the largest difference is observed in the same region for the sin 2ψ and the sin 4ψ parameters and does not exceed the value of 0.58.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Test #1 *a posteriori* variances on parameters. The variances on (a) isotropic parameters are readily deducted from the $C_{pp}$ matrix (eq. 10) whereas they are turned into uncertainties in (b) the amplitudes and (c) the azimuths of anisotropic perturbations. (b) and (c) are computed as detailed in Section 3.1.3</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Data set perturbations with respect to the PREM. For the fundamental mode and the first two overtones, the inverted phase velocity perturbations are computed as a function of period as defined in (eq. 12).</figcaption></figure>

# 4 APPLICATION TO REAL DATA SETS

We constructed global isotropic and anisotropic phase velocity maps from our fundamental and overtones Rayleigh wave data sets, measured by the roller-coaster algorithm (Beucler *et al.* 2003).

## 4.1 Preliminaries

The spheroidal mode-branch classification used hereafter is the same as Okal 's (1978). Only the first Rayleigh wave trains, with epicentral distances between 55° and 135°, are selected. As shown in Fig. 4, 141 FDSN stations (IRIS and GEOSCOPE networks) are processed to build the inverted data sets. 229 422 vertical seismograms were manually inspected and only 19 020 were picked and retained for phase velocity measurements. As explained in Beucler *et al.* (2003), all sufficiently close events are clustered into boxes of approximately $2° × 2°$, and each individual ray path belonging to the same box is considered to give equivalent results as a common ray path. Consequently, the amount of the integrated phase velocities is significantly reduced, compared to other studies. The ray path coverages corresponding to the 100 s fundamental mode of Rayleigh waves and the 150 s second higher mode are shown in Fig. 3.

The $C_d$ operators are diagonal and are only composed of *a posteriori* variances provided by the roller-coaster algorithm. Considering that our knowledge of the amplitude of isotropic and anisotropic perturbations does not differ with respect to the period or the mode branch, the same $C_p$ partitioned matrix is used for the inversions. The parameter variations are expected to be within 5 per cent and 1 per cent for the isotropic and the anisotropic components, respectively. The initial correlation lengths are set to 600 km for the fundamental mode and 1050 km and 1250 km for the first and the second overtones at 150 s, respectively.

According to the conclusion of the synthetic tests presented in the Appendix B, both $2\psi$ and $4\psi$ anisotropic parameters are taken into account in the inversions of our Rayleigh wave data sets. Since undulations in the topography of discontinuities have no more than a second-order effect on the azimuthal terms (Trampert & Woodhouse 2003), crustal corrections are not applied at this stage. For each calculation
<figure><img src="image_10.png" /></figure>

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

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

<figure><img src="image_13.png" /><figcaption>**Figure 9.** 0ψ inversion results of real Rayleigh wave data sets. Our 0ψ models are plotted on the left of TW2003's results, for the fundamental mode at (a) 60 s and (b) 100 s. Our 0ψ models for (c) the first and (d) the second overtones at 150 s are shown beside vHW1999's results. Plate boundaries are plotted in blue.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 10.** 2ψ (top) and 4ψ (bottom) models for the fundamental mode of Rayleigh waves at 60 s. Anisotropic parameters derived by the CLASH are displayed in the left column while TW2003's results are shown on the right side. The grey scale in the background corresponds to the peak-to-peak amplitude of anisotropy, expressed with respect to the PREM phase velocity. The black segments represent the fast-axis directions, which are also scaled to the amplitude shown in the background. Plate boundaries are plotted in white.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Same as in Fig. 10 but for the Rayleigh wave fundamental mode at 100 s.</figcaption></figure>
<table><caption>**Table 2.** Processing information and $\chi^2$ values for the inversions of real Rayleigh wave data sets. $\mathbf{p}_0$ and $\mathbf{p}_f$ refer to the starting and the final models, respectively.</caption><thead><tr><th></th><th>Amount of data</th><th>Max corr. length (km)</th><th>Mean corr. length (km)</th><th>χ<sup>2</sup>(p<sub>0</sub>)</th><th>χ<sup>2</sup>(p<sub>f</sub>)</th><th>Variance red. (per cent)</th><th>Figure(s) n°</th></tr></thead><tbody><tr><td>n = 0 – T = 60 s</td><td>9281</td><td>1050</td><td>623</td><td>2.9511</td><td>0.2682</td><td>90.9</td><td>9(a), 10, 12</td></tr><tr><td>n = 0 – T = 100 s</td><td>9288</td><td>1050</td><td>623</td><td>0.8516</td><td>0.1546</td><td>81.8</td><td>9(b), 11, 13</td></tr><tr><td>n = 1 – T = 150 s</td><td>1699</td><td>4150</td><td>1871</td><td>1.7588</td><td>0.8250</td><td>53.1</td><td>9(c)</td></tr><tr><td>n = 2 – T = 150 s</td><td>1831</td><td>2700</td><td>1460</td><td>1.1191</td><td>0.6063</td><td>45.8</td><td>9(d)</td></tr></tbody></table>

of phase velocity map, the starting model is the null vector, providing a good estimate of the data set variability through the $\chi^2(p_0)$ values shown in Table 2. As a final introducing remark, we used the same parametrization for all inversions. The Earth's surface is discretized in 2036 grid points (Fig. 1). Hence, the total amount of isotropic and anisotropic parameters is 10 180. Actually, since the neighbouring points are correlated to the $C_p$ operator, the amount of independent variables is smaller.

Before analysing and comparing models themselves, we can make a comment on the 1-D reference model. As other methods (e.g. van Heijst & Woodhouse 1997; Ekström et al. 1997), the roller-coaster algorithm (Beucler et al. 2003) measures phase velocities integrated along ray paths, and uses synthetic seismograms computed in PREM (Dziewonski & Anderson 1981). For a given period and a given mode branch, the discrepancy between all measured phase velocities, now considered as data, and the corresponding PREM value reflects how the PREM explains vertical seismograms. In order to evaluate this discrepancy, we compute, for each angular frequency,

$$ \frac{\delta C}{C_{\text{PREM}}} = \frac{\bar{C}(n, \omega) - C_{\text{PREM}}(n, \omega)}{C_{\text{PREM}}(n, \omega)} \times 100, \qquad (12) $$

where $\bar{C}$ is the median value of all measured phase velocities and $n$ refers to the radial order of the mode branch. The results are presented in Fig. 8, as a function of period, for the fundamental mode and the first two overtones. As pointed out by Trampert & Woodhouse (2001), we find that, for the fundamental mode, the longer the period, the better seismograms are explained by the PREM. But, with our data sets, we observe a small increase of the discrepancy for very long periods (larger than 250 s). We add that, for the first two overtones, completely different trends are observed. It might be noted that the velocity perturbations with respect to the PREM are nevertheless lower than 0.9 per cent and each curve exhibits coherent variations over the whole period range. Such behaviours may be related to the systematic exploration step of the roller-coaster algorithm, which provides a collection of starting models for further more detailed inversions. Consequently, the measured phase velocities can greatly differ from the PREM since the data signals impose it.

## 4.2 Isotropic models

The isotropic part of the constructed phase velocity maps, for the fundamental mode at 60 s and 100 s and for the first and the second overtones both at 150 s, are presented in Fig. 9. Our isotropic models are shown on the left side and are compared to the results obtained by Trampert & Woodhouse (2003, TW2003) and by van Heijst & Woodhouse (1999, vHW1999).

For both the 60 s and the 100 s fundamental mode inversions, the maximum value of laterally variable correlation lengths reaches 1050 km, at only two grid points (located around $\theta = -27^\circ$, $\phi = 218^\circ$). The necessary regularization, in order to preserve $C_p$ as a positive definite matrix, induces a average value of all correlation lengths of 623 km, which is, therefore, approximately the half the lateral resolution degree of our fundamental mode phase velocity maps. Such an observation tends to conclude that a set being composed of approximately 10 000 data, evenly distributed over the Earth's surface, is sufficient to construct phase velocity maps with a lateral resolution degree of 1000 km, say up to degree 20. For both the first and the second overtones, due to the ray path coverages, the correlation lengths are largely increased, as shown in Table 2. As a consequence, the lateral resolution degree of higher-mode phase velocity maps is degraded, mostly in the southern hemisphere. When comparing the value of the maximum correlation lengths for the first and the second higher mode in Table 2, it might be noted that only a difference of 130 ray paths can drastically increase the resolution degree.

Due to the shallow sensitivities to P and S body waves for the fundamental mode at 60 and 100 s, and kernels related to the $A_1$ parameters (Smith & Dahlen 1973; Montagner & Nataf 1986), it is well known that isotropic phase velocity perturbations are correlated to first order to the tectonic features. We can clearly observe as well this correlation in Figs 9(a) and (b). The low-velocity anomalies are mainly located in the tectonically active regions whereas high-velocity locations match more stable areas such as the cratons and the old part of oceanic pates. Concerning the stable regions it is of interest to note that the North Australian craton seems to extend off-shore into Papua New Guinea and towards the Indian Ocean, concurring with previous study results beneath Australia (Simons et al. 1999; Debayle et al. 2005). Many volcanic region locations are well correlated with low-velocity perturbations (e.g. Afar, Hawaii, Tuamotu archipelago, Azores, Iceland).

In order to compare our results, we plot TW2003's phase velocity maps on the right side of Figs 9(a) and (b). We find that our $0\psi$ models of fundamental mode at both 60 s and 100 s are very similar to TW2003's. Since the same colour scale is used for plotting maps at each period, it is essential to note that this good agreement concerns both the locations and the amplitudes of heterogeneities. We may observe however that, for both 60 s and 100 s phase velocity maps, our isotropic perturbations are more equally distributed around the corresponding PREM phase
<figure><img src="image_16.png" /><figcaption>**Figure 12.** 2ψ + 4ψ model for the fundamental mode of Rayleigh waves at 60 s. The grey scale in the background corresponds to the peak-to-peak amplitude of anisotropy, expressed with respect to the PREM phase velocity. The black segments represent the fast-axis directions, which are also scaled to the amplitude shown in the background. Plate boundaries are plotted in white.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Same as in Fig. 12 but for the Rayleigh wave fundamental mode at 100 s.</figcaption></figure>

velocity. As an example, for the 100 s isotropic model, the linear regression intercept and coefficient of our 0ψ perturbations are of 0.17 and −8.8 10⁻⁶, respectively, whereas they are of 1.08 and −2.9 10⁻⁵ for TW2003. The same feature is observed for the 0ψ perturbations at 60 s. This means that TW2003’s models are shifted of about +1 per cent with respect to the PREM phase velocity, which turns the resulting maps towards green colours. We may add that very few and small disparities can be noted, for instance in the regions surrounding Hawaii, Tuamotu archipelago or the Cap-verde islands, which raises the issue of the relationships between ridges and volcanic regions. The two different kinds
of parametrization make difficult any further quantitative comparison, which are beyond the scope of this article and will be performed in the future.

For both the first and the second overtones, our isotropic models are presented on the left side of Figs 9(c) and (d), facing the results of vHW1999. Since the sensitivities to *P* and *S* velocities for higher modes are more complicated than for the fundamental mode, the correlation between phase velocity perturbations and tectonic features is more hazardous. As pointed out by vHW1999, most first overtone high velocities are related to the cratonic areas of continents. We fairly agree with this observation except for the African continent, but it may be due to our ray path coverage (Fig. 3b). Concerning the latter point, for both the first and the second overtones, due to the laterally variable correlation lengths, a large discrepancy between northern and southern amplitude power spectra is observed. We add that a *posteriori* variances on both isotropic and anisotropic parameters (not shown here) clearly reflect uncertain measurements in the southern hemisphere (except around Australia) and for the whole African plate. Our best-resolved parts of the 0ψ model are the western part of the Pacific ocean up to central China, the Europe and the western part of North and South Americas. Due to the laterally variable correlation lengths, these regions are not affected by spurious anomalies related to Gibbs phenomenon. It should be mentioned that the first overtone low-velocity zone lying between India and the Kara sea (north of Ural mountains), exhibiting a northward stretching, is indisputably associated with large *a posteriori* variances.

When comparing our results with previous ones, we observe that, at long wavelengths, our anomaly locations in Figs 9(c) and (d) are broadly consistent with vHW1999. It appears that the heterogeneity amplitudes of phase velocity maps of vHW1999 are three times smaller than ours. This noticeable difference may be related, on one hand, to the damping factors or to the neglect of azimuthal anisotropy (as discussed in the Appendix B) for the vHW1999's maps and, on the other hand, to the relatively small amount of data for our maps. However, the best-resolved parts of our first and second overtone 0ψ models, such as the Pacific region, comprising West of Americas, Australia and East of Asia, match well those of vHW1999. We may add that both our and vHW1999's isotropic higher-mode anomalies are equally distributed around the corresponding PREM phase velocity. As example, for the first overtone, the linear regression intercept and coefficient of our 0ψ perturbations are of 0.18 and -4.9 10⁻⁵, respectively, while they are of 0.12 and -1.2 10⁻⁵ for vHW1999.

Finally, the variance reductions and the $\chi^2(p_f)$ values, shown in Table 2, tend to conclude that we can be confident in our 0$\psi$ models for the fundamental mode and it appears that the misfit between our overtone results and vHW1999 models might come from the data sets. Hence, the discrepancy between isotropic phase velocity maps is certainly more related to the different methods for measuring the higher-mode phase velocities (mode-branch stripping technique and roller-coaster algorithm) rather than to the CLASH itself.

## 4.3 Anisotropic models

When considering the ray path coverages and the variance reductions for both the first and the second overtones, we decide to only discuss the fundamental mode anisotropic models. This choice is also justified by the inability to compare higher-mode anisotropic maps with previous results. Thus, our 2$\psi$ and 4$\psi$ models for the fundamental mode at 60 s and 100 s are presented in Figs (10) and (11), respectively. They are compared with the results of TW2003. The comprehensive anisotropic perturbations, say 2$\psi$+4$\psi$, are computed as described in Section 3.1.2, and are shown in Figs 12 and 13.

As for the 0$\psi$ models, due to the sensitivities to *P* and *S* velocities for the fundamental mode, anisotropic anomalies can be related to shallow perturbations and hence to the tectonic features. Assuming a transversely isotropic medium, it is well known that the fast-axis directions of the 2$\psi$ Rayleigh terms match the directions of plate motions in the oceans (Tanimoto & Anderson 1985; Montagner & Tanimoto 1990). Qualitatively, such correlations can be observed indeed in our 2$\psi$ models (Figs 10 and 11), especially the fast-axis azimuths perpendicular to the oceanic ridges. At first order, we observe a good correlation between the 2$\psi$ fast-axis directions and the continental crust motions (e.g. India, Africa or Australia). Concerning the African continent, we observe that the 2$\psi$ fast-axis azimuths are aligned along a north-south direction, as found by Hadiouche *et al.* (1989). But, significant and rapid rotations of the 2$\psi$ azimuths are visible in the Kaapvaal craton, already observed by Freybourger *et al.* (2001). It might be noted that very low amplitudes of the 2$\psi$ terms are observed in the Afar region and in Iceland, which is coherent with the presence of large-scale reservoirs of partially melted material.

When comparing our 2$\psi$ models with TW2003's, shown on the right side of Figs 10 and 11, we globally observe the same patterns in the fast-axis directions. In both cases, our and TW2003's results, we note that the strength of the 2$\psi$ anisotropy at 60 s is greater than at 100 s, which is in agreement with the decreasing discrepancy between measured and the PREM phase velocities, as a function of period, as shown in Fig. 8. The 2$\psi$ amplitude perturbations are within the same range, say less than 1.3 per cent, but the locations of high-amplitude regions are not in good agreement. It must be, however, mentioned that TW2003's lateral resolution is of degree eight for the azimuthal part, which prevent any further quantitative comparison with our models, whose lateral resolution degree is of approximately 1200 km.

The Rayleigh 4$\psi$ perturbations are more difficult to interpret in terms of tectonic processes. As discussed in the Appendix B, they are neglected in most tomographic methods inverting Rayleigh wave data sets. However, our anisotropic models for the fundamental mode at 60 s and 100 s, shown in Figs 10 and 11, obviously display large 4$\psi$ anomalies. Even though the 2$\psi$ perturbations seem to be of larger extent than the 4$\psi$, amplitudes are of the same order of magnitude. Quantitatively, the amplitudes of the 4$\psi$ perturbations are 5 per cent smaller than the 2$\psi$ ones at 60 s, and they are 12 per cent smaller at 100 s. The reasons for the presence of Rayleigh 4$\psi$ anisotropy is not unique. It can either reflect oceanic complex mineralogical aggregates, comprising ortho- and/or clinopyroxenes, or different 2$\psi$ fossil fast-axis directions that may be encountered in regions such as old parts of oceanic plates or continents. In Figs 10 and 11, we observe indeed that the 4$\psi$ amplitudes are relatively large in many continental regions (e.g. Turkey, Tibet, central Australia, Greenland, North America, Antarctica) as well as in tectonically active areas (e.g. west coast of North America, Aleutians, west of the Pacific plate). It may be of
interest to add that all anisotropic models, between 40 s and 200 s, exhibit large 4ψ amplitudes stretched along the Hawaii hotspot track. Since the density and the azimuthal distribution of ray paths are checked during the computation of C<sub>p</sub>, it cannot exist biases due to the ray path coverages. Moreover, the fast-axis directions of these 4ψ anomalies are very coherent over the whole period range and hardly no 2ψ perturbations are visible in this region. Both petrological and thermal insights of such observations need to be further investigated. For all anisotropic models, computed at various periods, most locations of high 2ψ anisotropy do not match those of high 4ψ and the 2ψ fast-axis directions disagree with the 4ψ directions. Hence, we think that a significant part of the data signal cannot be only explained by the 2ψ terms. When considering 2ψ + 4ψ perturbations, represented in Figs 12 and 13, it appears nevertheless that the fast-axis directions are mostly dominated by the 2ψ anisotropy but some rapid rotations in specific locations show that the strength of the 4ψ anisotropy is definitely not negligible.

As previously mentioned, the difference of lateral resolution degree makes difficult any quantitative comparison with the 4ψ models of TW2003. We observe that the directions of fast-axes, at both 60 s and 100 s, are however in fairly good agreement (e.g. South of Alaska, Tonga-Kermadec zone, Tibet). The major difference between our 4ψ models and TW2003's concerns the evolution of amplitudes, as a function of period. In our case, as for the 2ψ models, the longer the period, the smaller our 4ψ anomaly amplitudes are, which confirms the trend observed in Fig. 8. Whereas for TW2003's models no significant variation is noted between the models at 60 s and 100 s. Finally, we observe that the 2ψ models at 60 s and 100 s are strongly correlated, as well as the 4ψ perturbations, which seems to indicate a strong azimuthal anisotropy in the shallower part of the mantle, as observed by TW2003.

# 5 DISCUSSION AND CONCLUDING REMARKS

We investigated the problem of the translation of horizontally propagating seismic wave velocities, integrated between the source and the receiver, into global isotropic and anisotropic perturbation velocity maps. The forward problem relies on the expression derived by Smith & Dahlen (1973), which is an extension of the Backus' decomposition applied to the surface wave case. As many other tomographic methods, models are calculated by using an iterative gradient least-squares algorithm. Ray paths are assumed to follow great circles, joining epicentres to stations, and we are aware that a more realistic propagation for surface waves has to be accounted for in the future by modifying the path sensitivity kernels. The relationships between grid points and data only depend on the crossing point locations between all ray paths, which prevents artefacts due to any preferential direction of ray bundles, usually met when dealing with real data sets. The inescapable uneven sampling of the medium, due to the locations of events and stations, is taken into account through the *a priori* covariance operator on parameters with laterally variable correlation lengths. The CLASH, is able to handle massive data sets in a reasonable computing time, with no simplification during matrix computations.

As mentioned in the Appendix B, several existing methods do not consider azimuthal anisotropy or only take into account the 2ψ or the 4ψ terms, based on synthetic experiment results and also on the theoretical dependence with depth between anisotropic coefficients and elastic parameters. In order to compare the relative weight of each anisotropic term on the total perturbation, Montagner & Nataf (1986) used realistic C<sub>ij</sub> matrices of elastic coefficients to simulate the azimuthal anisotropy in the lithosphere and in the asthenosphere. They found that the phase velocities of the Rayleigh fundamental mode have a significant 2ψ anisotropy and a negligible 4ψ whereas the *a contrario* trend is observed for the Love wave case. When they consider a realistic C<sub>ij</sub> matrix corresponding to a hartzburgite (Peselnick & Nicolas 1978), the same 2ψ dominant feature is found for the Rayleigh waves but concerning the Love wave anisotropy, they conclude that the amplitude of both 2ψ and 4ψ terms is quite small and somewhat dependent upon the C<sub>ij</sub> matrix. Other studies show that the effect of the azimuthal anisotropy is, in general, relatively small for Love waves (Schlue 1977; Kawakasi 1986; Montagner & Nataf 1986). And, Debayle *et al.* (1998) also showed that a negligible 4ψ Rayleigh and 2ψ Love wave anisotropy is regularly observed in fully, but weakly, anisotropic models. More recently, Trampert & van Heijst (2002) found that Love wave overtones are strongly sensitive to the anisotropic 2ψ terms.

On the other hand, Laske & Masters (1998) found that, although the contribution of the 4ψ terms is expected to be small for Rayleigh waves, this result is, however, not observed with their inverted data sets. They add that, if the 4ψ terms are included in the inversion, their spherical harmonic coefficients are almost as large as those for the 2ψ terms and the only restriction to retrieve reliable 4ψ terms comes from the ray path coverage. The recent results concerning the fundamental mode of Love and Rayleigh waves (Trampert & Woodhouse 2003) showed that Love waves are mostly sensitive to the 4ψ anisotropy indeed whereas Rayleigh wave data sets need both 2ψ and 4ψ terms.

We address here the question of the sensitivity of the anisotropic terms from two different points of view. We investigate, using synthetic data sets, the existence of biases when some of the anisotropic terms are neglected (Appendix B). It is clearly shown that no bias can be observed when azimuthal anisotropy is taken into account in the inversion process in its comprehensive form, say 2ψ and 4ψ terms. This observation holds even when perturbations to be recovered are only isotropic. Conversely, the neglect of one of the two kinds of anisotropy leads to strong artefacts in the resulting models.

According to this observation, we secondly present isotropic and both 2ψ and 4ψ anisotropic models, derived by inversions of real Rayleigh wave data sets. We observe that amplitudes of the 2ψ and the 4ψ terms are of the same order of magnitude. At various periods, a significant 4ψ signature is located along the Hawaii hotspot track and needs further investigations. Another interesting feature concerns the azimuthal anisotropic perturbations located over the continental lithospheres. For instance, at both 60 s and 100 s, large 4ψ terms are observed in the central and eastern parts of Australia (Figs 10 and 11). In this region, the 2ψ fast-axis directions agree well with the hotspot reference frame but, due to the considerable 4ψ heterogeneity, the comprehensive anisotropic fast-axis directions, shown in Figs 12 and 13, exhibit rapid and significant variations. Recent tomographic studies of the Australian lithosphere (Debayle & Kennett 2000; Simons *et al.* 2002; Simons &
van der Hilst 2003; Debayle *et al.* 2005) raised the problem of possible frozen anisotropy in the shallow lithosphere and rapid rotation of the fast-axis directions around 200 km depth. However, all shear-wave velocity models are inferred neglecting the 4ψ anisotropy and it appears that a certain part of the apparent discrepancy might be associated with the presence of large 4ψ anisotropic seismic heterogeneities.

Finally, concurring with the conclusions of Trampert & Woodhouse (2003), we think that Rayleigh wave data sets do need both 2ψ and 4ψ anisotropic perturbations in order to be correctly explained.

## ACKNOWLEDGMENTS

The CLASH acronym is given in memory of Joe Strummer who passed away in December 2002. Discussions with Albert Tarantola were really fruitful. Many thanks to Hendrik van Heijst and Jeannot Trampert for authorizing us to use their models. We would like to thank Jeannot Trampert and two anonymous reviewers, and Rob van der Hilst as editor, for very constructive reviews and critics. All figures are made by using the Generic Mapping Tools (Wessel & Smith 1998) but Fig. 2, realized with the freeware xfig. Data are provided by the IRIS and GEOSCOPE networks and we would like to warmly thank the corresponding staff for making such high-quality seismograms available. This work was partly supported by the program INSU-Dyeti.

## REFERENCES

* Backus, G.E., 1965. Possible forms of seismic anisotropy of the uppermost mantle under oceans, *J. geophys. Res.*, **70**, 3429–3439.

* Backus, G.E., 1970. A geometrical picture of anisotropic elastic tensors, *Rev. Geophys. Space Phys.*, **8**, 633–671.

* Barmin, M.P., Ritzwoller, M.H. & Levshin, A.L., 2001. A fast and reliable method for surface wave tomography, *Pure appl. geophys.*, **158**, 1351–1375.

* Beucler, É., Stutzmann, É. & Montagner, J.-P., 2003. Surface-wave higher mode phase velocity measurements using a roller coaster type algorithm, *Geophys. J. Int.*, **155**, 289–307.

* Debayle, É. & Kennett, B.L.N., 2000. Anisotropy in the Australasian upper mantle from love and Rayleigh waveform inversion, *Earth planet. Sci. Lett.*, **184**, 339–351.

* Debayle, É. & Sambridge, M., 2004. Inversion of massive surface wave data sets: model construction and resolution assessments, *J. geophys. Res.*, **109**, B02316, doi:10.1029/2003JB002652.

* Debayle, É., Lévêque, J.-J. & Maupin, V., 1998. Anisotropy in the Indian Ocean upper mantle from Rayleigh- and love-waveform inversion, *Geo-phys. J. Int.*, **133**, 529–540.

* Debayle, É., Kennett, B.L.N. & Priestley, K., 2005. Global azimuthal seismic anisotropy and the unique plate-motion deformation of Australia, *Nature*, **433**, 509–512, doi:10.1038/nature03247.

* Dziewonski, A.M., 1984. Mapping the lower mantle: Determination of lateral heterogeneity in P velocity up to degree and order 6, *J. geophys. Res.*, **89**, 5929–5952.

* Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference earth model, *Phys. of Earth and Planet. Int.*, **25**, 297–356.

* Ekström, G. & Dziewonski, A.M., 1998. The unique anisotropy of the Pacific upper mantle, *Nature*, **394**, 168–172.

* Ekström, G., Tromp, J. & Larson, E.W.F., 1997. Measurements and global models of surface wave propagation, *J. geophys. Res.*, **102**, 8137–8157.

* Fisher, R.A., 1925. Theory of statistical estimation, *Proc. Camb. Phil. Soc.*, **22**, 700–725.

* Fisher, R.A., 1953. Dispersion on a sphere, *Proc. Roy. Soc., A* **217**, 295.

* Freybourger, M., Gaherty, J.B., Jordan, T.H. & the Kaapvaal Seismic Group, 2001. Structure of the kaapvaal craton from surface waves, *Geophys. Res. Lett.*, **28**, 2489–2492.

* Grand, S.P., van der Hilst, R.D. & Widiyantoro, S., 1997. Global seismic tomography: a snapshot of convection in the earth, *GSA Today*, **7**, 1–7.

* Hadiouche, O., Jobert, N. & Montagner, J.-P., 1989. Anisotropy of the African continent inferred from surface waves, *Phys. of Earth and Planet. Int.*, **58**, 61–81.

* Kawakasi, I., 1986. Azimuthally anisotropic model of the oceanic upper mantle, *Phys. of Earth and Planet. Int.*, **43**, 1–21.

* Laske, G., 1995. Global observation of off-great-circle propagation of long period surface waves, *Geophys. J. Int.*, **123**, 245–259.

* Laske, G. & Masters, G., 1996. Constraints on global phase velocity maps from long-period polarization data, *J. geophys. Res.*, **101**, 16 059–19 075.

* Laske, G. & Masters, G., 1998. Surface-wave polarization data and global anisotropic structure, *Geophys. J. Int.*, **132**, 508–520.

* Levshin, A.L., Barmin, M.P., Ritzwoller, M.H. & Trampert, J., 2005. Minor-arc and major-arc global surface wave diffraction tomography, *Phys. of Earth and Planet. Int.*, **149**, 205–223.

* Li, X.D. & Romanowicz, B., 1995. Comparison of global waveform inversion with and without cross-branch coupling, *Geophys. J. Int.*, **121**, 695–709.

* Mégnin, C. & Romanowicz, B., 2000. The three-dimensional shear velocity structure of the mantle from the inversion of body, surface and higher-mode waveforms, *Geophys. J. Int.*, **143**, 709–728.

* Montagner, J.-P., 1986. Regional three-dimensional structures using long-period surface waves, *Annales Geophysicae*, **4**, 283–294.

* Montagner, J.-P. & Nataf, H.C., 1986. A simple method for inverting the azimuthal anisotropy of surface waves, *J. geophys. Res.*, **91**, 511–520.

* Montagner, J.-P. & Tanimoto, T., 1990. Global anisotropy in the upper mantle inferred from the regionalization of phase velocities, *J. geophys. Res.*, **95**, 4797–4819.

* Montagner, J.-P. & Tanimoto, T., 1991. Global upper mantle tomography of seismic velocities and anisotropies, *J. geophys. Res.*, **96**, 20 337–20 351.

* Nataf, H.C., Nakanishi, I. & Anderson, D.L., 1986. Measurements of mantle velocities and inversion for lateral heterogeneities and anisotropy III, *J. geophys. Res.*, **91**, 7261–7307.

* Nolet, G. & Montelli, R., 2005. Optimal parametrization of tomographic models, *Geophys. J. Int.*, **161**, 365–372.

* Okal, E., 1978. A physical classification of the Earth's spheroidal modes, *J. Phys. Earth*, **26**, 75–103.

* Peselnick, L. & Nicolas, A., 1978. Seismic anisotropy in an ophiolite peridotite, *J. geophys. Res.*, **83**, 1227–1235.

* Ritsema, J. & van Heijst, H.J., 2000. Seismic imaging of structural heterogeneity in earth's mantle: evidence for large-scale mantle flow, *Science Progress*, **83**, 243–259.

* Ritsema, J., van Heijst, H.J. & Woodhouse, J.H., 1999. Complex shear wave velocity structure imaged beneath Africa and Iceland, *Science*, **286**, 1925–1928.

* Sambridge, M., 1999. Geophysical inversion with a neighbourhood algorithm-i. searching a parameter space, *Geophys. J. Int.*, **138**, 479–494.

* Schlue, J.W., 1977. A physical model for surface wave azimuthal anisotropy, *Bull. Seismol. Soc. Am.*, **67**, 1515–1519.

* Simons, F.J. & van der Hilst, R.D., 2003. Seismic and mechanical anisotropy and the past and present deformation of the Australian lithosphere, *Earth planet. Sci. Lett.*, **211**(3-4), 271–286.

* Simons, F.J., Zielhuis, A. & van der Hilst, R.D., 1999. The deep structure of the Australian continent from surface wave tomography, *Lithos*, **48**, 17–43.
Simons, F.J., van der Hilst, R.D., Montagner, J.-P. & Zielhuis, A., 2002. Multi-mode Rayleigh wave inversion for heterogeneity and azimuthal anisotropy of the Australian upper mantle, *Geophys. J. Int.*, **151**, 738–754.

Smith, M.L. & Dahlen, F.A., 1973. The azimuthal dependence of love and Rayleigh wave propagation in a slightly anisotropic medium, *J. geophys. Res.*, **78**, 3321–3333.

Smith, M.L. & Dahlen, F.A., 1975. Correction to ‘The azimuthal dependence of Love and Rayleigh wave propagation in a slightly anisotropic medium’, *J. geophys. Res.*, **78**, 1923.

Spetzler, J. & Trampert, J., 2003. Implementing spectral leakage corrections in global surface wave tomography, *Geophys. J. Int.*, **155**, 532–538.

Spetzler, J., Trampert, J. & Snieder, R., 2001. Are we exceeding the limits of the great circle approximation in global surface wave tomography?, *Geophys. Res. Lett.*, **28**, 2341–2344.

Su, W.J., Woodward, R.L. & Dziewonski, A.M., 1994. Degree 12 model of shear velocity heterogeneity in the mantle, *J. geophys. Res.*, **99**, 6945–6980.

Tanimoto, T., 1990. Long wavelength s-velocity structure throughout the mantle, *Geophys. J. Int.*, **100**, 327–336.

Tanimoto, T. & Anderson, D.L., 1985. Lateral heterogeneity and azimuthal anisotropy of the upper mantle: Love and Rayleigh waves 100–250 sec., *J. geophys. Res.*, **90**, 1842–1858.

Tarantola, A., 1987. *Inverse problem theory*, Elsevier.

Tarantola, A. & Valette, B., 1982. Inverse problems=quest for information, *J. Geophys.*, **50**, 159–170.

Tarantola, A. & Valette, B., 1982. Generalized nonlinear inverse problems solved using the least squares criterion, *Rev. Geophys. Space Phys.*, **20**, 219–232.

Trampert, J. & Snieder, R., 1996. Model estimation biased by truncated expansions—possible artifacts in seismic tomography, *Science*, **271**, 1257–1260.

Trampert, J. & van Heijst, H.J., 2002. Global azimuthal anisotropy in the transition zone, *Science*, **296**, 1297–1299.

Trampert, J. & Woodhouse, J.H., 1995. Global phase velocity maps of love and Rayleigh waves between 40 and 150 seconds, *Geophys. J. Int.*, **122**,

675–690.

Trampert, J. & Woodhouse, J.H., 1996. High-resolution global phase velocity distribution, *Geophys. Res. Lett.*, **23**, 21–24.

Trampert, J. & Woodhouse, J.H., 2001. Assessment of global phase velocity models, *Geophys. J. Int.*, **144**, 165–174.

Trampert, J. & Woodhouse, J.H., 2003. Global anisotropic phase velocity maps for fundamental mode surface waves between 40 and 150 s, *Geophys. J. Int.*, **154**, 154–165.

van der Hilst, R.D., Widyantoro, S. & Engdahl, E.R., 1997. Evidence for deep mantle circulation from global tomography, *Nature*, **386**, 578–584.

van Heijst, H.J. & Woodhouse, J.H., 1997. Measuring surface-wave over-tone phase velocities using a mode-branch stripping technique, *Geophys. J. Int.*, **131**, 209–230.

van Heijst, H.J. & Woodhouse, J.H., 1999. Global high-resolution phase velocity distributions of overtone and fundamental-mode surface waves determined by mode-branch stripping, *Geophys. J. Int.*, **137**, 601–620.

Wang, Z. & Dahlen, F.A., 1995. Spherical-spline parametrization of three-dimensional earth, *Geophys. Res. Lett.*, **22**, 3099–3102.

Wang, Z., Tromp, J. & Ekström, G., 1998. Global and regional surface-wave inversions: A spherical-spline parametrization, *Geophys. Res. Lett.*, **25**, 207–210.

Wessel, P. & Smith, W.H.F., 1998. New, improved version of the generic mapping tools released, *EOS, Trans. Am. geophys. Un.*, **79**, 579.

Wong, Y.K., 1989. Upper mantle heterogeneity from phase and amplitude data of mantle waves, *Ph.D. thesis*, Harvard University, Cambridge, Massachusetts.

Woodhouse, J.H. & Dziewonski, A.M., 1984. Mapping the upper mantle: three-dimensional modeling of earth structure by inversion of seismic waveforms, *J. geophys. Res.*, **89**, 5953–5986.

Woodhouse, J.H. & Dziewonski, A.M., 1986. Three-dimensional mantle models based on mantle wave and long period body wave data, *EOS, Trans. Am. geophys. Un.*, **67**, 307.

Woodhouse, J.H. & Wong, Y.K., 1986. Amplitude, phase and path anomalies of mantle waves, *Geophys. J. R. astr. Soc.*, **87**, 753–773.

## APPENDIX A: QUANTITATIVE COMPARISONS IN THE MODEL SPACE

For all synthetic test results, in order to make some quantitative comparisons of both isotropic and anisotropic parameters, distances between input and output velocities are shown. For a given point, located on the Earth's surface, distance functions (equivalent to normalized differences for linear variables), denoted as $δ_l$, are defined as

$$ \forall l \in \{1, 2, 3, 4, 5\}, \forall i \in [1, \text{length}(\mathbf{pl})], \quad \delta_l(\theta_i, \phi_i) = \left\| \ln \left( \frac{C_l(\theta_i, \phi_i)^{\text{in}}}{C_l(\theta_i, \phi_i)^{\text{out}}} \right) \right\| \times 100, \qquad (A1) $$

where input and output velocities are noted as $C_l(\theta_i, \phi_i)^{\text{in}}$ and $C_l(\theta_i, \phi_i)^{\text{out}}$, respectively. They are defined as a function of the parameter subvector components,

$$ C_l(\theta_i, \phi_i) = C_0 (1 + \mathbf{pl}_i). \qquad (A2) $$

The subscript $l$ refers to the subvector type: $l=1$ for the 0$\psi$ parameters; $l=2$ and $3$ for the 2$\psi$ terms (cosine and sine, respectively) and $l=4$ and $5$ for the 4$\psi$ terms (cosine and sine, respectively). The quantitative comparison maps are plotted using Hammer projection.

## APPENDIX B: SYNTHETIC INVESTIGATIONS OF THE AZIMUTHAL ANISOTROPY

Considering that previously shown synthetic tests emphasize the reliability of the CLASH, we address the question to know whether our method can efficiently separate the 2$\psi$ from the 4$\psi$ anisotropic effects and whether this could have some consequences on the phase velocity models. As shown in eq. (6), the first-order effects of a small anisotropy on the propagation of surface waves are expressed as Fourier series in $\psi$, involving both 2$\psi$ and 4$\psi$ parameters. Assuming a transversely isotropic medium, several studies showed that the Rayleigh waves are only slightly sensitive to the 4$\psi$ terms and conversely for the Love waves with the 2$\psi$ terms (Smith & Dahlen 1973; Montagner & Nataf 1986; Debayle et al. 1998); and moreover conclude that most data signal can be explained by isotropic perturbations. Consequently some global tomographic methods neglect the 4$\psi$ anisotropy for Rayleigh waves, and conversely for Love waves, while others do not account for anisotropy at all. On the other hand, Trampert & Woodhouse (2003) introduced both 2$\psi$ and 4$\psi$ parameters for their inversions of Love and Rayleigh wave data sets, hence neglecting none of the anisotropic effects. They conclude that the Love wave fundamental mode data do
<figure><img src="image_18.png" /><figcaption>**Figure B1.** Test #2 quantitative comparisons between input and output models. For each subvector, referred to as *l*, the distances between input and output velocities, $\delta_l$, are computed as defined in eq. (A1). Contour lines are plotted each 0.5 distance unit.</figcaption></figure>

not require 2$\psi$ anisotropy indeed, whereas the Rayleigh waves need both 2$\psi$ and 4$\psi$ terms. Therefore, the question to know whether some azimuthal anisotropy form can be a *priori* neglected in the inversions is worth investigating and is attempted by using synthetic data sets. For each test presented in this appendix, the ray path coverage corresponds to one shown in Fig. 3(a).

For the Tests #3 and #4, the data set is synthetized in the same input model as used for the Test #2, including large isotropic and both 2$\psi$ and 4$\psi$ anisotropic perturbations (Table 1). Neglecting the 4$\psi$ anisotropy, the Test #3 inversion accounts for the 0$\psi$ and only the 2$\psi$ parameters. For the Test #4, only the 0$\psi$ parameters are introduced in the inversion process, which means that the whole synthetic data set has to be explained by isotropic perturbations. The quantitative differences, between input and output velocities, are presented in Figs B2(a) and B2(b).

As expected, compared to the Test #2 results, shown in Fig. B1, the agreement between input and output models is largely degraded. The misfit can be estimated in the data space through the $\chi^2(p_f)$ values or the variance reductions (Table 1). Concerning the Test #3, a large part of the 4$\psi$ input signal—which cannot be explained, since this kind of anisotropy is not introduced in the inversion—is transferred preferentially into the 0$\psi$ perturbations rather than into the 2$\psi$ parameters. The largest distance between input and output velocities (eq. A1) reaches the value of 2.98 for the 0$\psi$ parameters (located in central Pacific ocean) and they are of 1.39 and 1.17 for the cos 2$\psi$ and the sin 2$\psi$ terms, respectively. As predictable, the worst-case scenario corresponds to the Test #4, where the whole anisotropy of the input model is expected to be explained only by the isotropic anomalies. As shown in Fig. B2(b), the resulting model is strongly biased and the differences between input and output velocities reach the value of 3.97.

For the Tests #5 and #6, the data set is synthetized in the same input model as used for the Test #2, where now, for the purpose of the experiment, the 4$\psi$ anisotropic components are set to zero (Table 1). Hence, the model to be recovered comprises 0$\psi$ and only 2$\psi$ perturbations, which consequently decreases, in the data space, the $\chi^2(p_0)$ values from 4.5269 to 4.3117, with respect to the Test #2. Since the model resulting of the inversion accounting for 0$\psi$ and only 2$\psi$ parameters is almost identical to the input model, this test is not detailed
<figure><img src="image_19.png" /><figcaption>**Figure B2.** Same as in Fig. B1 but for the Tests #3 and #4. The input model contains 0ψ and both 2ψ and 4ψ perturbations. For the Test #3, (a) 0ψ and only 2ψ parameters are taken into account during inversion whereas for the Test #4, (b) only the 0ψ perturbations are introduced.</figcaption></figure>

here. Following the same strategy as for the Tests #3 and #4, the Test #5 inversion is performed with the 0ψ parameters and both 2ψ and 4ψ terms, whereas for the Test #6, only the 0ψ parameters are introduced. The distances between input and output velocities are presented in Figs B3(a) and B3(b).

Similarly to the Test #4, neglecting the anisotropic components with respect to the input model, as done for the Test #6, produces strong artefacts. As shown in Fig. B3(b), the velocity differences reach 2.93, although locations of heterogeneities of the derived model (not shown here) are still rather satisfactorily retrieved. The misfit can be also estimated through the $\chi^2(p_f)$ value, but it might be noted that the variance reduction seems nevertheless to imply good inversion results.

The most interesting feature, revealed by these synthetic tests, concerns the Test #5 output model. We remind that for this experiment, the three kinds of parameters, say 0ψ and both 2ψ and 4ψ, are taken into account during inversion, but only 0ψ and 2ψ parameters have to be recovered. The resulting model (not shown here) is almost identical to the input one, and the values of the output 4ψ parameters are accordingly negligible. As shown in Table 1, the $\chi^2(p_f)$ value reflects the very good agreement between the input and the output models. The quantitative differences for the 0ψ, the cos 2ψ and the sin 2ψ parameters, that is, for the **p**1, **p**2 and **p**3 subvectors, are presented in Fig. B3(a). The largest distances between input and output velocities reach 1.11 for the 0ψ terms while they reach values of 0.58 and 0.62 for the cos 2ψ and the sin 2ψ terms, respectively.

Hence, when the 4ψ anisotropy is not present in the data signal, the introduction of both 2ψ and 4ψ parameters does not affect the inversion results and consequently no 4ψ anisotropy is retrieved. We already showed with Test #1 that when heterogeneities to be recovered are purely isotropic, no output anisotropy is observed; the Test #5 result tends to show that the CLASH can efficiently separate the two kinds of anisotropy.

For the purpose of the Tests #7 and #8, the data set is now synthetized in the isotropic part of the input model used for the Test #2 (Table 1). Hence, the model to be recovered is purely isotropic, which consequently decreases the $\chi^2(p_0)$ values from 4.5269 to 4.2026, with respect to the Test #2. The Test #7 inversion is performed with the 0ψ and both the 2ψ and the 4ψ parameters whereas for the Test #8, the 0ψ and only the 2ψ parameters are introduced. The distances between input and output velocities are presented in Figs B4(a) and B4(b).
<figure><img src="image_20.png" /></figure>

**Figure B3.** Same as in Fig. B1 but for the Tests #5 and #6. The input model contains $0\psi$ and only $2\psi$ perturbations. For the Test #5, (a) $0\psi$ and both $2\psi$ and $4\psi$ parameters are taken into account during inversion whereas for the Test #6, (b) only the $0\psi$ perturbations are introduced.

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

**Figure B4.** Same as in Fig. B1 but for the Tests #7 and #8. The input model contains only $0\psi$ perturbations. For the Test #7, (a) $0\psi$ and both $2\psi$ and $4\psi$ parameters are taken into account during inversion whereas for the Test #8 (b) $0\psi$ and only $2\psi$ perturbations are introduced.
Strengthening the feature revealed by the Test #5, for both the Tests #7 and #8 the isotropic part of output models is almost identical to the input model. Consequently, no azimuthal anisotropy is recovered. The almost perfect matching between input and output models can be estimated through the $\chi^2(p_f)$ values. As shown in Figs B4(a) and B4(b), the quantitative differences between input and output velocities do not exceed the value of 0.9, say less than 1 per cent, for the both tests.

To summarize, we observe that when the 0ψ and both 2ψ and 4ψ parameters are taken into account during inversion—whatever the kind of anisotropy to be recovered—the CLASH retrieves very correctly the input models. This observation holds when data set is synthetized in a purely isotropic input model and, therefore, tends to conclude that an *a priori* neglect of at least one kind of the azimuthal anisotropy leads to artefacts.