# Surface wave tomography for azimuthal anisotropy in a strongly reduced parameter space

M.P. Panning, G. Nolet

► To cite this version:

M.P. Panning, G. Nolet. Surface wave tomography for azimuthal anisotropy in a strongly reduced parameter space. Geophysical Journal International, 2008, 174 (2), pp.629-648. 10.1111/j.1365-246X.2008.03833.x. hal-00407882

HAL Id: hal-00407882

https://hal.science/hal-00407882v1

Submitted on 2 Oct 2021

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Surface wave tomography for azimuthal anisotropy in a strongly reduced parameter space

Mark P. Panning<fnref n="1" /> and Guust Nolet<fnref n="1" />,<fnref n="2" />

Accepted 2008 April 25. Received 2008 April 24; in original form 2008 February 11

## SUMMARY

Large scale seismic anisotropy in the Earth's mantle is likely dynamically supported by the mantle's deformation; therefore, tomographic imaging of 3-D anisotropic mantle seismic velocity structure is an important tool to understand the dynamics of the mantle. While many previous studies have focused on special cases of symmetry of the elastic properties, it would be desirable for evaluation of dynamic models to allow more general axis orientation. In this study, we derive 3-D finite-frequency surface wave sensitivity kernels based on the Born approximation using a general expression for a hexagonal medium with an arbitrarily oriented symmetry axis. This results in kernels for two isotropic elastic coefficients, three coefficients that define the strength of anisotropy, and two angles that define the symmetry axis. The particular parametrization is chosen to allow for a physically meaningful method for reducing the number of parameters considered in an inversion, while allowing for straightforward integration with existing approaches for modelling body wave splitting intensity measurements. Example kernels calculated with this method reveal physical interpretations of how surface waveforms are affected by 3-D velocity perturbations, while also demonstrating the non-linearity of the problem as a function of symmetry axis orientation. The expressions are numerically validated using the spectral element method. While challenges remain in determining the best inversion scheme to appropriately handle the non-linearity, the approach derived here has great promise in allowing large scale models with resolution of both the strength and orientation of anisotropy.

**Key words:** Inverse theory; Surface waves and free oscillations; Seismic anisotropy; Seismic tomography; Theoretical seismology; Wave scattering and diffraction.

## 1 INTRODUCTION

Global tomographic modelling of the 3-D seismic velocity structure of the mantle is an important tool to allow us to image the dynamic processes of the Earth. Most tomographic models to date have modelled the isotropic *P* and/or *S* velocity. Under the assumption of relative chemical heterogeneity, these can be well interpreted as a snapshot of the current thermal state of the mantle. These thermal differences can then be interpreted as the driving force through thermal expansion for the dynamic processes of the mantle. The good agreement of isotropic structure at long wavelengths of many recent *S* velocity models (Grand 1997; Masters et al. 2000; Mégnin & Romanowicz 2000; Ritsema & van Heijst 2000; Gu et al. 2003) is a constraint on the current long-wavelength thermal structure in the mantle, while some finer scale isotropic *P* velocity models have imaged both cold, descending slabs (e.g. van der Hilst et al. 1997; Li et al. 2006) and hot, upwelling plumes (Montelli et al. 2004, 2006).

Such isotropic modelling has only an indirect connection to dynamics through the assumption that velocity heterogeneity can be explained primarily by thermal heterogeneity and can, therefore, be related to the flow-driving density contrasts through thermal expansion. Considering the anisotropy of seismic velocity, however, may allow for a more direct connection to the dynamic flow of the mantle. On a microscopic scale, most minerals that make up the Earth's mantle are elastically anisotropic, but random orientations of these crystals leads to isotropy on the scale of seismic wavelengths. The dynamic deformation processes of the mantle can produce seismically observable anisotropy either through preferentially aligning the crystalline axes (lattice preferred orientation or LPO) (e.g. Karato 1998) or alignment of materials with strongly contrasting elastic properties (shape preferred orientation or SPO) (e.g. Kendall & Silver 1996). In the relatively cold lithosphere, such signatures may be frozen in over geological timescales (Silver 1996), but they likely require dynamic support at greater depths (Vinnik et al. 1992).

This requirement of dynamic support suggests that observed seismic anisotropy, if sufficiently resolved, has the potential to serve as a direct proxy to the dynamic mantle flow field, although we require further information to understand how to relate flow to seismic anisotropy. Theoretical studies of predicted anisotropy due to mantle flow based on kinematic theory (e.g. Becker et al. 2003, 2006, 2008) suggest that at

<footnote n="1">¹Department of Geosciences, Princeton University, Princeton, NJ 08544, USA. E-mail: mpanning@princeton.edu</footnote>

<footnote n="2">²Geosciences Azur, Université de Nice/Sophia Antipolis, Observatoire de la Côte d'Azur, Université Pierre et Marie Curie, IRD, CNRS, France</footnote>
least upper-mantle anisotropy may be well modelled using an anisotropic material with hexagonal symmetry, with the fast axis oriented roughly parallel with the long axis of the finite strain ellipse. Detailed anisotropic modelling could then be compared with mantle flow calculations in order to better constrain the dynamics of the mantle. It is valuable then to attempt to directly resolve both strength and orientation of anisotropy.

There are already several available models of seismic anisotropy. Models of seismic anisotropy based on surface waves, though, often only deal with anisotropy with a vertical axis of symmetry (radial anisotropy) (e.g. Ekström & Dziewonski 1998; Beghein & Trampert 2004; Panning & Romanowicz 2006; Kustowski et al. 2008). Conversely, body wave studies such as SKS splitting measurements have typically been modelled assuming a horizontal axis of symmetry for the material properties (e.g. Vinnik et al. 1992; Fouch & Fischer 1996; Silver 1996). Recently, however, methods have been developed for using splitting measurements in tomographic inversions for more general symmetry axis orientations (Chevrot 2006; Abt & Fischer 2008; Long et al. 2008). There have also been surface wave studies which model azimuthal dependence, in particular, of fundamental mode Rayleigh wave velocity through linearized dependence on a horizontal fast axis (e.g. Montagner & Tanimoto 1991; Simons et al. 2002; Debayle et al. 2005).

In general, allowing for arbitrary symmetry axis orientation and anisotropic strength would lead to 21 independent elastic coefficients, which is a very large number for use in an inverse approach. In order to best use surface wave measurements to constrain the strength and anisotropy through a stable inversion process, we require a method that allows for a reasonable way to reduce the number of parameters for which we invert. An important method for modelling surface wave sensitivity to anisotropic structure through 13 ‘natural’ linear combinations of the elastic constants (Montagner & Nataf 1986; Montagner & Jobert 1988) has been proposed, and arguments based on the relative sensitivity to these different parameters and the expected amplitudes of these parameters in the Earth have led to inversions for a smaller subset of these parameters (e.g. Montagner & Tanimoto 1991; Simons et al. 2002). However, recent numerical studies (Sieminski et al. 2007a) suggest that Rayleigh waves may exhibit significant sensitivity to a large number of these parameters. Additionally, if combination with body wave measurements is desired, such phases exhibit sensitivity to a different subset of linear combinations of elastic parameters (Sieminski et al. 2007b). In order to have a physically meaningful method to reduce the number of parameters, it is advantageous to consider a parametrization explicitly based upon an arbitrarily oriented hexagonal medium (Chevrot 2006), which has five elastic parameters plus two orientation angles to describe the symmetry axis. Additional reduction in parameters can then be obtained as both theoretical studies of material deformation due to mantle convection (Becker et al. 2006) and observations of mantle xenoliths (Montagner & Anderson 1989) suggest significant correlation between different elastic parameters.

In order to maximize the resolution of an anisotropic model, we would also like to combine surface wave and shear wave splitting measurements (e.g. Marone & Romanowicz 2007). Such data sets are complementary in the sense that the surface wave measurements have the potential for excellent depth resolution, while the shear wave splitting measurements allow for greater lateral resolution and determination of the orientation of the fast axis in the horizontal plane. Since such measurements are made in very different frequency bands, though, combination of these data sets could be greatly aided by the use of kernels calculated using finite-frequency theory (e.g. Dahlen et al. 2000; Zhao et al. 2000; Zhou et al. 2004).

In this paper, we derive 3-D finite-frequency kernels for surface waveforms based on the Born approximation with a model parametrized by an arbitrarily oriented hexagonal symmetry (Chevrot 2006). Such an approach allows for a physically meaningful way to reduce the number of anisotropic model parameters, as well as providing a complementary approach to invert surface wave data in combination with SKS splitting intensity measurements (Chevrot 2006; Long et al. 2008) in a consistent theoretical framework. Such an approach is, however, non-linear and we suggest some possible approaches to deal with such non-linearity.

## 2 SURFACE WAVE SENSITIVITY TO ANISOTROPIC STRUCTURE

Chevrot (2006) gives the following coordinate-free expression of the perturbed elastic stiffness tensor for an anisotropic perturbation with hexagonal symmetry:

$$ \begin{aligned} \delta c_{jmnk} = & \delta \lambda \delta_{jm} \delta_{nk} + \delta \mu (\delta_{jn} \delta_{mk} + \delta_{jk} \delta_{mn}) \\ & + 2 \rho \alpha^2 \epsilon (\hat{s}_j \hat{s}_m \hat{s}_n \hat{s}_k - \delta_{jm} \hat{s}_n \hat{s}_k - \delta_{nk} \hat{s}_j \hat{s}_m) \\ & + \rho \alpha^2 \delta (\delta_{jm} \hat{s}_n \hat{s}_k + \delta_{nk} \hat{s}_j \hat{s}_m - 2 \hat{s}_j \hat{s}_m \hat{s}_n \hat{s}_k) \\ & + 2 \rho \beta^2 \gamma (2 \delta_{jm} \hat{s}_n \hat{s}_k + 2 \delta_{nk} \hat{s}_j \hat{s}_m - \delta_{jk} \hat{s}_m \hat{s}_n - \delta_{mn} \hat{s}_j \hat{s}_k - \delta_{mk} \hat{s}_j \hat{s}_n - \delta_{jn} \hat{s}_m \hat{s}_k), \end{aligned} \quad (1) $$

where $\hat{s}$ is the unit vector describing the orientation of the symmetry axis, $\lambda$ and $\mu$ are the Lamé parameters of the reference isotropic medium, and $\epsilon$, $\delta$ and $\gamma$ are based on an extension of the Thomsen parameters widely in use in the seismic exploration literature (Thomsen 1986; Mensch & Rasolofosaon 1997). While the above expression is true for arbitrary orientation of the symmetry axis, it is convenient to define the parameters $\delta\lambda$, $\delta\mu$, $\epsilon$, $\delta$ and $\gamma$ with respect to a locally defined coordinate system where the $x_3$ axis is defined to be along the axis of symmetry. We can then define

$$ \epsilon = (C_{11} - C_{33})/2(\lambda + 2\mu) \quad (2) $$

$$ \delta = (C_{13} - C_{33} + 2C_{44})/(\lambda + 2\mu) \quad (3) $$

$$ \gamma = (C_{66} - C_{44})/2\mu \quad (4) $$

$$ \delta\lambda = C_{11} - 2C_{66} - \lambda \quad (5) $$
$$ \delta\mu = C_{66} - \mu, \qquad (6) $$

where the $C_{l,l}$ represent elements of the Voigt matrix (see Appendix A) in the locally defined coordinate system. We should also note here that the form of eq. (1) specifically results from the parametrization choices in eqs (2)–(6). If, for example, $\delta\mu$ were instead defined relative to $C_{44}$, a different form of eq. (1) would result.

Although these anisotropic and isotropic perturbation terms are defined with respect to a specific choice of coordinate system relative to the symmetry axis, it is important to note that eq. (1) is correct for arbitrary axis orientation. The coordinate system used to define $\epsilon$, $\delta$ and $\gamma$ as above need only be defined locally.

We shall wish to develop scattering coefficients for $\epsilon$, $\delta$ and $\gamma$. We adopt the normalization and Fourier sign convention of Zhou et al. (2004), who write the Green's function for surface waves as:

$$ G^{rs} = \mathbf{p}_s \cdot \mathbf{p}_r^* P(\theta), \qquad (7) $$

where $r$ denotes receiver, $s$ the source, $\mathbf{p}$ is the polarization vector:

$$ \mathbf{p} = U(r)\hat{\mathbf{r}} - iV(r)\hat{\theta} + iW(r)\hat{\phi}, \qquad (8) $$

where $\hat{\mathbf{r}}$ is the radial unit vector, $\hat{\theta}$ is the unit vector on the surface of the sphere along the great-circle path in the direction of propagation, and $\hat{\phi} = \hat{\mathbf{r}} \times \hat{\theta}$. $P(\theta)$ describes the propagation effect over a path of length $\theta$ radians,

$$ P(\theta) = \frac{e^{-i(\nu\theta + \pi/4 - n\pi/2)}}{|8\pi\nu \sin\theta|^{\frac{1}{2}}}, \qquad (9) $$

for $\nu = \ell + \frac{1}{2}$ and Maslov index $n$. If we choose our coordinate system such that the earthquake is at the pole $\theta = 0$, then $\theta$ is the spherical colatitude. We adopt this system, with $\phi$ as the longitudinal coordinate.

The first order (Born) approximation to the wave generated by a unit force in the $l$-direction at $s$ is:

$$ \delta G_{il}^{rs} = - \int \partial_m G_{ij}^{rx} \delta c_{jmnk} \partial_n G_{kl}^{xs} d\mathbf{x}. \qquad (10) $$

We need the gradients of the Green's tensor:

$$ \partial_m G_{ij}^{rx} = \partial_m [p_{xj}^{*} p_{ri} P(\theta'')] \qquad (11) $$

$$ \partial_n G_{kl}^{xs} = \partial_n [p_{sl}^* p_{sk} P(\theta')], \qquad (12) $$

where a single prime denotes the incoming wave at $x$, and a double prime the outgoing wave. We can evaluate the derivatives in (11) and (12) with the expressions

$$ \partial_m P(\theta'') = i\nu''r^{-1}P(\theta'')\hat{\theta}_m'' \qquad (13) $$

$$ \partial_n P(\theta') = -i\nu'r^{-1}P(\theta')\hat{\theta}_n' \qquad (14) $$

$$ \partial_m p_{xj}'' = \hat{r}_m(\hat{r}_j \dot{U}'' + i\hat{\theta}_j'' \dot{V}'' - i\hat{\phi}_j'' \dot{W}'') \qquad (15) $$

$$ \partial_n p_{xk}' = \hat{r}_n(\hat{r}_k \dot{U}' - i\hat{\theta}_k' \dot{V}' + i\hat{\phi}_k' \dot{W}') \qquad (16) $$

The unit vectors in expressions (15) and (16) are defined such that $\hat{\theta}'$ is the unit vector pointing along the azimuth of the incoming wave at a scattering point in the Born integral (eq. 10), and $\hat{\theta}''$ is a unit vector along the azimuth of the outgoing wave (Fig. 1). $\hat{\phi}'$ and $\hat{\phi}''$ are then defined as the cross-products $\hat{\mathbf{r}} \times \hat{\theta}'$ and $\hat{\mathbf{r}} \times \hat{\theta}''$. We note that since we are considering surface waves in expressions (11) and (12), we neglect terms like $\partial_m \hat{\theta}$ arising from the covariance of the unit vectors on the surface of the sphere which are proportional to $r^{-1}$, as they will be small in comparison to the terms which also include the wavenumber $\nu$, which is assumed large for the domain relevant for surface waves (Snieder & Nolet 1987). A full derivation including these covariant terms can be accommodated using a derivation based on surface wave strain tensors (e.g. Dahlen & Tromp 1998, expression 11.39). This can be seen in the supplementary online material for this paper, but the expressions are more complicated, and the simpler expressions derived here give excellent results when numerically validated (Section 3.2).

With that, the $i$th component of the scattered signal is, due to a unit force in the $L$-direction:

$$ \delta u_i = \delta G_{il}^{rs} = - \int p_{ri} P(\theta'') P(\theta') p_{sl}^* \Omega d\mathbf{x}, \qquad (17) $$

with the scattering coefficient

$$ \Omega = [\partial_m p_{xj}'' + i\nu''r^{-1}\hat{\theta}_m'' p_{xj}''][\partial_n p_{xk}' - i\nu'r^{-1}\hat{\theta}_n' p_{xk}']\delta c_{jmnk} \\ = O_{mj} \mathcal{I}_{nk} \delta c_{jmnk}, \qquad (18) $$

where we define the factors for outgoing and incoming waves:

$$ O_{mj} = \hat{r}_m(\hat{r}_j \dot{U}'' + i\hat{\theta}_j'' \dot{V}'' - i\hat{\phi}_j'' \dot{W}'') + i\nu''r^{-1}\hat{\theta}_m''(\hat{r}_j U'' + i\hat{\theta}_j'' V'' - i\hat{\phi}_j'' W'') \qquad (19) $$

$$ \mathcal{I}_{nk} = \hat{r}_n(\hat{r}_k \dot{U}' - i\hat{\theta}_k' \dot{V}' + i\hat{\phi}_k' \dot{W}') - i\nu'r^{-1}\hat{\theta}_n'(\hat{r}_k U' - i\hat{\theta}_k' V' + i\hat{\phi}_k' W') \qquad (20) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Diagram demonstrating the geometry of the source (asterisk), receiver (triangle) and scattering point (square) as points projected on the surface of a unit sphere. Solid lines indicate great-circle paths, and dashed lines are lines of longitude. The unit vectors $\hat{\theta}'$ and $\hat{\theta}''$ are shown. At the source, the mode take-off angle $\sigma'$ is shown, and at the receiver, the azimuths of the radial component, $\xi$, and of the mode arriving from the scattering point, $\xi''$, are shown. In the inset, the azimuths $\psi'$ and $\psi''$ as well as the scattering angle $\eta = \psi'' - \psi'$.</figcaption></figure>

Following the derivation in Zhou *et al.* (2004), we can use the Green's tensor expression to calculate the moment tensor response, which can be written as an implied mode sum,

$$u = S \times P(\theta) \times \mathcal{R}, \qquad (21)$$

where $S$ is defined in eq. (2.22) of Zhou *et al.* (2004), and is a function of the moment tensor elements, source take-off angle, and $\mathbf{p}_s$, and $\mathcal{R} = \mathbf{p}_r \cdot \hat{\mathbf{v}}$, where $\hat{\mathbf{v}}$ is the unit vector describing the polarization of the seismometer at the receiver.

Similarly, the expression in eq. (17) can be used to calculate the perturbed moment tensor response,

$$\delta u = - \int S' \mathcal{R}'' P(\theta'') P(\theta') \Omega \, d\mathbf{x}, \qquad (22)$$

where $S'$ is $S$ evaluated for the incoming mode at take-off angle appropriate for the scattering point, and $\mathcal{R}''$ is defined as in Zhou *et al.* (2004) eq. (2.31), and accounts for the polarization of the outgoing mode at the receiver location with a sign difference due to the sign convention used for $\Omega$. Correcting some typographical errors from Zhou *et al.* (2004) eq. (2.22), we write

$$\begin{aligned} S' = & -i\omega^{-1} \left[ M_{rr} \dot{U}_s' + (M_{\theta\theta} + M_{\phi\phi}) r_s^{-1} \left( U_s' - \frac{1}{2} \nu' V_s' \right) \right] \\ & + \omega^{-1} (-1)^n (\dot{V}_s' - r_s^{-1} V_s' + \nu' r_s^{-1} U_s') (M_{r\phi} \sin \sigma' + M_{r\theta} \cos \sigma') \\ & + i\omega^{-1} \nu' r_s^{-1} V_s' \left[ M_{\theta\phi} \sin 2\sigma' + \frac{1}{2} (M_{\theta\theta} - M_{\phi\phi}) \cos 2\sigma' \right] \\ & + \omega^{-1} (-1)^n (\dot{W}_s' - r_s^{-1} W_s') (M_{r\theta} \sin \sigma' - M_{r\phi} \cos \sigma') \\ & + i\omega^{-1} \nu' r_s^{-1} W_s' \left[ \frac{1}{2} (M_{\theta\theta} - M_{\phi\phi}) \sin 2\sigma' - M_{\theta\phi} \cos 2\sigma' \right], \end{aligned} \qquad (23)$$

where $M_{rr}$, $M_{\theta\theta}$, $M_{\phi\phi}$, $M_{r\theta}$, $M_{r\phi}$ and $M_{\theta\phi}$ are the six independent elements of the moment tensor, $r_s$ is the source radius, the s subscript on radial eigenfunctions means that they are evaluated at the source radius, $\sigma'$ is the azimuth from south of the outgoing mode (Fig. 1), and $\omega$ is the angular frequency. This is defined for the displacement response but can be converted by multiplying a factor of $i\omega$ for velocity or $-\omega^2$ for acceleration.

The receiver polarization term can be written as

$$\mathcal{R}'' = \begin{cases} U_r'' & \text{vertical component} \\ -i V_r'' \cos(\xi'' - \xi) - i W_r'' \sin(\xi'' - \xi) & \text{radial component} \\ i W_r'' \cos(\xi'' - \xi) - i V_r'' \sin(\xi'' - \xi) & \text{transverse component,} \end{cases} \qquad (24)$$

where $\xi''$ is the azimuth of the incoming scattered mode, and $\xi$ is the azimuth of the radial component (Fig. 1).
For a parameter $X$, which in this parametrization can be $\delta\lambda$, $\delta\mu$, $\epsilon$, $\delta$ or $\gamma$, we can define the kernel as

$$K_X(\bar{\mathbf{x}}) = -S'\mathcal{R}''P(\theta'')P(\theta')\Omega_X(\bar{\mathbf{x}}), \quad (25)$$

where $\bar{\mathbf{x}}$ is a position vector. In this expression, we have used the notation $\Omega = \Omega_X\delta X$, where $\delta X$ is the perturbation to the elastic tensor contributed by a perturbation of parameter $X$. Similarly, once a model parametrization basis is decided, we can define the partial derivatives for use in an inversion as

$$\frac{\partial u_i}{\partial m_h} = \int K_X(\bar{\mathbf{x}}) f_h(\bar{\mathbf{x}}) d\bar{\mathbf{x}}, \quad (26)$$

where $f_h$ is a model basis function.

Explicit expressions of the interaction coefficients $\Omega_X$ are given in Table 1. A detailed derivation of each coefficient is provided in the

<table><caption>**Table 1.** Non-zero interaction coefficients for surface waves for perturbations to $\delta\lambda$, $\delta\mu$, $\epsilon$, $\delta$ and $\gamma$.</caption><tbody><tr><td>Rayleigh → Rayleigh</td><td>Ω<sub>δλ</sub>:<br>U''U' - r<sup>-1</sup>(v'U''V' + v''V''U') + v''v'r<sup>-2</sup>V''V'</td></tr><tr><td>Rayleigh → Rayleigh</td><td>Ω<sub>δμ</sub>:<br>2U''U' + cos ηr<sup>-1</sup>(v'V''U' + v''U''V') + cos ηV''V'<br>+ v''v'r<sup>-2</sup>[cos ηU''U' + 2 cos<sup>2</sup> ηV''V']</td></tr><tr><td>Rayleigh → Love</td><td>sin η[v'r<sup>-1</sup>W''U' + W''V' + 2v''v'r<sup>-2</sup> cos ηW''V']</td></tr><tr><td>Love → Rayleigh</td><td>- sin η[v'r<sup>-1</sup>U''W' + V''W' + 2v''v'r<sup>-2</sup> cos ηV''W']</td></tr><tr><td>Love → Love</td><td>cos ηW''W' + v''v'r<sup>-2</sup> cos(2η)W''W'</td></tr><tr><td colspan="2">Ω<sub>ε</sub> = 2ρα<sup>2</sup> (Ω<sub>ε</sub><sup>(1)</sup> + Ω<sub>ε</sub><sup>(2)</sup> + Ω<sub>ε</sub><sup>(3)</sup>)</td></tr><tr><td colspan="2">For Ω<sub>ε</sub><sup>(1)</sup>:</td></tr><tr><td>Rayleigh → Rayleigh</td><td>s<sub>r</sub><sup>2</sup> A''A' + is<sub>r</sub>r<sup>-1</sup>[s<sub>θ</sub><sup>2</sup>v''A''A' - s'<sub>θ</sub>v'A''A']<br>+ v''v'r<sup>-2</sup>s'<sub>θ</sub>s''A''A'</td></tr><tr><td>Rayleigh → Love</td><td>-is<sub>φ</sub><sup>2</sup> s''A'W'' + s<sub>r</sub>r<sup>-1</sup>[s''<sub>θ</sub>s''v''A'W'' - s'<sub>θ</sub>s''v'A'W'']<br>- i v''v'r<sup>-2</sup>s'<sub>θ</sub>s''<sub>φ</sub>A'W'</td></tr><tr><td>Love → Rayleigh</td><td>is<sub>r</sub><sup>2</sup>s'<sub>φ</sub>'W''A'' + s<sub>r</sub>s<sub>φ</sub>'r<sup>-1</sup>[s'<sub>θ</sub>v'W'A'' - s''<sub>θ</sub>v''W'A'']<br>+ i v''v'r<sup>-2</sup>s'<sub>θ</sub>s<sub>φ</sub>'s''W'A''</td></tr><tr><td>Love → Love</td><td>s<sub>θ</sub><sup>2</sup>s''<sub>φ</sub>'W''W'' + is<sub>r</sub>s<sub>φ</sub>'s''<sub>φ</sub>'r<sup>-1</sup>[s''<sub>θ</sub>v''W'W'' - s'<sub>θ</sub>v'W''W'']<br>+ v''v'r<sup>-2</sup>s<sub>φ</sub>'s''<sub>φ</sub>'s''<sub>θ</sub>s''<sub>θ</sub>W''W''</td></tr><tr><td colspan="2">For Ω<sub>ε</sub><sup>(2)</sup>:</td></tr><tr><td>Rayleigh → Rayleigh</td><td>-(U'' - v''r<sup>-1</sup>V'')(s<sub>r</sub>A' - i v'r<sup>-1</sup>s'<sub>θ</sub>A')</td></tr><tr><td>Love → Rayleigh</td><td>-(U'' - v''r<sup>-1</sup>V'')(is<sub>φ</sub>'s<sub>r</sub>W' + v'r<sup>-1</sup>s'<sub>θ</sub>s''<sub>φ</sub>'W')</td></tr><tr><td colspan="2">For Ω<sub>ε</sub><sup>(3)</sup>:</td></tr><tr><td>Rayleigh → Rayleigh</td><td>-(U' - v'r<sup>-1</sup>V')(s<sub>r</sub>A'' + i v''r<sup>-1</sup>s''<sub>θ</sub>A'')</td></tr><tr><td>Rayleigh → Love</td><td>-(U' - v'r<sup>-1</sup>V')(-is<sub>φ</sub>''s<sub>r</sub>W'' + v''r<sup>-1</sup>s''<sub>φ</sub>s''<sub>θ</sub>W'')</td></tr><tr><td colspan="2">Ω<sub>δ</sub> = 2ρα<sup>2</sup> (Ω<sub>δ</sub><sup>(1)</sup> + Ω<sub>δ</sub><sup>(2)</sup> + Ω<sub>δ</sub><sup>(3)</sup>)</td></tr><tr><td colspan="2">Ω<sub>δ</sub><sup>(1)</sup> = -1/2Ω<sub>ε</sub><sup>(2)</sup></td></tr><tr><td colspan="2">Ω<sub>δ</sub><sup>(2)</sup> = -1/2Ω<sub>ε</sub><sup>(3)</sup></td></tr><tr><td colspan="2">Ω<sub>δ</sub><sup>(3)</sup> = -Ω<sub>ε</sub><sup>(1)</sup></td></tr><tr><td colspan="2">Ω<sub>γ</sub> = -2ρβ<sup>2</sup> Σ<sub>i=1</sub><sup>6</sup> Ω<sub>γ</sub><sup>(i)</sup></td></tr><tr><td colspan="2">Ω<sub>γ</sub><sup>(1)</sup> = 2Ω<sub>ε</sub><sup>(2)</sup></td></tr><tr><td colspan="2">Ω<sub>γ</sub><sup>(2)</sup> = 2Ω<sub>ε</sub><sup>(3)</sup></td></tr><tr><td colspan="2">For Ω<sub>γ</sub><sup>(3)</sup>:</td></tr><tr><td>Rayleigh → Rayleigh</td><td>s<sub>r</sub><sup>2</sup>[U''U' + cos ηV''V'] - i v'r<sup>-1</sup>s<sub>r</sub>s''<sub>θ</sub>[U''U' + V''V' cos η]<br>+i v''r<sup>-1</sup>s<sub>r</sub>s''<sub>θ</sub>[U''U' + V''V' cos η]<br>+v''v'r<sup>-2</sup>s''<sub>θ</sub>s''<sub>θ</sub>[U''U' + V''V' cos η]</td></tr><tr><td>Rayleigh → Love</td><td>sin η[s<sub>r</sub><sup>2</sup>W''V' - i v'r<sup>-1</sup>s<sub>r</sub>s''<sub>θ</sub>W''V' + i v''r<sup>-1</sup>s<sub>r</sub>s''<sub>θ</sub>W''V'<br>+ v''v'r<sup>-2</sup>s''<sub>θ</sub>s''<sub>θ</sub>W''V']</td></tr><tr><td>Love → Rayleigh</td><td>- sin η[s<sub>r</sub><sup>2</sup>W''W' - i v'r<sup>-1</sup>s<sub>r</sub>s''<sub>θ</sub>W''W' + i v''r<sup>-1</sup>s<sub>r</sub>s''<sub>θ</sub>V''W'<br>+v''v'r<sup>-2</sup>s''<sub>θ</sub>s''<sub>θ</sub>V''W']</td></tr><tr><td>Love → Love</td><td>cos η[s<sub>r</sub><sup>2</sup>W''W' - i v'r<sup>-1</sup>s<sub>r</sub>s''<sub>θ</sub>W''W' + i v''r<sup>-1</sup>s<sub>r</sub>s''<sub>θ</sub>W''W'<br>+ v''v'r<sup>-2</sup>s''<sub>θ</sub>s''<sub>θ</sub>W''W']</td></tr></tbody></table>
<table><caption>**Table 1.** (Continued.)</caption><thead><tr><th></th><th>For Ω<sub>γ</sub><sup>(4)</sup>:</th></tr></thead><tbody><tr><td>Rayleigh → Rayleigh</td><td>s<sub>r</sub><sup>2</sup>[U''U' + v''v'r<sup>-2</sup> cos ηU''U'] - is<sub>r</sub>s'<sub>θ</sub>[U''V' + v''v'r<sup>-2</sup> cos ηU''V']<br>+s''<sub>θ</sub>s''<sub>θ</sub>[V''V' + v''v'r<sup>-2</sup> cos ηV''V']<br>+is<sub>r</sub>s''<sub>θ</sub>[V''U' + v''v'r<sup>-2</sup> cos ηV''U']</td></tr><tr><td>Rayleigh → Love</td><td>-s''<sub>φ</sub>s''<sub>θ</sub>[W''V' + v''v'r<sup>-2</sup> cos ηW''V']<br>-is''<sub>φ</sub>s''<sub>r</sub>[W''U' + v''v'r<sup>-2</sup> cos ηW''U']</td></tr><tr><td>Love → Rayleigh</td><td>is<sub>r</sub>s''<sub>φ</sub>[U''W' + v''v'r<sup>-2</sup> cos ηU''W']<br>-s''<sub>θ</sub>s''<sub>φ</sub>[V''W' + v''v'r<sup>-2</sup> cos ηV''W']</td></tr><tr><td>Love → Love</td><td>s''<sub>φ</sub>s''<sub>φ</sub>[W''W' + v''v'r<sup>-2</sup> cos ηW''W']</td></tr></tbody><thead><tr><th></th><th>For Ω<sub>γ</sub><sup>(5)</sup>:</th></tr></thead><tbody><tr><td>Rayleigh → Rayleigh</td><td>s<sub>r</sub><sup>2</sup>[U''U' + v''r<sup>-1</sup> cos ηU''V']<br>-iv'r<sup>-1</sup>s<sub>r</sub>s''<sub>θ</sub>[U''U' + v''r<sup>-1</sup> cos ηU''V']<br>+is<sub>r</sub>s''<sub>θ</sub>[V''U' + v''r<sup>-1</sup> cos ηV''V']<br>+v'r<sup>-1</sup>s''<sub>θ</sub>s''<sub>θ</sub>[V''U' + v''r<sup>-1</sup> cos ηV''V']</td></tr><tr><td>Rayleigh → Love</td><td>-is''<sub>φ</sub>s''<sub>r</sub>W''U' - iv''r<sup>-1</sup>s''<sub>φ</sub>s''<sub>r</sub> cos ηW''V'<br>-v'r<sup>-1</sup>s''<sub>φ</sub>s''<sub>θ</sub>(W''U' + v''r<sup>-1</sup> cos ηW''V')</td></tr><tr><td>Love → Rayleigh</td><td>-iv''r<sup>-1</sup>s''<sub>r</sub>s''<sub>θ</sub> sin ηV''W' + iv''v'r<sup>-2</sup>s''<sub>r</sub>s''<sub>θ</sub> sin ηU''W'<br>-v''v'r<sup>-2</sup>s''<sub>θ</sub>s''<sub>θ</sub> sin ηV''W' - v''r<sup>-1</sup>s''<sub>r</sub> sin ηU''W'</td></tr><tr><td>Love → Love</td><td>iv''r<sup>-1</sup>s''<sub>φ</sub>s''<sub>r</sub> sin ηW''W' + v''v'r<sup>-2</sup>s''<sub>θ</sub>s''<sub>φ</sub> sin ηW''W'</td></tr></tbody><thead><tr><th></th><th>For Ω<sub>γ</sub><sup>(6)</sup>:</th></tr></thead><tbody><tr><td>Rayleigh → Rayleigh</td><td>s<sub>r</sub><sup>2</sup>[U''U' + v''r<sup>-1</sup> cos ηU''V']<br>-iv'r<sup>-1</sup>s<sub>r</sub>s''<sub>θ</sub>[U''U' + v''r<sup>-1</sup> cos ηU''V']<br>+is<sub>r</sub>s''<sub>θ</sub>[V''U' + v''r<sup>-1</sup> cos ηV''V']<br>+v'r<sup>-1</sup>s''<sub>θ</sub>s''<sub>θ</sub>[V''U' + v''r<sup>-1</sup> cos ηV''V']</td></tr><tr><td>Rayleigh → Love</td><td>-is''<sub>φ</sub>s''<sub>r</sub>W''U' - iv''r<sup>-1</sup>s''<sub>φ</sub>s''<sub>r</sub> cos ηW''V'<br>-v'r<sup>-1</sup>s''<sub>φ</sub>s''<sub>θ</sub>(W''U' + v''r<sup>-1</sup> cos ηW''V')</td></tr><tr><td>Love → Rayleigh</td><td>-iv''r<sup>-1</sup>s''<sub>r</sub>s''<sub>θ</sub> sin ηV''W' + iv''v'r<sup>-2</sup>s''<sub>r</sub>s''<sub>θ</sub> sin ηU''W'<br>-v''v'r<sup>-2</sup>s''<sub>θ</sub>s''<sub>θ</sub> sin ηV''W' - v''r<sup>-1</sup>s''<sub>r</sub> sin ηU''W'</td></tr><tr><td>Love → Love</td><td>iv''r<sup>-1</sup>s''<sub>φ</sub>s''<sub>r</sub> sin ηW''W' + v''v'r<sup>-2</sup>s''<sub>θ</sub>s''<sub>φ</sub> sin ηW''W'</td></tr></tbody></table>

supplementary online material for this paper. In addition to the quantities already defined, the terms in Table 1 also depend on the scattering angle, η, defined by the outgoing azimuth at the scattering point, ψ", minus the incoming azimuth, ψ' (see Fig. 1). These expressions then give us a means to define the kernels for use in an inversion. In the table, we also use the terms

$$s_r = (\hat{s} \cdot \hat{r})$$

$$s'_{\theta} = (\hat{s} \cdot \hat{\theta}')$$

$$s''_{\phi} = (\hat{s} \cdot \hat{\phi}')$$

$$s''_{\theta} = (\hat{s} \cdot \hat{\theta}'')$$

$$s''_{\phi} = (\hat{s} \cdot \hat{\phi}'')$$

$$A' = s_r U' - is'_{\theta} V'$$

$$A'' = s_r U'' + is''_{\theta} V''.$$

(27)

The expressions for the scattered waveform δu derived here can be used directly to compute kernels for waveform inversion studies, as follows in Section 3, or be applied to dispersion data using the relationship between the phase perturbation δφ and the perturbation (Zhou et al. 2004):

$$\delta\phi = \operatorname{Im}\left(\frac{\delta u}{u}\right),$$

and similar expressions for amplitude and arrival angle (see also Nolet 2008).

We note that the surface wave mode approach used in this study is comparable with a complete normal mode sum approach for periods shorter than approximately 250 s. For shorter periods, the surface wave mode approach has advantages over normal modes in the application of Born scattering, as coupling need only be considered within the relatively small number of surface wave modes at a specific frequency, rather than across the entire set of normal modes.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Subsets of scattering and total $\delta\mu/\mu$ (A) and $\gamma$ kernels (B) for the time point shown in the bottom seismogram (2800 s) during the fundamental mode Rayleigh wave train recorded on the *L* component for slices at a depth of 150 km. A vertical symmetry axis is assumed for this calculation. The source has the mechanism shown and is at a depth of 33 km, and the receiver is at a distance of 100°. For each model parameter, Rayleigh to Rayleigh (RR), Rayleigh to Love (RL), Love to Rayleigh (LR) and Love to Love (LL) scattering is shown, with the total summed contribution at the bottom. For each plot, the maximum value for the colour scale is shown in parentheses.</figcaption></figure>

### 3 KERNEL CALCULATION

The derived expressions in Table 1 and eq. (25) can be directly implemented and the kernels can then be examined in order to gain further insight into the finite-frequency sensitivity of surface waveforms. For the following figures, kernel calculations were made for seismograms bandpass filtered between 0.005 and 0.01 Hz (periods from 100 to 200 s), with cosine tapers with corners at 0.006 and 0.009 Hz. A complete set of Love and Rayleigh surface wave modes at 128 frequencies up to 0.01667 Hz (60 s period) with phase velocities between 2.5 and 250 km s⁻¹ were used for the calculations. These modes were calculated in a simplified version of PREM (Dziewonski & Anderson 1981), with the crust and 220 km discontinuity removed, as in the numerical tests in Section 3.2. The seismograms are rotated so that its three components are vertical (Z), and the horizontal components along (L) and perpendicular (T) to the great-circle path from source to receiver. Analysing a waveform kernel for a time point within the fundamental mode Rayleigh wave train recorded on the *L* component through a model with a vertical symmetry axis (Fig. 2) shows a clear connection to the relatively simple path-average approximation (Woodhouse & Dziewonski 1984), commonly used to model surface waves. This approximation assumes a surface wave can be modelled by the average structure along the great-circle path between source and receiver. While such an approach does not require the calculation of a volumetric kernel, it can be visualized as a constant valued averaging kernel along the infinitesimal path. The Rayleigh to Rayleigh scattering that dominates the sensitivity within the fundamental mode Rayleigh wave approaches such a kernel, although there is some along-path pulsing which is a function of the coupling with higher-mode Rayleigh energy (e.g. Zhou *et al.* 2004). The other scattering terms for this portion of the seismogram are much smaller, but still illustrate some interesting physics. Note the asymmetry of the kernels, which is due to the source mechanism, with stronger Rayleigh to Love scattering north of the source where the maximum of the Rayleigh source radiation occurs, while the Love to Rayleigh and Love to Love terms are larger to the east, where the Love radiation maximum occurs. For this case, where the symmetry axis is vertical, the $\gamma$ kernels are the opposite sign of the $\delta\mu/\mu$ kernels and are roughly twice the amplitude. We can gain some additional physical insight into the significance of these kernels by considering the Christoffel matrix ($\delta B_{jl}$ Dahlen & Tromp 1998, p. 82), which has the property of
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Same as Fig. 2, but on the transverse component at a time point during the fundamental Love wave train (2275 s). Slices in this figure are shown at a depth of 50 km.</figcaption></figure>

having eigenvectors that define the polarity of elastic waves with eigenvalues related to the squared velocity of the particular elastic wave. It is calculated by the expression

$$ \rho \delta B_{jl} = \delta c_{ijkl} \hat{k}_i \hat{k}_j, \qquad (28) $$

where $\hat{k}$ is the unit vector in the direction of propagation. We can substitute $\delta c$ from eq. (1) in the above expression for the case of a vertical symmetry axis. In this case, the eigenvalue for the principal axis perpendicular to $\hat{k}$ in the vertical plane is perturbed by a value of $(\delta\mu - 2\rho\beta^2\gamma)/\rho$. If we convert this to a perturbation $\delta V_{SV}/V_{SV}$, we obtain $\frac{1}{2}\delta\mu/\mu - \gamma$, consistent with Fig. 2.

The $\delta\mu/\mu$ kernel for a time point within the fundamental Love wave train (Fig. 3) is a little more complicated. It is dominated by the Love to Love coupling, as might be expected, but there is more along path variation due to coupling with higher modes than in the case of the fundamental Rayleigh wave. This is not surprising due to the smaller separation in terms of group velocity between the fundamental mode and overtones of Love waves compared to Rayleigh waves. The $\gamma$ kernels are much smaller in amplitude, and are actually dominated by off-path contributions due to Rayleigh to Rayleigh scattering. Once again, we can explain this by using a Christoffel matrix calculation and determining in this case the $\delta V_{SH}/V_{SH}$ from the perturbation to the eigenvalue related to the principal axis perpendicular to $\hat{k}$ in the horizontal plane. For the case of horizontal propagation with a vertical symmetry axis, this is simply equal to $\frac{1}{2}\delta\mu/\mu$, and has no dependence on $\gamma$. Finally for a time point during the Love overtones (Fig. 4), the $\delta\mu/\mu$ kernels have even more along-path variation, as might be expected for overtones, which can be considered as interference patterns of multiply free-surface reflected S phases, while the $\gamma$ kernels still demonstrate the importance of Rayleigh to Rayleigh and Rayleigh to Love scattering due to the stronger $\gamma$ influence on $V_{SV}$, and therefore, Rayleigh waves. Examples of kernels for all five elastic parameters for the time points in Figs 2–4 are available in the online supplementary material.

## 3.1 Non-linearity of kernels

As discussed further in Section 4, the kernels for this parametrization lead to a non-linear inversion since the values of $\Omega_X$ depend on the orientation of the symmetry axis $\hat{s}$ in the current model. We can examine how this non-linearity is expressed, by seeing how the kernels discussed above change as we vary the orientation angles of the symmetry axis: $\theta_s$, the angle from vertical and $\phi_s$, the azimuthal orientation.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Same as Fig. 3 at a time point during the Love overtones wave train (1500 s), also shown at a depth of 50 km.</figcaption></figure>

We consider values of $\theta_s$ of 0°, 45° and 90° (vertical, intermediate and horizontal), and values of $\phi_s$ of -45°, 0° and 45°, measured counter-clockwise from south. For the source-receiver configuration in Figs 2–6, the path is oriented roughly southwest-northeast, although it varies from more north-south near the source to more east-west near the receiver, which means the three $\phi_s$ values are approximately parallel, 45° oblique, and perpendicular to the path.

For a $\gamma$ kernel within the fundamental Rayleigh wave (Fig. 5), the kernel pattern is quite similar for all angles. If the axis of symmetry is oriented parallel with the path, the kernel is effectively unchanged for all values of $\theta_s$, while the amplitude of the kernel is strongly reduced for a horizontal axis with a southeast orientation, and the pattern becomes slightly more complicated as other forms of scattering besides Rayleigh to Rayleigh become more important. Once again, this behaviour is consistent with predictions from a Christoffel matrix calculation, which would predict no $\gamma$ contribution to $V_{SV}$ for a horizontal axis perpendicular to the path. This variation of sensitivity as a function of azimuth is also consistent with the theoretical expectation that Rayleigh wave sensitivity should be strongest for 2$\Psi$ terms (Smith & Dahlen 1973), where $\Psi$ is the path azimuth.

For Love waves, on the other hand, kernel amplitude depends strongly on the value of $\theta_s$ (Fig. 6). As discussed before, the kernel is near zero and quite complicated for a vertical axis, but the $\gamma$ kernel amplitude increases strongly as $\theta_s$ goes to horizontal. The pattern approaches that of the $\delta\mu/\mu$ kernel (Fig. 3), but at twice the amplitude and the opposite sign for all values of $\phi_s$. As for the Rayleigh waves, the azimuthal dependence is as expected from a Christoffel calculation, for which the $\gamma$ contribution to $V_{SH}$ is the same for paths parallel and perpendicular to the symmetry axis. This is also consistent with the theoretical prediction that 2$\Psi$ terms should be small for Love waves (Smith & Dahlen 1973), but does not clearly demonstrate the 4$\Psi$ sensitivity predicted in that study.

## 3.2 Numerical validation

In order to ensure the accuracy of the derivations in Section 2, we performed a numerical validation of the five elastic kernels for a simple test case (Fig. 7). For the test case we assumed a simple uniform perturbation in the whole mantle and evaluated the differential seismograms with respect to seismograms in the initial reference model using a spectral element method approach which also includes coupling to a normal mode solution for the core for computational efficiency (CSEM, Capdeville et al. 2002; Chaljub et al. 2003). The reference model for the
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Figure showing non-linear dependence of $\gamma$ kernel for *L* component during fundamental Rayleigh wave train. The top panel shows the 150 km slice for a vertical axis, as in Fig. 2. The middle panels show kernels for the same time point for the case of an axis dipping 45° from vertical ($\theta_s = \pi/4$) to the southwest (left-hand panel), south (middle panel) and southeast (right-hand panel), $\phi_s = -\pi/4, 0$ and $+\pi/4$, respectively. The bottom panels show the same azimuths of $\phi_s$ for a horizontal axis $\theta_s = \pi/2$.</figcaption></figure>

test case is the simplified version of PREM, as in the previous sections, which is used for numerical efficiency. As in the example kernel calculations, all seismograms shown are bandpassed between 0.005 and 0.01 Hz (periods from 100 to 200 s), with cosine tapers with corners at 0.006 and 0.009 Hz. Velocity seismograms shown in Fig. 7 are at a distance of 40° (4450 km), with both source and receivers located along the equator. A source mechanism with a vertical dip-slip plane oriented 45° from the path azimuth was chosen so as to be away from both maximums and nodes in the Love and Rayleigh wave radiation patterns.

When we compare the 1-D CSEM synthetics and the surface wave mode synthetics within the arrival window of the high amplitude surface wave energy, they are in quite good agreement. However, there is visible amplitude mismatch in the ScS2 phase arriving just at the end of the *T* component trace. This mismatch is progressively larger for higher multiple ScS phases (not shown). This is due to the limitation of the surface wave approach, primarily due to the breakdown of the far field approximation to $P(\theta)$ in eq. (9) for modes with small wavenumber $v$, as typical for near-vertical arriving body phases. This is an important limitation to be aware of when applying these kernels to actual data, although it is somewhat exaggerated by this test, as compared to real data, due to the neglect of attenuation in this test.

For the perturbed seismograms, we assumed a uniform symmetry axis orientation with $\theta_s = \pi/4$ and $\phi_s = \pi/3$. This orientation was chosen so as to be oblique to the propagation direction and to validate the kernels for non-zero values of $s_r, s_\theta$, and $s_\phi$. Uniform 0.5 per cent perturbations throughout the mantle of $\delta\lambda/\lambda, \delta\mu/\mu, \epsilon, \delta$ and $\gamma$ are considered. As shown in Fig. 7, there is quite good agreement between kernel predictions and numerical results for most cases. There is some mismatch noticeable for the $\epsilon$ and $\delta$ kernels on the *T* component, and to a lesser extent on the *L* component. This mismatch occurs during the arrival of ScS, which is during the fundamental mode for this distance in this reference model, and there is also a mismatch apparent in the differential traces at the very end of the displayed records where ScS2 energy is arriving. Once again, this is a function of the breakdown of the surface wave mode approach when dealing with near vertically arriving long period body wave phases. As this mismatch is greatest for the *T* component $\epsilon$ and $\delta$ kernels, which have a relatively small impact on Love waves (note the amplification factors in Fig. 7), such kernels should be sufficiently accurate for inversion purposes.
<figure><img src="image_7.png" /><figcaption>Figure 6. Same as Fig. 5 for 50 km depth slices and for the T component of data during the fundamental Love wave train.</figcaption></figure>

## 4 NON-LINEAR INVERSION STRATEGIES

In terms of inversion strategies, we may choose to either simultaneously invert for isotropic perturbations or to fix the isotropic structure and only invert for the three anisotropic parameters as well as the orientation of the symmetry axis. In the first case, we can simply also use the scattering coefficients due to perturbations in $\lambda$ and $\mu$, or alternatively $P$ and $S$ velocity ($\alpha$ and $\beta$) (e.g. Zhou et al. 2006). In the second case the $\delta\lambda$ and $\delta\mu$ terms are determined by the combination of $\epsilon$, $\delta$ and $\gamma$ and the definition of the isotropic average. If we choose the Voigt average, which is a logical choice for a parametrization based on an elastic stiffness tensor, we obtain additional scattering coefficients for the anisotropic parameters (Appendix B). Alternatively, we can choose to simply define the isotropic velocity by a fixed $\delta\lambda$ and $\delta\mu$ in eq. (1) (e.g. Calvet et al. 2006). This is advantageous in the sense that it does not introduce any additional terms for the anisotropic scattering coefficients that would act to correlate the anisotropic kernels with the isotropic kernels. However, it is equivalent to fixing the values of $C_{11}$ and $C_{66}$, which relate to elastic properties measured perpendicular to the symmetry axis, rather than fixing some average over all directions, such as the Voigt average, which may be physically more desirable.

It may not be intuitively obvious from examination of eq. (1) why the choice of the definition of the isotropic average matters or why perturbing only $\epsilon$, $\delta$ and $\gamma$ could affect the isotropic average velocity. It does, however, matter because the expression in eq. (1) defines $\delta\lambda$ and $\delta\mu$ relative to specific $C_{II}(C_{66}$ for $\delta\mu$ and and $C_{11} - 2C_{66}$ for $\delta\lambda)$. This means that perturbing, for example, $\epsilon$, while holding $\delta\lambda$ and $\delta\mu$ to 0, would fix the $P$ velocity along one axis perpendicular to the symmetry axis (related to $C_{11}$) to the original isotropic value, while perturbing it in all other directions. In general, most definitions of an isotropic average which measure elastic properties along other axes would then be changed. For mathematical completeness, we need to add in $\delta\lambda$ and $\delta\mu$ terms if we really want the Voigt average velocities fixed (Appendix B). For a different choice of isotropic average (e.g. Reuss or Voigt-Reuss-Hill), we would require different $\delta\lambda$ and $\delta\mu$ terms to hold the desired isotropic average fixed.

Due to the inclusion of the orientation of the symmetry axis in the expressions for $\Omega_X$, the kernels defined in this approach are non-linear in the sense that they depend on the model parameter $\hat{s}$. We can of course proceed to invert in an iterative fashion for the isotropic and anisotropic perturbations and then for the orientation angles of the symmetry axis, $\theta_s$ and $\phi_s$, using the expressions

$$ \frac{\partial u_i}{\partial \theta_s} = \epsilon \left( \frac{\partial K_{\epsilon}}{\partial s_j} \right) \left( \frac{\partial s_j}{\partial \theta_s} \right) + \delta \left( \frac{\partial K_{\delta}}{\partial s_j} \right) \left( \frac{\partial s_j}{\partial \theta_s} \right) + \gamma \left( \frac{\partial K_{\gamma}}{\partial s_j} \right) \left( \frac{\partial s_j}{\partial \theta_s} \right) \quad (29) $$
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Left-hand column shows comparison of SEM calculated velocity traces (solid black) and surface wave mode sum (dashed red) for vertical (Z, top), as well as horizontal component along (L, middle) and perpendicular to the great-circle path (T, bottom), scaled to the same peak to peak amplitude. The following columns show differential seismograms (perturbed minus reference) due to a uniform 0.5 per cent mantle perturbation of δλ/λ, δμ/μ, ε, δ, and γ from SEM calculation (solid black) and kernel prediction (dashed red). Each column is multiplied by the factor specified in parentheses relative to the 1-D traces in the left column. All traces shown are from 200 to 1900 s.</figcaption></figure>

$$ \frac{\partial u_i}{\partial \phi_s} = \epsilon \left( \frac{\partial K_{\epsilon}}{\partial s_j} \right) \left( \frac{\partial s_j}{\partial \phi_s} \right) + \delta \left( \frac{\partial K_{\delta}}{\partial s_j} \right) \left( \frac{\partial s_j}{\partial \theta_s} \right) + \gamma \left( \frac{\partial K_{\gamma}}{\partial s_j} \right) \left( \frac{\partial s_j}{\partial \phi_s} \right), \quad (30) $$

where $s_j$ is an implied summation of the components of $\hat{s}$ in the $\hat{r}$, $\hat{\theta}$ and $\hat{\phi}$ directions. Note that $\delta\lambda$ and $\delta\mu$ terms do not enter the expression, as these kernels do not depend on $\hat{s}$. The partial derivatives of the kernels, $K_X$, with respect to the components of $\hat{s}$ can be evaluated directly from the expressions in Table 1, using the relations

$$ s'_{\theta} = \cos \psi's_{\theta} + \sin \psi's_{\phi} \quad (31) $$

$$ s'_{\phi} = -\sin \psi's_{\theta} + \cos \psi's_{\phi}, \quad (32) $$

where $\psi'$ is the azimuth of the incoming wave measured counter-clockwise from south (Fig. 1). The same expressions hold for $s''_{\theta}$ and $s''_{\phi}$ by replacing all the accents. The partial derivatives of the components of $\hat{s}$ with respect to the orientation angles are

$$ \frac{\partial s_r}{\partial \theta_s} = -\sin \theta_s $$

$$ \frac{\partial s_r}{\partial \phi_s} = 0 $$

$$ \frac{\partial s_{\theta}}{\partial \theta_s} = \cos \theta_s \cos \phi_s $$

$$ \frac{\partial s_{\theta}}{\partial \phi_s} = -\sin \theta_s \sin \phi_s $$

$$ \frac{\partial s_{\phi}}{\partial \theta_s} = \cos \theta_s \sin \phi_s $$

$$ \frac{\partial s_{\phi}}{\partial \phi_s} = \sin \theta_s \cos \phi_s. \quad (33) $$

However, like all non-linear inversions, this strategy may fail if the starting model is not sufficiently close to the true model. Its also important to repeat that the non-linearity of the problem, and the iterative approach, proposed, only accounts for the dependence of the kernels upon the orientation of the axis in the previous model iteration. The surface wave mode calculation, however, is only performed in the 1-D reference
model. This is not a fully non-linear inversion, in which all kernels are iteratively updated to reflect the 3-D model, which would require a computationally expensive numerical calculation.

We may also wish to reduce the number of parameters using scaling relationships. Most surface wave-based shear velocity models choose to scale perturbations in *P* velocity (to which the surface waves have less sensitivity) to those in *S* velocity, and there is ample empirical evidence for the choice for such scaling. It is less obvious how to go about scaling anisotropic parameters. There is evidence for scaling of anisotropic ratios in samples taken of natural peridotites (Montagner & Anderson 1989), and these ratios have been used in previous anisotropic global inversions (Panning & Romanowicz 2006). The problem with this approach is that it is difficult to know if the samples are representative of the mantle in general. Another option is to predict scaling based on geodynamic flow calculations and mineral physics based deformation mechanisms (Becker *et al.* 2006), which show good linear correlation between ε, δ and γ, as well as supporting the choice of using hexagonal symmetry to represent most of the mantle anisotropy. This suggests a parametrization for modelling surface wave and splitting data using just μ, or *V*<sub>S</sub> to more closely match the actual sensitivity of the data, as well as γ and the orientation angles of the symmetry axis, θ<sub>s</sub> and φ<sub>s</sub>.

Another important consideration for such a non-linear inversion is computational efficiency. As derived here, an iterative inversion would require the relatively time-consuming evaluation of the volume integral in eq. (26) after each update of the orientation of the symmetry axis. However, as discussed in Appendix A, we can limit the evaluation of this integral to one time in the reference model for the 21 independent elastic coefficients, *C*<sub>II</sub>, and then we can assemble the desired ε, δ and γ kernels for each updated symmetry axis orientation as linear combinations of the *C*<sub>II</sub> kernels.

# 5 LINEAR INVERSION STRATEGIES

With this in mind, a more desirable strategy may be to define an approximate linear inversion to obtain an adequate starting point for the non-linear inversion. There are many previous studies that define a series of linear kernels that represent the sensitivity of surface waves to a general anisotropic model as the sum of terms related to a radially anisotropic model (hexagonal symmetry with a vertical axis of symmetry) as well as other terms related to sinusoidal 2Ψ and 4Ψ dependence, where Ψ is an angle in the $\hat{\theta} - \hat{\phi}$ plane between the path azimuth and an axis of the coordinate system (e.g. Montagner & Nataf 1986; Sieminski *et al.* 2007a). In the full expressions, this in general leads to a large number of linear parameters [13 in the case of Montagner & Nataf (1986)], although in practice most studies have neglected the 4Ψ terms, and only considered the radial anisotropy terms and two terms ($G_c$ and $G_s$) related to the cos 2Ψ and sin 2Ψ terms, respectively (e.g. Simons *et al.* 2002). One potential approach then would be to use these terms to invert for the best-fitting combination of the strength of radial anisotropy and the strength and orientation of the 2Ψ anisotropy. This then could be used to define a best-fitting symmetry axis by considering the linear radial anisotropy and 2Ψ terms as the projections of the symmetry axis on the vertical and horizontal axes, respectively. This approximate symmetry could then be used as a starting model for an iterative non-linear inversion.

For convenience, we also derive expressions for the surface wave sensitivity kernels for the cases of fixed vertical or horizontal symmetry axis (Appendices C and D) using the parametrization defined by eq. (1). As expected, the case of a horizontal axis of symmetry can be defined by a linear dependence on azimuthally independent terms and the sinusoids of 2Ψ and 4Ψ. However, there are complications inherent in using the path azimuth, Ψ, in the definition of the kernels in a finite-frequency Born scattering formalism, where structure off the geometric path between source and receiver is taken into account assuming a scattering which changes the azimuth of energy propagation. For this reason, we instead choose to use the angle φs in our derivations in Appendix D, which defines the azimuthal orientation of the anisotropic symmetry axis, rather than the path azimuth. Of course an azimuthal dependence that depends on a linear combination of sin 2Ψ and cos 2Ψ can be equivalently represented by a different linear combination of sin 2φs and cos 2φs, and likewise for 4Ψ terms. Using these expressions, we could choose to linearly invert for a model that is simply the sum of sensitivity to a vertical axis of symmetry and a horizontal axis of symmetry. Using the same scaling relationships as above, this would reduce us to one isotropic parameter (e.g. *V*<sub>S</sub>), plus 6 anisotropic parameters (γ<sub>v</sub>, from the vertical symmetry axis, plus five terms from the horizontal symmetry axis, relating to azimuthally independent terms plus the sines and cosines of 2φs and 4φs). In comparison with the non-linear approach, this neglects all terms which contain both *s*<sub>*p*</sub> as well as *s*<sub>*θ*</sub> or *s*<sub>*φ*</sub>. Also the number of parameters is likely still too large to reasonably resolve. However, as in previous studies, we can choose to only model the 2φs sensitivity of the horizontal medium, leaving us with a total of three anisotropic parameters (e.g. γ<sub>v</sub>, γ<sub>c</sub> and γ<sub>s</sub>, which relate to the vertical axis, cos 2φs and sin 2φs terms, respectively). Future experiments with synthetic as well as real data will be necessary to resolve whether such an approach which neglects many terms is sufficient to develop a close enough starting model for the non-linear inversion for strength and orientation of anisotropy, which would include all those terms.

# 6 DISCUSSION AND CONCLUSIONS

The compact sensitivity kernels derived in this study offer a number of distinct advantages for the inversion of surface wave data in terms of the most important—and resolvable—anisotropic characteristics of the Earth. The model parametrization is reduced to five elastic parameters and the direction of the symmetry axis. If one seeks only an anisotropic perturbation on top of a first inversion attempt for an isotropic model, the number of elastic parameters is reduced to the remaining three anisotropic parameters ε, γ and δ.

The number of anisotropic parameters can be further reduced (to one) through the predicted and observed correlations between ε, γ and δ (Montagner & Anderson 1989; Becker *et al.* 2006). Parameter reduction does not need to be obtained by throwing out large numbers of parameters assumed to generally have small amplitudes in the real Earth, or little influence on most data, which may not always be justified
when considering the finite-frequency sensitivity of the data (e.g. Sieminski *et al.* 2007a,b). Also, this parametrization gives great latitude for incorporating other *a priori* information about the expected anisotropy, such as constraining $\gamma$ to be positive or negative, corresponding to the two families of fabrics identified in Becker *et al.* (2006), or constraining the axis to be horizontal or vertical.

Finally, this approach gives us a consistent theoretical framework for combining surface wave and SKS splitting data. These kernels and two proposed approaches for modelling shear splitting intensity (Chevrot 2006; Long *et al.* 2008) are based on equivalent formulations of the elastic perturbation and allow for the calculation of 3-D finite-frequency kernels. Such data sets have already shown great promise when combined in modelling efforts (Marone & Romanowicz 2007). Using finite-frequency kernels for both data sets in a consistent theoretical framework should improve the modelling, and when combined with appropriate error estimates should greatly simplify the process of relative weighting of the different data types.

The approach, however, comes with a possible disadvantage, due to the non-linearity of the kernels with respect to orientation of the symmetry axis. While the examples shown in Figs 5 and 6 indicate some stability of kernel patterns, the amplitudes can vary greatly depending on orientation. If there is not adequate data coverage, there could certainly be significant tradeoffs for example between strength of anisotropy and axis orientation. However, the differences in these tradeoffs between Love and Rayleigh waves, and for different azimuths of Rayleigh wave propagation suggest that such difficulties can be overcome with adequate data. Shear splitting measurements provide yet another complementary constraint. However, any non-linear inversion has the potential to be strongly dependent on starting model. This suggests that any approach can benefit from considering several starting models. Deriving these starting models from both linearized seismic inversion as well as families of geodynamic flow models would seem to be a hopeful avenue for constraining model stability and confidence. Further tests with both synthetic and real data will be necessary to determine what can be gained by this approach.

The models of anisotropic seismic velocity structure, in particular the magnitude of anisotropy inferred and the direction of the symmetry axis, are important for the understanding of dynamic processes in the mantle. By reducing the parameter space to a small number of parameters using empirical correlations, we bring the woefully underdetermined problem for full anisotropy back to a more feasible inverse problem while incorporating prior knowledge about the likely nature of anisotropy. By formulating the problem using a finite-frequency approach based on the Born approximation, we enhance resolvability by making use of the frequency-dependent width of the sensitivity. Chevrot & Zhao (2007) point out that such enhanced resolution is only obtained if the kernels are sampled sufficiently densely, which makes the computation more expensive, but not out of reach of a moderate Beowulf cluster. We provide source code for the computation of the kernels in the electronic Appendix.

Application of the theoretical work presented in this paper should show how difficult it is to overcome the non-linearity of the kernels, but this is reserved for future work. In this paper we proposed a number of inversion strategies that may allow for stable results, as well as practical approaches to reduce the number of parameters.

## ACKNOWLEDGMENTS

All figures except Fig. 1 were made using GMT (Wessel & Smith 1998). This work was funded through NSF grants EAR 0309298 and EAR 0105387. M.P. was also funded through the Princeton University Council on Science and Technology.

## REFERENCES

* Abt, D. & Fischer, K., 2008. Resolving three-dimensional anisotropic structure with shear-wave splitting tomography, *Geophys. J. Int.*, **173**, 859–886.

* Babuska, V. & Cara, M., 1991. *Seismic Anisotropy in the Earth*, Kluwer Academic Press, Boston, Massachusetts.

* Becker, T., Kellogg, J., Ekström, G, & O'Connell, R., 2003. Comparison of azimuthal seismic anisotropy from surface waves and finite strain from global mantle-circulation models, *Geophys. J. Int.*, **155**, 696–714.

* Becker, T., Chevrot, S., Schulte-Pelkum, V & Blackman, D., 2006. Statistical properties of seismic anisotropy predicted by upper mantle geodynamic models, *J. geophys. Res.*, **111**, doi:10.1029/2005JB004095.

* Becker, T., Kustowski, B., & Ekström, G., 2007. Radial seismic anisotropy as a constraint for upper mantle rheology, *Earth planet. Sci. Lett.*, **267**, 213–227.

* Beghein, C. & Trampert, J., 2004. Probability density functions for radial anisotropy for fundamental mode surface wave data and the Neighbourhood Algorithm, *Geophys. J. Int.*, **157**, 1163–1174.

* Calvet, M., Chevrot, S., & Souriau, A., 2006. P-wave propagation in transversely isotropic media I. Finite-frequency theory, *Phys. Earth planet. Inter.*, **156**, 12–20.

* Capdeville, Y., Chaljub, E., Vilotte, J.-P., & Montagner, J.-P., 2002. Coupling the spectral element method with a modal solution for elastic wave propagation in global Earth models, *Geophys. J. Int.*, **152**, 34–66.

* Chaljub, E., Capdeville, Y., & Vilotte, J.-P., 2003. Solving elastodynamics

in a fluid-solid heterogeneous sphere: a parallel spectral element approximation on non-conforming grids, *J. Comp. Phys.*, **187**, 457–491.

* Chevrot, S., 2006. Finite-frequency vectorial tomography: a new method for high-resolution imaging of upper mantle anisotropy, *Geophys. J. Int.*, **165**, 641–657.

* Chevrot, S. & Zhao, L., 2007. Multiscale finite frequency Rayleigh wave tomography of the Kaapvaal craton, *Geophys. J. Int.*, **169**, 201–215.

* Dahlen, F. & Tromp, J., 1998. *Theoretical Global Seismology*, Princeton University Press, Princeton, New Jersey.

* Dahlen, F., Hung, S.-H., & Nolet, G., 2000. Frechet kernels for finite-frequency traveltimes—I. Theory, *Geophys. J. Int.*, **141**, 157–174.

* Debayle, E., Kennett, B., & Priestley, K., 2005. Global azimuthal seismic anisotropy and the unique plate-motion deformation of Australia, *Nature*, **433**, 509–512.

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

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

* Fouch, M. & Fischer, K., 1996. Mantle anisotropy beneath northwest Pacific subduction zones, *J. geophys. Res.*, **101**, 15 987–16 002.

* Grand, S., 1997. Global seismic tomography: a snapshot of convection in the Earth, *GSA Today*, **7**, 1–7.

* Gu, Y., Dziewonski, A., & Ekström, G., 2003. Simultaneous inversion for mantle shear velocity and topography of transition zone discontinuities, *Geophys. J. Int.*, **154**, 559–583.
Karato, S.-I., 1998. Some remarks on the origin of seismic anisotropy in the D'layer, *Earth Planets Space*, **50**, 1019–1028.

Kendall, J.-M. & Silver, P., 1996. Constraints from seismic anisotropy on the nature of the lowermost mantle, *Nature*, **381**, 409–412.

Kustowski, B., Ekström, G., & Dziewonski, A., 2008. The anisotropic shear-wave velocity structure of the Earth's mantle: a global model, *J. geophys. Res.*, doi:10.1029/2007JB004937, in press.

Li, C., van der Hilst, R, & Toksöz, M., 2006. Constraining P-wave velocity variation in the upper mantle beneath Southeast Asia, *Phys. Earth planet. Inter.*, **154**, 180–195.

Long, M., de Hoop, M., & van der Hilst, R., 2008. Wave-equation shear wave splitting tomography, *Geophys. J. Int.*, **172**, 311–330.

Marone, F. & Romanowicz, B., 2007. The depth distribution of azimuthal anisotropy in the continental upper mantle, *Nature*, **447**, 198–201.

Masters, G., Laske, G., Bolton, H, & Dziewonski, A., 2000. The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle: implications for chemical and thermal structure, in *Earth's Deep Interior*, Vol. AGU Monograph 117, eds Karato, S.-I., Forte, A., Liebermann, R., Masters, G. & Stixrude, L., Amer. Geophys. Union, Washington, DC.

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

Mensch, T. & Rasolofosaon, P., 1997. Elastic-wave velocities in anisotropic media of arbitrary symmetry—generalization of Thomsen's parameters $\epsilon, \delta$, and $\gamma$, *Geophys. J. Int.*, **128**, 43–64.

Montagner, J.-P. & Anderson, D., 1989. Petrological constraints on seismic anisotropy, *Phys. Earth planet. Inter.*, **54**, 82–105.

Montagner, J.-P. & Jobert, N., 1988. Vectorial tomography: II. application to the Indian Ocean, *Geophys. J. R. astr. Soc.*, **94**, 295–307.

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., 1991. Global upper mantle tomography of seismic velocities and anisotropies, *J. geophys. Res.*, **96**, 20337–20351.

Montelli, R., Nolet, G., Dahlen, F., Masters, G., Engdahl, E., & Hung, S.-H., 2004. Finite-frequency tomography reveals a variety of plumes in the mantle, *Science*, **303**(5656), 338–343.

Montelli, R., Nolet, G., Dahlen, F., & Masters, G., 2006. A catalogue of deep mantle plumes: new results from finite-frequency tomography, *Geochem. Geophys. Geosys.* (G3), **7**, Q11007.

Nolet, G., 2008. *A Breviary for Seismic Tomography*, Cambridge Univ. Press, Cambridge, UK.

Panning, M. & Romanowicz, B., 2006. A three-dimensional radially anisotropic model of shear velocity in the whole mantle, *Geophys. J. Int.*, **167**, 361–379.

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

Sieminski, A., Liu, Q., Trampert, J., & Tromp, J., 2007a. Finite-frequency sensitivity of surface waves to anisotropy based upon adjoint methods, *Geophys. J. Int.*, **168**, 1153–1174.

Sieminski, A., Liu, Q., Trampert, J., & Tromp, J., 2007b. Finite-frequency sensitivity of body waves to anisotropy based upon adjoint methods, *Geo-phys. J. Int.*, **171**, 368–389.

Silver, P., 1996. Seismic anisotropy beneath the continents: probing the depths of geology, *Annu. Rev. Earth planet. Sci.*, **24**, 385–432.

Simons, F., van der Hilst, R., Montagner, J.-P. & Zielhuis, A., 2002. Multi-mode Rayleigh wave inversion for heterogeneity and azimuthal anisotropy of the australian upper mantle, *J. geophys. Res.*, **151**, 738–754.

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

Snieder, R. & Nolet, G., 1987. Linearized scattering of surface waves on a spherical Earth, *J. Geophys.*, **61**, 55–63.

Thomsen, L., 1986. Weak anisotropy, *Geophysics*, **51**, 1954–1966.

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

Vinnik, L., Makeyeva, L., Milev, A, & Usenko, Y., 1992. Global patterns of azimuthal anisotropy and deformations in the continental mantle, *Geo-phys. J. Int.*, **111**, 433–447.

Wessel, P. & Smith, W., 1998. New, improved version of Generic Mapping Tools released, *EOS Trans. Am. geophys. Un.*, **79**, 579.

Woodhouse, J. & Dziewonski, A., 1984. Mapping the upper mantle: Three dimensional modeling of Earth structure by inversion of seismic wave-forms, *J. geophys. Res.*, **89**, 5953–5986.

Zhao, L., Jordan, T, & Chapman, C., 2000. Three-dimensional Frèchet differential kernels for seismic delay times, *Geophys. J. Int.*, **141**, 558–576.

Zhou, Y., Dahlen, F., & Nolet, G., 2004. Three-dimensional sensitivity kernels for surface wave observables, *Geophys. J. Int.*, **158**, 142–168.

Zhou, Y., Nolet, G., Dahlen, F., & Laske, G., 2006. Global upper-mantle structure from finite-frequency surface-wave tomography, *J. geophys. Res.*, **111**, B04304, doi:10.1029/2005JB003677.

## APPENDIX A: PRACTICAL CONSIDERATIONS FOR NON-LINEAR INVERSIONS

Aside from the cases of fixed symmetry axis discussed in Appendices C and D, all of the above expressions are non-linear because the sensitivity to perturbations in $\epsilon, \delta$ and $\gamma$ depend upon the orientation of the symmetry axis $\hat{s}$. However, the calculation of the scattering terms over the whole scattering volume can be a computationally expensive step that would have to be repeated after each perturbation of the symmetry axis in the 3-D model, even though the radial eigenfunctions $U, V$ and $W$ would still be defined in the spherically symmetric reference model. This is a computationally undesirable route for the non-linear inversion approaches, as a non-linear inversion will almost certainly take multiple iterations to converge.

We can, however, make a non-linear inversion for the anisotropic parameters from eq. (1) and symmetry axis orientation more feasible, by defining the kernels as linear combinations of kernels which only need to be calculated once. A reasonable choice is to initially calculate the sensitivity to perturbations to the 21 independent coefficients of the symmetric Voigt matrix $C_{IJ}$, where $I$ and $J$ are indices that vary from 1 to 6 which represent two indices in the fourth order elastic tensor, and can be defined in spherical geometry as

$$1 = \theta\theta$$

$$2 = \phi\phi$$

$$3 = rr$$

$$4 = \phi r$$

$$5 = \theta r$$

$$6 = \theta\phi.$$

(A1)
Once we define the desired kernels (e.g. $K_\epsilon$, $K_\delta$ and $K_\gamma$) as a linear combination of the kernels with respect to $C_{IJ}$, $K_{IJ}$, we can speed up the inversion process by only updating the coefficients of that linear combination after each iteration.

Using eq. (1), we can define relationships between the perturbed $\delta C_{IJ}$ and the parameters $\delta\lambda$, $\delta\mu$, $\epsilon$, $\delta$ and $\gamma$ for a given symmetry axis $\hat{s}$,

$$ \delta C_{11} = \delta\lambda + 2\delta\mu + 2\rho\alpha^2\epsilon(s_{\theta}^4 - 2s_{\theta}^2) + \rho\alpha^2\delta(s_{\theta}^2 - 2s_{\theta}^4) \quad (A2) $$

$$ \delta C_{22} = \delta\lambda + 2\delta\mu + 2\rho\alpha^2\epsilon(s_{\phi}^4 - 2s_{\phi}^2) + \rho\alpha^2\delta(s_{\phi}^2 - 2s_{\phi}^4) \quad (A3) $$

$$ \delta C_{33} = \delta\lambda + 2\delta\mu + 2\rho\alpha^2\epsilon(s_r^4 - 2s_r^2) + \rho\alpha^2\delta(s_r^2 - 2s_r^4) \quad (A4) $$

$$ \delta C_{44} = \delta\mu + 2\rho\alpha^2\epsilon(s_{\phi}^2s_r^2 - 2\rho\alpha^2\delta(s_{\phi}^2s_r^2) - 2\rho\beta^2\gamma(s_{\phi}^2 + s_r^2)) \quad (A5) $$

$$ \delta C_{55} = \delta\mu + 2\rho\alpha^2\epsilon(s_{\theta}^2s_r^2 - 2\rho\alpha^2\delta(s_{\theta}^2s_r^2) - 2\rho\beta^2\gamma(s_{\theta}^2 + s_r^2)) \quad (A6) $$

$$ \delta C_{66} = \delta\mu + 2\rho\alpha^2\epsilon(s_{\theta}^2s_{\phi}^2 - 2\rho\alpha^2\delta(s_{\theta}^2s_{\phi}^2) - 2\rho\beta^2\gamma(s_{\theta}^2 + s_{\phi}^2)) \quad (A7) $$

$$ \delta C_{12} = \delta\lambda + 2\rho\alpha^2\epsilon(s_{\theta}^2s_{\phi}^2 - s_{\theta}^2 - s_{\phi}^2) + \rho\alpha^2\delta(s_{\theta}^2 + s_{\phi}^2 - 2s_{\theta}^2s_{\phi}^2) + 4\rho\beta^2\gamma(s_{\theta}^2 + s_{\phi}^2) \quad (A8) $$

$$ \delta C_{13} = \delta\lambda + 2\rho\alpha^2\epsilon(s_{\theta}^2s_r^2 - s_{\theta}^2 - s_r^2) + \rho\alpha^2\delta(s_{\theta}^2 + s_r^2 - 2s_{\theta}^2s_r^2) + 4\rho\beta^2\gamma(s_{\theta}^2 + s_r^2) \quad (A9) $$

$$ \delta C_{14} = 2\rho\alpha^2\epsilon(s_{\theta}^2s_{\phi}s_r - s_{\phi}s_r) + \rho\alpha^2\delta(s_{\phi}s_r - 2s_{\theta}^2s_{\phi}s_r) + 4\rho\beta^2\gamma s_{\phi}s_r \quad (A10) $$

$$ \delta C_{15} = 2\rho\alpha^2\epsilon(s_{\theta}^3s_r - s_{\theta}s_r) + \rho\alpha^2\delta(s_{\theta}s_r - 2s_{\theta}^3s_r) \quad (A11) $$

$$ \delta C_{16} = 2\rho\alpha^2\epsilon(s_{\theta}^3s_{\phi} - s_{\theta}s_{\phi}) + \rho\alpha^2\delta(s_{\theta}s_{\phi} - 2s_{\theta}^3s_{\phi}) \quad (A12) $$

$$ \delta C_{23} = \delta\lambda + 2\rho\alpha^2\epsilon(s_{\phi}^2s_r^2 - s_{\phi}^2 - s_r^2) + \rho\alpha^2\delta(s_{\phi}^2 + s_r^2 - 2s_{\phi}^2s_r^2) + 4\rho\beta^2\gamma(s_{\phi}^2 + s_r^2) \quad (A13) $$

$$ \delta C_{24} = 2\rho\alpha^2\epsilon(s_{\phi}^3s_r - s_{\phi}s_r) + \rho\alpha^2\delta(s_{\phi}s_r - 2s_{\phi}^3s_r) \quad (A14) $$

$$ \delta C_{25} = 2\rho\alpha^2\epsilon(s_{\phi}^2s_{\theta}s_r - s_{\theta}s_r) + \rho\alpha^2\delta(s_{\theta}s_r - 2s_{\phi}^2s_{\theta}s_r) + 4\rho\beta^2\gamma s_{\theta}s_r \quad (A15) $$

$$ \delta C_{26} = 2\rho\alpha^2\epsilon(s_{\phi}^3s_{\theta} - s_{\phi}s_{\theta}) + \rho\alpha^2\delta(s_{\phi}s_{\theta} - 2s_{\phi}^3s_{\theta}) \quad (A16) $$

$$ \delta C_{34} = 2\rho\alpha^2\epsilon(s_r^3s_{\phi} - s_r s_{\phi}) + \rho\alpha^2\delta(s_r s_{\phi} - 2s_r^3s_{\phi}) \quad (A17) $$

$$ \delta C_{35} = 2\rho\alpha^2\epsilon(s_r^3s_{\theta} - s_r s_{\theta}) + \rho\alpha^2\delta(s_r s_{\theta} - 2s_r^3s_{\theta}) \quad (A18) $$

$$ \delta C_{36} = 2\rho\alpha^2\epsilon(s_r^2s_{\theta}s_{\phi} - s_{\theta}s_{\phi}) + \rho\alpha^2\delta(s_{\theta}s_{\phi} - 2s_r^2s_{\theta}s_{\phi}) + 4\rho\beta^2\gamma s_{\theta}s_{\phi} \quad (A19) $$

$$ \delta C_{45} = 2\rho\alpha^2\epsilon s_r^2s_{\theta}s_{\phi} - 2\rho\alpha^2\delta s_r^2s_{\theta}s_{\phi} - 4\rho\beta^2\gamma s_{\theta}s_{\phi} \quad (A20) $$

$$ \delta C_{46} = 2\rho\alpha^2\epsilon s_{\theta}^2s_{\phi}s_r - 2\rho\alpha^2\delta s_{\theta}^2s_{\phi}s_r - 4\rho\beta^2\gamma s_{\theta}s_r \quad (A21) $$

$$ \delta C_{56} = 2\rho\alpha^2\epsilon s_{\theta}^2s_{\phi}s_r - 2\rho\alpha^2\delta s_{\theta}^2s_{\phi}s_r - 4\rho\beta^2\gamma s_{\phi}s_r \quad (A22) $$

Eqs (A2)–(A22) are then all of the form

$$ \delta C_{IJ} = c_{IJ}^{(\delta\lambda)}\delta\lambda + c_{IJ}^{(\delta\mu)}\delta\mu + c_{IJ}^{(\epsilon)}\epsilon + c_{IJ}^{(\delta)}\delta + c_{IJ}^{(\gamma)}\gamma \quad (A23) $$

For example, $c_{11}^{(\epsilon)}$ is defined as $2 \rho\alpha^2 (s_{\theta}^4 - 2s_{\theta}^2)$ from eq. (A2).

A local structural perturbation affects the first order scattered wavefield through the integrand $\sum_{IJ} K_{IJ} \delta C_{IJ}$, where $K_{IJ}$ is the sensitivity kernel for perturbations to $C_{IJ}$. We can then equivalently write

$$ \sum_{IJ} K_{IJ} \delta C_{IJ} = \sum_{IJ} K_{IJ} \left[ c_{IJ}^{(\delta\lambda)} \delta\lambda + c_{IJ}^{(\delta\mu)} \delta\mu + c_{IJ}^{(\epsilon)} \epsilon + c_{IJ}^{(\delta)} \delta + c_{IJ}^{(\gamma)} \gamma \right] \quad (A24) $$

$$ = K_{\delta\lambda}\delta\lambda + K_{\delta\mu}\delta\mu + K_{\epsilon}\epsilon + K_{\delta}\delta + K_{\epsilon}\delta. \quad (A25) $$

Equating like terms leads to the expression

$$ K_{\delta\lambda, \delta\mu, \epsilon, \delta, \gamma} = \sum_{IJ} c_{IJ}^{(\delta\lambda, \delta\mu, \epsilon, \delta, \gamma)} K_{IJ}, \quad (A26) $$

with the $c_{IJ}$ defined in eqs (A2)–(A22). This then defines the coefficients that will need to be updated after each iteration of the non-linear inversion for anisotropic parameters and symmetry orientation.

To determine the $K_{IJ}$ in eq. (A26), we need to define the scattering coefficients $\Omega_{\delta C_{IJ}}$ as in eq. (18), and these explicit expressions are derived and tabulated in the supplementary material.
## APPENDIX B: FIXED ISOTROPIC MODEL

If an iterative approach is desired where the isotropic structure is fixed and we invert for the anisotropic parameters, $\delta\lambda$ and $\delta\mu$ can then be considered only a function of $\epsilon$, $\delta$ and $\gamma$. Given a definition of the isotropic average of an anisotropic structure, this then defines extra scattering terms for $\epsilon$, $\delta$ and $\gamma$.

With an arbitrarily oriented hexagonally symmetric material, we can locally define an orthogonal coordinate system with the $x_3$ axis along the axis of symmetry, as used to define the parameters in eqs (2)–(6). With that choice of local coordinate system, we can define the Voigt average shear and bulk modulus (Babuska & Cara 1991, p. 190) as

$$ \mu^V = \frac{C_{11} + C_{33} - 2C_{13} + 6C_{44} + 5C_{66}}{15} \qquad (B1) $$

$$ \kappa^V = \frac{4C_{11} + C_{33} + 4C_{13} - 4C_{66}}{9}, \qquad (B2) $$

or, equivalently, we can use the Lamé parameters with $\mu^V$ defined as above and

$$ \lambda^V = \kappa^V - \frac{2}{3}\mu^V = \frac{6C_{11} + C_{33} + 8C_{13} - 4C_{44} - 10C_{66}}{15}. \qquad (B3) $$

As in the definitions of terms in eq. (1), this particular coordinate system need only be defined locally, and the expressions for the isotropic average are correct for arbitrary orientation of the symmetry axis.

We can then express the Voigt average $\lambda^V$ and $\mu^V$ in terms of $\lambda$ and $\mu$ in the reference model, as well as the perturbations $\delta\lambda$, $\delta\mu$, $\epsilon$, $\delta$ and $\gamma$ using the definitions in eqs (2)–(6) as

$$ \mu^V = \mu + \delta\mu + \frac{2}{15}(\lambda + 2\mu)(\epsilon - \delta) - \frac{4}{3}\mu\gamma \qquad (B4) $$

$$ \lambda^V = \lambda + \delta\lambda - \frac{6}{5}(\lambda + 2\mu)\epsilon + \frac{8}{15}(\lambda + 2\mu)\delta + \frac{8}{3}\mu\gamma. \qquad (B5) $$

These expressions clearly show that even if we constrain $\delta\lambda$ and $\delta\mu$ to be zero during the inversion, the isotropic velocity structure is still perturbed by non-zero anisotropic parameters. If we fix $\mu^V$ and $\lambda^V$ to be constant during the inversion, this defines $\delta\mu$ and $\delta\lambda$ as

$$ \delta\mu = \frac{4}{3}\mu\gamma - \frac{2}{15}(\lambda + 2\mu)(\epsilon - \delta) \qquad (B6) $$

$$ \delta\lambda = \frac{6}{5}(\lambda + 2\mu)\epsilon - \frac{8}{15}(\lambda + 2\mu)\delta - \frac{8}{3}\mu\gamma. \qquad (B7) $$

If these are substituted into the scattering terms for $\delta\lambda$ and $\delta\mu$, we obtain two additional scattering terms for $\epsilon$, $\delta$ and $\gamma$.

For $\epsilon$, these terms are

$$ \Omega_{\epsilon}^{(\delta\lambda)} = \frac{6}{5}(\lambda + 2\mu)\Omega_{\delta\lambda} \qquad (B8) $$

$$ \Omega_{\epsilon}^{(\delta\mu)} = -\frac{2}{15}(\lambda + 2\mu)\Omega_{\delta\mu}. \qquad (B9) $$

For $\delta$, we have

$$ \Omega_{\delta}^{(\delta\lambda)} = -\frac{8}{15}(\lambda + 2\mu)\Omega_{\delta\lambda} \qquad (B10) $$

$$ \Omega_{\delta}^{(\delta\mu)} = \frac{2}{15}(\lambda + 2\mu)\Omega_{\delta\mu}. \qquad (B11) $$

For $\gamma$, we have

$$ \Omega_{\gamma}^{(\delta\lambda)} = -\frac{8}{3}\mu\Omega_{\delta\lambda} \qquad (B12) $$

$$ \Omega_{\gamma}^{(\delta\mu)} = \frac{4}{3}\mu\Omega_{\delta\mu}. \qquad (B13) $$

## APPENDIX C: VERTICAL SYMMETRY AXIS

It is useful to consider the scattering matrix terms for some end member cases, as these may be desirable in order to set up a linear inversion to obtain a best-fitting starting model prior to a non-linear inversion for orientation of the axis of symmetry and strength of anisotropy.

The easiest case to derive is the vertical axis of symmetry, commonly referred to as Vertical Transverse Isotropy (VTI) or radial anisotropy. In this case, the inner products with the orientation vector $\hat{s}$ are defined such that $s_r$ is 1 and all other inner products are zero.

For $\epsilon$, only the Rayleigh to Rayleigh scattering terms are non-zero, and are defined by

$$ \Omega_{\epsilon}^{(1)} = \dot{U}''\dot{U}' \qquad (C1) $$
$$ \Omega_{\epsilon}^{(2)} = -(\dot{U}'' - \nu''r^{-1}V'')\dot{U}' \qquad (C2) $$

$$ \Omega_{\epsilon}^{(3)} = -\dot{U}''(\dot{U}' - \nu'r^{-1}V') \qquad (C3) $$

$$ \begin{aligned} \Omega_{\epsilon} &= 2\rho\alpha^2(\Omega_{\epsilon}^{(1)} + \Omega_{\epsilon}^{(2)} + \Omega_{\epsilon}^{(3)}) \\ &= 2\rho\alpha^2[r^{-1}(\nu''V''\dot{U}' + \nu'\dot{U}''V') - \dot{U}''\dot{U}'] \end{aligned} \qquad (C4) $$

Likewise for $\delta$ we obtain

$$ \begin{aligned} \Omega_{\delta} &= -\rho\alpha^2(2\Omega_{\epsilon}^{(1)} + \Omega_{\epsilon}^{(2)} + \Omega_{\epsilon}^{(3)}) \\ &= -\rho\alpha^2 r^{-1}(\nu''V''\dot{U}' + \nu'\dot{U}''V') \end{aligned} \qquad (C5) $$

The terms $\Omega_{\gamma}^{(3-6)}$ still include all forms of scattering with a vertical symmetry axis, and so we can write the full scattering contribution from $\gamma$ as

$$ \Omega_{\gamma} = 2\rho\beta^2 \left( 2\Omega_{\epsilon}^{(2)} + 2\Omega_{\epsilon}^{(3)} - \sum_{i=3}^{6} \Omega_{\gamma}^{(i)} \right) \qquad (C6) $$

$$ R \to R := 2\rho\beta^2[2r^{-1}(\nu''V''\dot{U}' + \nu'\dot{U}''V') - 8\dot{U}''\dot{U}' - \cos\eta(\dot{V}''\dot{V}' + \nu''\nu'r^{-2}U''U' + r^{-1}(\nu''U''\dot{V}' + \nu'\dot{V}''U'))] $$

$$ R \to L := -2\rho\beta^2[\sin\eta(\dot{V}''\dot{V}' + \nu'r^{-1}\dot{W}''U')] $$

$$ L \to R := 2\rho\beta^2[\sin\eta(\dot{V}''\dot{W}' + \nu''r^{-1}U''\dot{W}')] $$

$$ L \to L := -2\rho\beta^2 \cos\eta \dot{W}''\dot{W}' \qquad (C7) $$

As in the case of isotropic perturbations, all of the coupling between Love and Rayleigh waves has a $\sin\eta$ dependence and would not be predicted by ray theory along the great-circle path.

## APPENDIX D: HORIZONTAL SYMMETRY AXIS

The full derivations for these terms are available in the supplementary material. For these derivations, it is useful to define some additional angles. $\phi_s$ is the angle defining the orientation of the horizontal symmetry axis, measured counter-clockwise from south. For convenience of expressions, we also define an angle $\zeta = \psi'' + \psi'$. These angles allow us to make simple trigonometric expressions to define the relevant components of the orientation vector $\hat{s}$ to derive the scattering terms.

### D1 Scattering from $\epsilon$

We can combine the derived terms for $\Omega_{\epsilon}$ and group them in terms of the azimuthal dependence. First we can write out the azimuthally independent ($\Omega_{\epsilon}^{(0\phi_s)}$) terms as

$$ \begin{aligned} R \to R := & \rho\alpha^2 \left( \frac{1}{4} \cos 2\eta - \frac{3}{2} \right) \nu''\nu'r^{-2}V''V' + r^{-1}(\nu'\dot{U}''V' + \nu''V''\dot{U}') \\ R \to L := & \frac{1}{4}\rho\alpha^2 \sin 2\eta\nu''\nu'r^{-2}W''V' \\ L \to R := & \frac{1}{4}\rho\alpha^2 \sin 2\eta\nu''\nu'r^{-2}V''W' \\ L \to L := & \frac{1}{4}\rho\alpha^2 \cos 2\eta\nu''\nu'r^{-2}W''W' \end{aligned} \qquad (D1) $$

The terms proportional to $\cos 2\phi_s(\Omega_{\epsilon}^{(\cos 2\phi_s)})$ can be written as

$$ \begin{aligned} R \to R := & \cos 2\phi_s \rho\alpha^2[r^{-1}(\cos 2\psi''V''\dot{U}' + \cos 2\psi'\nu'\dot{U}''V') \\ & -\frac{1}{2}(\cos 2\psi'' + \cos 2\psi')\nu''\nu'r^{-2}V''V'] \\ R \to L := & \cos 2\phi_s \rho\alpha^2 \sin 2\psi''(\nu''r^{-1}W''\dot{U}' - \frac{1}{2}\nu''\nu'r^{-2}W''V') \\ L \to R := & \cos 2\phi_s \rho\alpha^2 \sin 2\psi'\left(\nu'r^{-1}\dot{U}''W' - \frac{1}{2}\nu''\nu'r^{-2}V''W'\right) \\ L \to L := & 0. \end{aligned} \qquad (D2) $$
Likewise for $\Omega_{\epsilon}^{(\sin 2\phi_s)}$, we have

$$R \to R := \sin 2\phi_s \rho \alpha^2 [r^{-1}(\sin 2\psi'' v'' V'' \dot{U}' + \sin 2\psi' v' \dot{U}'' V') \\ -\frac{1}{2}(\sin 2\psi'' + \sin 2\psi') v'' v' r^{-2} V'' V']$$

$$R \to L := \sin 2\phi_s \rho \alpha^2 \cos 2\psi'' \left( \frac{1}{2} v'' v' r^{-2} W'' V' - v'' r^{-1} W'' \dot{U}' \right)$$

$$L \to R := \sin 2\phi_s \rho \alpha^2 \cos 2\psi' \left( \frac{1}{2} v'' v' r^{-2} V'' W' - v' r^{-1} \dot{U}'' W' \right)$$

$$L \to L := 0. \qquad (D3)$$

For $\Omega_{\epsilon}^{(\cos 4\phi_s)}$, we obtain

$$R \to R := \frac{1}{4} \cos 4\phi_s \cos 2\zeta \rho \alpha^2 v'' v' r^{-2} V'' V'$$

$$R \to L := \frac{1}{4} \cos 4\phi_s \sin 2\zeta \rho \alpha^2 v'' v' r^{-2} W'' V'$$

$$L \to R := \frac{1}{4} \cos 4\phi_s \sin 2\zeta \rho \alpha^2 v'' v' r^{-2} V'' W'$$

$$L \to L := -\frac{1}{4} \cos 4\phi_s \cos 2\zeta \rho \alpha^2 v'' v' r^{-2} W'' W'. \qquad (D4)$$

Finally for the $\Omega_{\epsilon}^{(\sin 4\phi_s)}$ terms, we have

$$R \to R := \frac{1}{4} \sin 4\phi_s \sin 2\zeta \rho \alpha^2 v'' v' r^{-2} V'' V'$$

$$R \to L := -\frac{1}{4} \sin 4\phi_s \cos 2\zeta \rho \alpha^2 v'' v' r^{-2} W'' V'$$

$$L \to R := -\frac{1}{4} \sin 4\phi_s \cos 2\zeta \rho \alpha^2 v'' v' r^{-2} V'' W'$$

$$L \to L := -\frac{1}{4} \sin 4\phi_s \sin 2\zeta \rho \alpha^2 v'' v' r^{-2} W'' W'. \qquad (D5)$$

## D2 Scattering from $\delta$

For $\Omega_{\delta}^{(0\phi_s)}$, the azimuthally independent terms, we derive

$$R \to R := \rho \alpha^2 \left( \frac{1}{2} - \frac{1}{4} \cos 2\eta \right) v'' v' r^{-2} V'' V' - \frac{1}{2} r^{-1} (v' \dot{U}'' V' + v'' V'' \dot{U}')$$

$$R \to L := -\frac{1}{4} \rho \alpha^2 \sin 2\eta v'' v' r^{-2} W'' V'$$

$$L \to R := -\frac{1}{4} \rho \alpha^2 \sin 2\eta v'' v' r^{-2} V'' W'$$

$$L \to L := -\frac{1}{4} \rho \alpha^2 \cos 2\eta v'' v' r^{-2} W'' W'. \qquad (D6)$$

Note that these are all proportional but of opposite sign to the $\Omega_{\epsilon}^{(0\phi_s)}$ terms, with the exception of the Rayleigh to Rayleigh scattering term, which has different leading factors.

$\Omega_{\delta}^{(\cos 2\phi_s)}$ can be written as

$$R \to R := -\frac{1}{2} \cos 2\phi_s \rho \alpha^2 r^{-1} (\cos 2\psi'' v'' V'' \dot{U}' + \cos 2\psi' v' \dot{U}'' V')$$

$$R \to L := -\frac{1}{2} \cos 2\phi_s \rho \alpha^2 \sin 2\psi'' v'' r^{-1} W'' \dot{U}'$$

$$L \to R := -\frac{1}{2} \cos 2\phi_s \rho \alpha^2 \sin 2\psi' v' r^{-1} \dot{U}'' W'$$

$$L \to L := 0. \qquad (D7)$$

For $\Omega_{\delta}^{(\sin 2\phi_s)}$, we have

$$R \to R := -\frac{1}{2} \sin 2\phi_s \rho \alpha^2 r^{-1} (\sin 2\psi'' v'' V'' \dot{U}' + \sin 2\psi' v' \dot{U}'' V')$$

$$R \to L := -\frac{1}{2} \sin 2\phi_s \rho \alpha^2 \cos 2\psi'' v'' r^{-1} W'' \dot{U}'$$

$$L \to R := -\frac{1}{2} \sin 2\phi_s \rho \alpha^2 \cos 2\psi' v' r^{-1} \dot{U}'' W'$$

$$L \to L := 0. \qquad (D8)$$
Finally,

$$ \Omega_{\delta}^{(\cos 4\phi_s)} = -\Omega_{\epsilon}^{(\cos 4\phi_s)} \qquad (D9) $$

$$ \Omega_{\delta}^{(\sin 4\phi_s)} = -\Omega_{\epsilon}^{(\sin 4\phi_s)}. \qquad (D10) $$

## D3 Scattering from $\gamma$

For the case of a horizontal symmetry axis for $\Omega_{\gamma}$, we can sum up all terms and once again group by $\phi_s$ dependence. For $\Omega_{\gamma}^{(0\phi_s)}$, the azimuthally independent terms, we obtain

$$ \begin{aligned} R \rightarrow R: &= \rho\beta^2[2r^{-1}(v'\dot{U}''V' + v''\dot{V}''\dot{U}') - \cos\eta r^{-1}(v'\dot{V}''U' + v''U''\dot{V}') - \cos\eta v''v'r^{-2}U''U' \\ &\quad -4v''v'r^{-2}(1 + \cos^2\eta)V''V' - \cos\eta\dot{V}''\dot{V}'] \\ R \rightarrow L: &= -\rho\beta^2[4v''v'r^{-2}\cos\eta\sin\eta W''V' + \sin\eta\dot{W}''\dot{V}' + v'r^{-1}\sin\eta\dot{W}''U'] \\ L \rightarrow R: &= \rho\beta^2[4v''v'r^{-2}\cos\eta\sin\eta V''W' + \sin\eta\dot{V}''\dot{W}' + v''r^{-1}\sin\eta U''\dot{W}'] \\ L \rightarrow L: &= \rho\beta^2[v''v'r^{-2}(\sin^2\eta + \sin\eta - \cos^2\eta - \cos\eta)W''W' - \cos\eta\dot{W}''\dot{W}'] \end{aligned} \qquad (D11) $$

For $\Omega_{\gamma}^{(\cos 2\phi_s)}$, we obtain

$$ \begin{aligned} R \rightarrow R: &= \cos 2\phi_s \rho\beta^2[2r^{-1}(v' \cos 2\psi' \dot{U}''V' + v'' \cos 2\psi'' \dot{V}''\dot{U}') - 4v''v'r^{-2}(\cos 2\psi'' + \cos 2\psi')V''V' \\ &\quad - \cos \zeta(v''v'r^{-2}U''U' + v''r^{-1}U''\dot{V}' + v'r^{-1}\dot{V}''U' + \dot{V}''\dot{V}')] \\ R \rightarrow L: &= \cos 2\phi_s \rho\beta^2[2 \sin 2\psi''v'r^{-1}(W''\dot{U}' - 2v'r^{-1}W''V') - \sin \zeta(v'r^{-1}\dot{W}''U' + \dot{W}''\dot{V}')] \\ L \rightarrow R: &= \cos 2\phi_s \rho\beta^2[2 \sin 2\psi'v'r^{-1}(\dot{U}''W' - 2v''r^{-1}V''W') - \sin \zeta(v''r^{-1}U''\dot{W}' + \dot{V}''\dot{W}')] \\ L \rightarrow L: &= \cos 2\phi_s \rho\beta^2 \cos \zeta \dot{W}''\dot{W}'. \end{aligned} \qquad (D12) $$

Likewise for $\Omega_{\gamma}^{(\sin 2\phi_s)}$, we derive

$$ \begin{aligned} R \rightarrow R: &= \sin 2\phi_s \rho\beta^2[2r^{-1}(v' \sin 2\psi' \dot{U}''V' + v'' \sin 2\psi'' \dot{V}''\dot{U}') - 4v''v'r^{-2}(\sin 2\psi'' + \sin 2\psi')V''V' \\ &\quad - \sin \zeta(v''v'r^{-2}U''U' + v''r^{-1}U''\dot{V}' + v'r^{-1}\dot{V}''U' + \dot{V}''\dot{V}')] \\ R \rightarrow L: &= \sin 2\phi_s \rho\beta^2[2 \cos 2\psi''v'r^{-1}(2v'r^{-1}W''V' - W''\dot{U}') + \cos \zeta(v'r^{-1}\dot{W}''U' + \dot{W}''\dot{V}')] \\ L \rightarrow R: &= \sin 2\phi_s \rho\beta^2[2 \cos 2\psi'v'r^{-1}(2v''r^{-1}V''W' - \dot{U}''W') + \cos \zeta(v''r^{-1}U''\dot{W}' + \dot{V}''\dot{W}')] \\ L \rightarrow L: &= \sin 2\phi_s \rho\beta^2 \sin \zeta \dot{W}''\dot{W}'. \end{aligned} \qquad (D13) $$

For both the $\cos 2\phi_s$ and $\sin 2\phi_s$ terms, the $V''V'$, $V''W'$ and $W''V'$ terms were further simplified using the trigonometric identities

$$ 2 \cos \zeta \cos \eta = \cos 2\psi'' + \cos 2\psi' $$

$$ 2 \cos \zeta \sin \eta = \sin 2\psi'' - \sin 2\psi' $$

$$ 2 \sin \zeta \cos \eta = \sin 2\psi'' + \sin 2\psi' $$

$$ 2 \sin \zeta \sin \eta = \cos 2\psi' - \cos 2\psi'' $$

(D14)

$\Omega_{\gamma}^{(\cos 4\phi_s)}$ and $\Omega_{\gamma}^{(\sin 4\phi_s)}$ are both zero.

## SUPPLEMENTARY MATERIAL

The following supplementary material is available for this article:

**Appendix S1.** The appendix contains derivations of terms in the main body of the paper. It also contains a code fragment used in calculating the scattering coefficient terms used for the example kernels and numerical validation of the main paper. Two additional figures showing the kernels for all 5 elastic parameters are also presented (PDF format).

This material is available as part of the online article from: http://www.blackwell-synergy.com/doi/abs/10.1111/j.1365-246X.2008.03833.x (this link will take you to the article abstract).

Please note: Blackwell Publishing are not responsible for the content or functionality of any supplementary materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.