Testing observables for teleseismic shear-wave splitting inversions: ambiguities of intensities, parameters, and waveforms

We assess the capabilities of different observables for the inversion of core-refracted shear waves (XKS phases) to uniquely resolve the anisotropic structure of the upper mantle. For this purpose, we perform full-waveform calculations for relatively simple, canonical models of upper-mantle anisotropy. The models are characterized by two and four domains of different anisotropic properties. Specifically, we assume hexagonal symmetry with arbitrarily chosen strength of the anisotropy and orientation of the horizontal fast axis. XKS waveforms, generated from plane-wave initial conditions, traverse through anisotropic models and are recorded at the surface by a single station (in case of vertical variations) and by a dense station profile across the laterally and vertically varying structure. In addition to waveforms, we consider the effects of anisotropic variations on apparent splitting parameters and splitting intensity. The results show that, generally, it is not possible to fully resolve the anisotropic parameters of a given model, even if complete waveforms (under noise-free conditions and for the complete azimuthal range) are considered. This is because waveforms for significantly different anisotropic models can be indistinguishable. However, inversions of both waveforms and apparent splitting parameters lead to similar models that exhibit systematic variations of anisotropic parameters. These characteristics may be exploited to better constrain the inversions. The results also show that splitting intensity holds some significant drawbacks: First, even from measurements over a wide range of back-azimuth, there is no characteristic signature that would indicate depth variations of anisotropy. Secondly, identical azimuthal variations of splitting intensity for different anisotropic structures do not imply that the corresponding split waveforms are also similar. Thus, fitting of observed and calculated splitting intensities could lead to anisotropic models that are incompatible with the observed waveforms. We conclude that (band-limited) XKS-splitting inversions and related tomographic schemes, even if based on complete waveforms, are not sufficient to fully resolve the heterogeneous anisotropic structures of the upper mantle and that combinations with alternative methods, based on e.g., receiver-function splitting, P-wave travel-time deviations, or surface waves, are required.


Introduction
Since the early work of Vinnik et al. [1984] and Silver and Chan [1988] there has been considerable progress in the development of inversion approaches to infer the anisotropic properties of the upper mantle from core-refracted phases such as SKS, SKKS or PKS (denoted XKS in the following).The methods are based on linking the anisotropic material properties and their changes with surface observables including waveforms and their proxies, such as (apparent) splitting parameters (e.g., Silver and Savage, 1994] and splitting intensities [e.g., Chevrot, 2000].These approaches are supported by increasingly dense, permanent and temporary, seismic networks that allow for a more precise characterization of waveform variations due to upper-mantle anisotropy [e.g., Ryberg et al., 2005;Mondal and Long, 2020].However, the question arises as to what extent these inversions are able to infer the laterally and vertically varying anisotropic parameters of the upper mantle, and also which type of observable is best suited for the task. Effects of shear-wave splitting have been conventionally described in terms of two splitting parameters: the polarization direction of the fast shear wave with respect to North and the accumulated delay time between the fast and slow shear waves [e.g., Savage, 1999;Liu and Gao, 2013].For a single anisotropic layer, these parameters can be directly related to the anisotropic properties of the medium, where it is generally assumed that the fast polarization is indicative of the (horizontal) orientation of the a-axis of olivine crystals in the upper mantle, whereas the delay time scales with the strength of the anisotropy and/or the extent of the anisotropic domain.In cases of two or more anisotropic layers the two splitting parameters exhibit a characteristic 90°-periodicity as function of initial polarization (back-azimuth) and are only "apparent" in the sense that the obtained values are not directly representative of the anisotropic parameters in the layers [e.g., Silver and Savage, 1994;Rümpker and Silver, 1998;Levin et al., 1999].It is well known that the apparent splitting parameters are non-unique and may show similar variations as functions of back-azimuths even for significantly different anisotropic layering [e.g., Park and Levin, 2002;Abt and Fischer 2008;Latifi et al., 2018].
Splitting intensity is a single parameter that is representative of the relative energy of the XKS-phase on the transverse component of the seismogram (for a given back-azimuth).It was first utilized by Chevrot [2000] in a multi-channel analysis to derive the two splitting parameters of a single anisotropic layer from its variation as function of back-azimuth.The splitting intensity further provides a measure to relate the sensitivity of the split waveforms to anisotropic structures in the mantle [e.g., Favier and Chevrot, 2003].It has also been used in studies of shear-wave splitting tomography [e.g., Long et al. 2008;Sieminski et al., 2008;Chevrot and Monteiller, 2009;Mondal and Long, 2020].
Generally, "full" waveforms provide the most complete information about the anisotropic properties encountered by the wavefield along its way to the receiving station [e.g., Menke and Levin, 2003].However, the limitations of waveform(-based) inversions of XKS splitting to infer complex anisotropic structures in the mantle have not been investigated in detail.While finite-frequency effects are considered by the tomographic schemes that relate splitting intensity and splitting sensitivity to anisotropic variations in the mantle, their resolving power is certainly limited.
For example, Chevrot [2006] showed that vectorial tomography based on splitting intensity is able to resolve fast-axis variations if the strength of anisotropy is assumed to be known a priori.More recently, Mondal and Long [2019], concluded that their tomographic approach is robust when constraining one of the parameters (strength of anisotropy or symmetry axis) and by conducting the inversion for anomalies of the other quantity.
To assess the suitability of the different observables and their effects on shear-wave splitting inversions, the finite-frequency nature of XKS waveforms must be considered.A more direct approach to examine the resolving power of XKS splitting observations is to consider the full wavefield generated for different anisotropic models.
Previous full-waveform forward modeling and observations showed that vertical and lateral variations of anisotropy can leave a complex and non-unique signature on the XKS wavefield traversing the upper mantle [e.g., Rümpker et al., 2003;Zhao et al., 2008;Kaviani et al., 2011;Wölbern et al., 2014;Reiss et al., 2018].Which observables are best suited to unravel the waveform modifications is subject to the investigations presented in this study.We employ a "brute-force" inversion approach by sampling the complete model space and by directly comparing full waveforms and waveform proxies.To this end, we consider different canonical models of anisotropic structure in the upper mantle and focus on accessing the suitability of apparent splitting parameters, splitting intensity, and full waveforms for the inversions to infer the underlying model parameters.

Vertically varying anisotropic medium
Here, we first assume a single anisotropic layer with a horizontal symmetry axis (oriented at an angle , with respect to North, Figure 1a) where a vertically incident shear wave is split into fast and slow components, polarized horizontally, perpendicular to each other.Generally, the two quasi-shear waves travel at different speeds such that they become separated by the delay time, , over the thickness of the layer.For teleseismic phases in a radially symmetric Earth, the polarization of the incident wave, , (unit vector) is given by the back-azimuth, , measured with respect to North and defined by the azimuth of the earthquake with respect to the receiving station (see Figure 1a).Mathematically, the effect can be described by the following system of equations.Similar formulations have been presented before [e.g., Silver and Savage, 1994;Rümpker and Silver, 1998;Montagner et al., 2000;Romanowicz and Yuan, 2012].Here we use a modified self-consistent formulation that can be easily adopted and generalized to multi-layered media.We consider a vertically propagating, horizontally polarized shear wave incident from below.
With the waveform in the frequency domain, , the displacement can be given as . (1) Georg Rümpker et al. 4 In view of Figure 1a, the initial polarization may be expressed in terms of its geographic (east, north) components .
(2) Thus, we may write where the geographic components of the displacement are given by , .
Upon entering the anisotropic layer, the displacement is split into a fast and a slow component.The corresponding unit vectors are given by (see Figure 1a) or, in terms of the displacement where  corresponds to horizontal symmetry axis (or fast axis) in the layer.After passing through the anisotropic layer, fast and slow waves are offset from each other by the delay time, .The offset may be distributed symmetrically, which leads to the displacement components (in fast-slow coordinates) .
In view of eq. ( 5), the corresponding displacement in geographic coordinates is given by To summarize, we express the shear-wave displacement at the top of the anisotropic layer as where  and  are matrices defined by the splitting parameters  (i.e., the fast-polarization direction or fast axis) and  (delay time), respectively, of the anisotropic layer, and the product matrix  represents the corresponding splitting operator.
The generalization of eq. ( 10) to account for the effects of a second anisotropic layer, is straightforward where the indices 1 and 2 refer to the anisotropic parameters of the first and second anisotropic layer, respectively.
Further details and the generalization to multiple anisotropic layers are discussed in, e.g., Silver and Savage [1994] and Rümpker and Silver [1998].
Shear-wave splitting may also be calculated directly in the time domain.In this case the diagonal matrix  can be replaced by a time-shift operator that accounts for the advance and delay of the time-domain waveform components expressed in fast-slow coordinates.However, the generalization to multiple layers is more cumbersome as the multiplication (as applied in eq.11) of the time-domain operators is not defined and the operations can only be performed successively.
One may note that the formalism described above is "complete" in the sense that there is no approximation with respect to the frequency content of the vertically propagating wavefield in a laterally homogenous medium.
For multiple anisotropic layers, effects of internal reverberations are not considered.However, this limitation is not significant if the overall (isotropic background) velocity remains near constant as it is the case within the upper mantle.Effects of lateral variations of material properties are considered in the next section.

Laterally and vertically varying medium
The wave equation for anisotropic elastic media has not been discussed explicitly in the previous section, it nevertheless forms the basis for the derivation of the presented equations.It may be written as [Aki and Richards, 1980] (12) with displacement   , stress tensor   , and density .The comma denotes the partial derivative and indices ,  range from 1 to 3, and the summation convention applies.In the following, we consider 2-dimensional variations of elastic parameters (in the  1 , 3 -plane, where  3 points vertically downwards and  1 along the direction of material changes) and further assume that all derivatives of the material properties as well as of the displacement field with respect to  2 vanish (i.e. 2 = 0).
The full wave equation for the three displacement components can then be written in the form [Ryberg et al., 2002;Zhao et al., 2008] , , , where we introduce the crystallographic notation by replacing the indices as follows, , , The 6 components of the stress tensor,   , are then given by .( 15) Georg Rümpker et al. 6 The elastic constants,   , may be given in matrix form , ( 16) where we assume an elastic tensor of orthorhombic symmetry with, at this point, symmetry axes (a,b,c) oriented parallel to the coordinate axes.The matrix is symmetric and elements  2 do not contribute because of the reduced dimensionality.Rotations of the 4 th order elastic tensor can be performed by the transformation We now assume the a and c axes in the horizontal  1 , 2 -plane and only account for rotations of the elastic tensor with respect to the vertical  3 axis, such that the rotation matrix,   , can be given by ( 18) with rotation angle .In view of this, the application of eq. ( 17) yields an elastic tensor with a general orientation of the orthogonal a-c axes in the horizontal plane, such that its matrix representation takes the form The elastic constants of the vertically and laterally varying model are more specifically defined in appendix A2.
The displacement components are obtained by solving the wave equation (eq.13) numerically using a finite difference method.The corresponding synthetic seismograms are complete in the sense that they include finitefrequency effects, scattering and diffractions due to material inhomogeneities, and conversions between wave types.As indicated in the formulation given above, the numerical simulations are performed with respect to a 2D cartesian coordinate system.Changes in material properties and wavefield variations are restricted to the  1 -  3 plane.For the numerical evaluation, we adopted the explicit second-order isotropic finite-difference formulation of Kelly et al. [1976] and generalized it to the anisotropic case.Free-surface conditions are specified at the upper boundary of the computational domain.Numerical grid dispersion is limited by choosing sufficiently small step sizes in space and time.The computational domain is chosen large enough to suppress any influence of spurious reflections originating from the grid boundaries.Effects of wavefront curvature are not considered.

Splitting parameters and splitting intensity
To access the suitability of different observables for anisotropic inversion schemes we analyze the effects of seismic anisotropy not only by waveforms, as calculated by the schemes described above, but also in terms of

Observables for XKS-splitting inversion
splitting parameters [e.g., Silver and Savage, 1994] and splitting intensity [Chevrot, 2000].We will investigate the ability of these "secondary" characteristics or waveform proxies to discriminate between different anisotropic structures and to resolve the anisotropic properties within Earth's upper mantle.
The splitting parameters for a single anisotropic layer are simply given by the fast axis and delay time of that layer.They can be obtained from the horizontal component waveforms by multiplication with the inverse splitting matrix, , that best linearizes the resulting waveform (see eq. 10).However, in cases of two or more anisotropic layers the linearization may lead to an incident waveform that is not aligned with the back-azimuth, as it would be required for XKS phases in a radially symmetric Earth (before entering the anisotropic domain).It is therefore convenient to first express the waveform in terms of radial-transverse coordinates and to obtain the splitting parameters from the inverse of that best minimizes the transverse waveform component.Explicit expression for are given in the appendix (A1).
In vertically varying media, the splitting parameters obtained by this procedure are only "apparent" as the contributing individual layer parameters cannot be resolved (in view of the relatively long periods of XKS phases).
Approximated closed form solutions for 2-layer apparent splitting parameters (appendix A1) exhibit the well-known -periodicity as functions of back-azimuth (Silver and Savage, 1994).This periodicity also applies to general multi-layered anisotropic media [Rümpker and Silver, 1998].
The splitting intensity for XKS phases relates to the relative amplitude of the transverse component [Chevrot, 2000].It depends on the angle between the back azimuth and the fast axis and on the delay time .Alternatively, the splitting intensity can be measured by projecting the transverse-component waveform on the time-derivative of the radial-component waveform [Monteiller and Chevrot, 2010] , where  denotes time and the integration limits are chosen to enclose the complete XKS phase.Note that  is given in units of time, usually in seconds.For a single anisotropic layer, the splitting parameters can be determined from the azimuthal dependence of the splitting intensity.However, it is well known, that, in contrast to splitting parameters, the azimuthal variations of the splitting intensity are not indicative of vertical variations of anisotropy, e.g., due to two or more anisotropic layers [Kong et al., 2015].
For a single anisotropic layer, analytical expressions for the splitting intensity have been derived which are valid at relatively low frequencies [Chevrot, 2000; see also appendix A1]. (21) This expression can be generalized to two or more anisotropic layers by superposition of the splitting intensities for the individual layers.By using the approach of Rümpker et al. (2014) for the derivation of closed-form expressions of the Ps-phase moveout due to multiple anisotropic layers above the converting interface, we find here that the 2-layer splitting intensity can be expressed as (22) with "effective" layer parameters , .
From eq. ( 22) it is obvious that the 2-layer splitting intensity is not indicative of layered anisotropy; furthermore, it is commutative, i.e., it is independent of the order of the layers.This contrasts with the apparent splitting parameters and with the corresponding waveforms, as we will see in the examples given below.

Reference model
We consider a layered medium consisting of two anisotropic layers defined by four parameters.We choose  1 = 80°,  1 = 0,5 s, for the lower layer and  2 = 10°,  2 = 0,8 s, for the upper layer (Figure 1b).Fast axes for this medium differ by 70° and delay times are relatively large such that both layers are expected to have a significant influence on the traversing shear wavefront and the resulting splitting.This is also shown by the corresponding transverse-component waveforms which exhibit characteristic variations as a function of back-azimuth, whereas the radial components are almost unchanged (Figure 2).For the modeling, we assume a near-sinusoidal incident waveform with a dominant period of 9 s (Figure S1), which is considered typical for XKS waves.

Apparent splitting parameters
Based on eqs.(A7, A8) we calculate apparent splitting parameters for the reference model at 36 equally distributed values of back-azimuth between 0° and 350°.The values represent rather ideal observables and exhibit the characteristic 90°-periodicity for layered anisotropic media.Discontinuities of apparent fast polarization occur near back-azimuths of 25°, 115° etc. and are accompanied by increasing apparent delay times.The corresponding transverse-component waveforms near these back-azimuths are characterized by vanishingly small amplitudes (see Figure 2).
To investigate the ambiguities of the apparent splitting parameters with respect to anisotropic model parameters, we calculate apparent splitting parameters for a complete range of different 2-layer models (also based on eqs.A7, A8).We assume and at intervals of 5° and 0.1 s, respectively.This results in 518.400 different test models, for which the apparent splitting parameters are calculated.We compare the apparent

Observables for XKS-splitting inversion
splitting parameters for these models with those obtained for the reference model described above by evaluating their squared and summed differences.We account for the different units of the 2 splitting parameters by multiplication with appropriate weights (1/30°) and (1/0.5 s).It is found that the 20 best models provide a reasonable spread in model parameters while the corresponding observables are nearly indistinguishable (in view of error estimates for observation based on real data).This approach can be considered a "brute-force" inversion by forward modeling and data fitting.Results for the 20 best-fitting models are shown in Figure 3a.Whereas the apparent splitting parameters as functions of back-azimuth are remarkably similar, the corresponding models (Figure 3b) vary significantly in both fast axes and delay times.For the lower layer, possible delay times can differ by up to 0.9 s (between 0.5 s and 1.4) and fast axes by up to 40° (between 40° and 105°) from the corresponding parameters of the reference model.In the upper layer the possible variations for the fast axis are even larger with differences (w.r.t. the reference model) of up to 50°, and 1.2 s for the delay time (Figure 3b).Following the possible values of the fast axis in the lower layer from 40° to 105°, we find that the corresponding fast axis in the upper layer varies from 140° over 180° to 20°.Thus, the inversion yields a wide range of possible solutions for the fast axes in the upper and lower layer with a systematic trade-off, whereas their difference can vary by up to 20° (between 65° and 85°).

Splitting intensity
Based on the numerical approach (eq.20), we further calculate splitting intensities for the 518.400 anisotropic models and again select the 20 models for which the results best fit the splitting intensities for the reference model.
The resulting splitting intensities are indistinguishable and may be considered equivalent (Figure 4a), especially, in view of a real-data analysis.For comparison, we also show splitting intensities calculated using the new analytical expressions (eq.23).The results are in good agreement, but it should be noted that individual delay times in the reference model are relatively small compared to the dominant period, which facilitates the approximation.It is interesting to note that the splitting intensities are not in any way indicative of depth variations of the anisotropic structure beneath the profile, in agreement with Kong et al. [2015].
While the splitting intensities are very similar, the corresponding models significantly disagree (Figure 4b), even more than the previous models derived from the comparison of the apparent splitting parameters (see Figure 3b).
Possible models cover a wide range of fast axes in both the upper and lower layer while delay times range from 0.1 s

Georg Rümpker et al.
to up to 1.4 s.Also, in difference to the previous models, there is no obvious systematic tradeoff between fast axes and delay times.However, models with an upper-layer or lower-layer fast axis of 25° to 35° exhibit either small delay times or near-perpendicular fast axes in the other layer.These parameters can be explained by the effective one-layer parameters for the 2-layer splitting intensity as function of back-azimuth.For the reference model, the value for the "effective" fast axis,  1,2 , is 29° and the corresponding "effective" delay time is  1,2 = 0.53 s (see eq. 23).

Waveforms
Full ("finite frequency") waveforms are thought to provide the most complete information about layered anisotropy.We use the transverse-component waveforms to obtain the 20 best-fitting anisotropic models (of 518.400 tested) since, for relatively long periods ( ), the radial components remain almost unaffected by the anisotropic layers.Figure 5 shows that the transverse-component waveforms for the 20 different models are almost identical and agree well with those of the reference model.The corresponding fast axes and delay times are quite similar to those obtained by fitting of the apparent splitting parameters (see Figure 3b).Fast axes in the lower layer range from 40° to 105° with delay times between 0.4 s and 1.2 s.Fast axes in the upper layer vary from 130° over 175° to 20° with delay times between 0.4 s and 1.6 s.
It is instructive to compare the transverse-component waveforms for these models (Figure 5b) with those obtained from the apparent splitting parameters (models shown in Figure 3b) and from splitting intensities (Figure 4b).
We find that the latter are more variable and partially incompatible with those of the reference model (Figure 6).
These are also different from the models obtained from the fit of apparent splitting parameters.

Model setup
In addition to the previous setup, we now include lateral variations of anisotropy by adding a second layered structure to the reference model (Figure 7).The new model consists of 4 anisotropic domains, where the anisotropic

Observables for XKS-splitting inversion
parameters on the left-hand side are identical to those of the reference model.On the right-hand side we choose parameters for which we expect that the splitting is nearly identical.From the previous comparison of best-fitting waveforms and splitting parameters we select a model with  1 = 50°,  1 = 0.7 s, for the lower layer and  2 = -15°, (or equivalently 165°),  2 = 0.5 s, for the upper layer.We choose this model to show that different anisotropic structures can lead to remarkably similar waveform and splitting effects, even if significant lateral and vertical variations occur.
Conversely, this would suggest that the inversion of the waveforms or splitting parameters is ambiguous with respect to both lateral and vertical changes of anisotropic parameters in the mantle.
For the waveform calculations we revert to eq. ( 13) which also requires to explicitly set up elastic constants which are consistent with the anisotropic parameters in the four domains.In this case, we chose elastic constants for transverse isotropic (TI) media (appendix A2) which allows for a straightforward scaling between delay times and the corresponding percentage of anisotropy.The anisotropic domains are 100 km thick each and extend

Georg Rümpker et al.
from 35 km down to 235 km depth.Laterally, the blocks extend over the complete computational domain.They are embedded between an isotropic crust above and an isotropic mantle below.To generate the waveforms, we define a plane horizontal XKS wavefront entering the region from below and the pulse shape of the incident wave is given by a near-sinusoidal function.In view of the anisotropic parameters and to account for splitting in all domains, we set the initial polarization to 90° (parallel to the  1 axis).The orientation of the station profile also follows this axis.The lateral boundary between the two anisotropic domains occurs at the (arbitrarily chosen) range of 950 km.

Wavefronts and waveforms
In Figure 8 we show snapshots of the  1 and  2 displacement fields at three different times ( in the laterally varying model.The incoming wavefront is initially polarized along  1 , which corresponds to a back-azimuth of 90°.In this specific case, the two components can be interpreted as radial and transverse components also. On entering the anisotropic domains, part of the energy is transferred from  1 to  2 .The wavefront is differently affected by the lateral changes in anisotropy.However, after the wavefront has passed through the entire mantle, the wavefronts show little lateral variations, except for some minor diffraction effects induced by the lateral changes in anisotropy.Also visible, near the bottom of the cross-section, are faint reflections from the horizontal boundary between the two layers.The corresponding radial and transverse waveform components (Figure 9) are nearly identical for all stations along the profile.This is to be expected and is consistent with the findings from the vertically varying models presented in the previous section.
We further investigate the influence of the initial polarization on the waveform similarity.Radial and transverse component XKS waveforms for 18 evenly distributed back-azimuths between 0° and 170° and recorded at distance ranges between 660 km and 1240 km (at intervals of 20 km) are shown in Figure 10.The radial component waveforms at the different positions are almost indistinguishable and the transverse components exhibit only minor variations, which is due to diffraction effects near the center of the profile.In comparison to the vertically varying modeling, the transverse waveforms are slightly more and dispersed.This is inherent to the more complete finite- When passing through the anisotropic domains, part of the energy is transferred from  1 to  2 .In this specific case, the two components correspond to the radial and the transverse components.Range (horizontal axis) and depth (vertical axis) are given in km.Lateral and horizonal lines denote boundaries between anisotropic domains in the mantle (see Figure 7).Georg Rümpker et al.
difference modeling approach, which also accounts for reflections and scattering from the internal boundaries and inhomogeneities.

Splitting parameters
Directly from the waveforms we also determine (numerically calculated) splitting parameters based on minimizing the transverse component energy [e.g., Link et al., 2022].The results shown in Figure 11 agree well with those for the 2-layer reference model (see Figure 3a) which have been obtained from analytical expressions for the apparent splitting parameters (eq.A7-8).Waveform diffraction, as visible in Fig. 9b, plays a role and affects some of the splitting parameters obtained from waveforms at the transition between the anisotropic domains.However, the results confirm that the splitting along the profile is almost uniform, especially in view of real-data observations where noise and other sources of uncertainty (e.g., sensor orientation) can have a considerable influence.Observables for XKS-splitting inversion 15

Splitting intensity
The previous findings are further supported by splitting intensities as function of back-azimuth which we obtain by repeating the finite-difference waveform calculations for 37 values of initial polarization from 0° and 360° (Figure 12).The splitting intensities at different receiver positions show minor variation even for stations located close to the center of the profile, where the diffraction effects are also more pronounced.Therefore, the superimposed curves from the separate locations can be approximated well by the single curve based on effective layer parameters,  = - 0.54 s and  = 27.3°,obtained from the analytical expression (eqs.22, 23) using the parameters on the right-hand side of the model shown in Figure 7. Again, in view of real data, the variations along the profile as observed in this numerical example, are probably not significant.

Discussion
The results presented above show that it is not possible to fully resolve the vertical and lateral variations of anisotropic parameters in the mantle from shear-wave splitting observations of teleseismic phases.This is certainly true if splitting intensities are used to quantify the splitting, but it also holds for apparent splitting parameters and even full waveforms.However, the latter two observables exhibit systematic variations between parameters of the upper and lower layers that can be exploited to better constrain the anisotropic properties.It is clear from Figs. 3 and 5 that by fixing the fast axis in one of the layers (based on, e.g., constraints from absolute plate motion) the three remaining parameters can be obtained without further ambiguity.
Abrupt lateral and vertical changes, as studied here, may not occur in nature.On the other hand, the effective Fresnel zones of XKS waves are relatively large such that effects of sharp (lateral) contrasts on the waveforms are smoothed out in any case.Another aspect is the purely vertical incidence of the XKS wavefronts considered in our modeling.In reality, the angles of incidence are close to about 15° and this could potentially help in resolving heterogenous anisotropic structure [Huang and Chevrot, 2021;Löberich and Bokelmann, 2022].Tilted-axis anisotropy can have a significant influence on the waveforms [e.g., Šílený and Plomerová, 1996].In our examples, however, we deliberately chose a relatively simple type of the anisotropy which can be described by only two parameters.Additional unknowns, such as tilted axes, would further complicate the inversion and probably lead to an even wider range of possible models.
The question may arise if higher-frequency waveforms could be used to better discriminate between the vertically and laterally varying models.We tested this by calculating transverse waveforms, characterized by dominant periods of 4 s, for the 20 best-fitting models obtained from the fitting of longer-period waveforms (see Figure 5).The results (Figure S2) show that differences become apparent but may not be significant enough for a reliable discrimination (also considering the longer period nature of XKS phases).However, higher frequency SKS waveforms have been used previously to discriminate between complex models of anisotropy [Menke and Levin, 2003;Rümpker et al., 2003] and we will further explore their potential in future studies.
In view of the results for laterally and vertically varying model we may consider different inversion strategies.
We exclude splitting intensities since they do not provide any insight into depth variations and may lead to models not compatible with the observed waveforms.If we observed splitting intensities such as those given in Figure 12, we could only assume a uniform layer without any variations of anisotropy.The inversion in terms of this model leads to a fast axis of 27° and a delay time of about 0.6 s, values that correspond to the effective layer parameters given in eq. ( 23), but without any overlap with the real 8 anisotropic parameters of the model (see Figure S3).Our findings regarding the limitations of splitting intensity to infer complex anisotropy are in agreement with a study by Kong et al. [2015] where differences between splitting parameters and splitting intensity are analyzed in detail.
Measurements of apparent splitting parameters (Figure 11) can provide clear evidence for depth variations of anisotropy.However, in our case, the lateral variations are not obvious.From the observation one would, therefore, assume a model with vertical variations only and perform the inversion accordingly.To simulate this situation, we perform inversion for the averaged curves of the apparent splitting parameters obtained by stacking the splitting parameters for all stations along the profile.Anisotropic model parameters obtained from the 20 best fitting curves (based on analytical expressions) are given in Figure S4.Not surprisingly, the parameters for the individual models are comparable to those already shown in Figure 3.
In a combined inversion approach, azimuthal variations of apparent splitting parameters can provide first evidence for depth variations of anisotropy such that the following waveform inversion focuses on obtaining model parameters of (at least) 2 anisotropic layers.However, one may wonder about the possibility to directly discriminate between single and two-layer anisotropy using waveforms only.To investigate this, we first apply an inversion of the split waveforms that accounts for a single anisotropic layer.This is done by rotating the recorded waveforms into a trial fast-slow coordinate system and by subsequent time shifting of ± /2 with the objective to minimize the transverse component energy.
Results of applying the approach to the waveforms in Figure 10 are shown in Figure 13.As before, we ignore the vanishingly small lateral variations and use azimuthally varying waveforms as input for the inversion obtained from averaging over all receiver locations.We note that the minimization is certainly not sufficient.The corresponding single-layer models are all grouped near a fast axis of 25° and delay time of 0.6 s.Interestingly, the results are similar to those obtained from the inversion of the splitting intensities (Figure S3).However, in the case of using waveforms, the remaining significant transverse energy provides clear evidence for the influence of anisotropic layering.This essential information is lost when using splitting intensity as an observable.

Observables for XKS-splitting inversion
In a more general approach, we apply a 2-layer splitting inversion to the waveforms (by applying a second rotation and time shift).In turn, this leads to a significant reduction of the transverse amplitudes for the 20 models that best minimize the energy (Figure 14).The model parameters are again similar to those obtained from the inversion of apparent splitting parameters, but 2 of the models fully agree with the input models used in the finitedifference modeling.This is certainly due to using the complete waveforms as observables.
One may also be tempted to use the radial-component waveform as a proxy for the unsplit incoming waveform and to minimize the differences between calculated and observed transverse waveforms in the inversion.We found, however, that the energy minimization is more stable as it is applicable even if the radial components are significantly modified due to the splitting.

Conclusions
We have applied different approaches of waveform modeling to investigate effects of heterogeneous anisotropic structures on shear wave splitting and to assess the resolving power of different observables such as waveforms, apparent splitting parameters, and splitting intensities.For 2-layer anisotropic media we confirm earlier inferences which have shown that azimuthal variations of the (two) apparent splitting parameters are ambiguous and do not allow to uniquely identify the four anisotropic parameters of the underlying structure.Relevant examples for the limitations of inversions based on splitting parameters are also discussed in a recent study by Lamarque and Agostinetti [2020].Ambiguities are even more significant when the splitting intensity is used to quantify the azimuthal waveform variations.While the (single parameter) splitting intensity may seem a convenient measure for characterizing the waveform splitting, it also holds significant drawbacks: First, even from measurements over a wide range of back-azimuth, there is no characteristic signature that would indicate depth variations of anisotropy [see Kong et al., 2015].Secondly, identical azimuthal variations of splitting intensity for different anisotropic structures do not imply that the corresponding split waveforms are also similar.This can have profound consequences for any inversion scheme that is based on fitting of observed and calculated splitting intensities, as it could lead to anisotropic models that are not compatible with the observed waveforms.
However, as our modeling shows, even by employing the complete ("full") XKS waveforms in the inversion we cannot remedy the general problem of non-uniqueness.Waveforms for significantly different anisotropic models, often, are indistinguishable.This remains true even under ideal (and noise-free) conditions with waveforms from the complete range of back-azimuths.The problem is partly related to the band-limited nature of XKS waveforms with relatively long dominant periods that mask the distinct shear-wave arrivals from layered structures.Therefore, higher-frequency waveform information should be included in the analysis whenever possible.Our results also show that inversions of both waveforms and apparent splitting parameters lead to similar models that exhibit systematic variations of anisotropic parameters.These characteristics can be exploited to better constrain the anisotropic structures.
We conclude that complete waveforms or splitting parameters are more appropriate observables for XKSsplitting inversions and related tomographic schemes.In view of our results, the use of splitting intensities seems questionable due to their insensitivity to vertically varying anisotropic structures.However, (band-limited) XKSsplitting inversions, alone, are not sufficient to resolve general heterogeneous anisotropic structures of the upper mantle and, therefore, combinations of methods based on, e.g., receiver-function splitting [e.g., Park and Levin, 2016;Link et al., 2020]; P-wave travel-time deviations (e.g., Babuška and Plomerová, 2020), or surface waves [e.g., Endrun et al., 2011] are required.Additionally, geodynamic constraints on mantle flow [e.g., Conrad et al., 2007] can also provide essential information to better resolve the anisotropic structures.

Figure 1 .
Figure 1.(a) Relationships between coordinates and angles used to describe the shear-wave splitting.The initial polarization points towards the source.and denote geographic North and East directions.and correspond to the fast and slow axes in the anisotropic layer.For the radial and transverse directions, and , respectively, we follow the seismological conventions (e.g., Gubbins, 1990).Angles  and  denote back-azimuth and fast axis with respect to North.(b) Schematic view of the 1-dimensional 2-layer reference model with anisotropic parameters of the lower and upper layers (fast axes with respect to North and delay times between the fast and the slow shear-waves).The inverted triangle denotes the seismic station.Note, that for long-period XKS phases, as considered in this study, individual split shear waves are not fully separated.

Figure 2 .
Figure 2. (a) Radial and (b) transverse-component XKS waveforms as functions of back-azimuth after the shear wave has passed through the 2-layer anisotropic reference model (Fig. 1b).Waveforms in (b) are enlarged a factor 1.5 with respect to waveforms in (a).

Figure 3 .
Figure 3. (a) 20 best-fitting apparent splitting parameters as function of back-azimuth.The green circles mark the apparent splitting parameters of the reference model at intervals of 10°.(b) Corresponding models of 2-layer anisotropy.Lower and upper-layer parameters for one model are marked by the same (colored) symbol.Parameters of the reference model are denoted by the red (filled) and black circles.The fit deceases from red, orange/yellow to green and finally blue, while symbols (from circle, square, triangle to inverted triangle) are repeated.Note that the red circle is overprinted by a (light) blue triangle in the upper panel.

Figure 4 .
Figure 4. (a) 20 best-fitting splitting intensities as function of back-azimuth.The green circles mark the splitting intensities of the reference model at intervals of 10°.The pink crosses indicate splitting intensities obtained from the closed-form (approximate) expressions for the effective anisotropic 1-layer parameters given by eq.(23) with  1,2 = 29° and  1,2 = 0.53 s.(b) Corresponding models of 2-layer anisotropy.Lower and upper-layer parameters for one model are marked by the same (colored) symbol.Parameters of the reference model are denoted by the red (filled) and black circles.

Figure 5 .Figure 6 .
Figure 5. (a) 20 best-fitting transverse-component waveforms as function of back-azimuth.The green circles correspond to values for the waveforms of the reference model.(b) Corresponding models of 2-layer anisotropy.Lower and upper-layer parameters for one model are marked by the same (colored) symbol.Parameters of the reference model are denoted by the red (filled) and black circles.

Figure 7 .
Figure 7. Laterally and vertically varying anisotropic reference model characterized by 4 distinct anisotropic domains in the depth range between 35 km and 235 km.In the numerical model, corresponds to  1 and to  2 .

Figure 8 .
Figure 8. Snapshots of the  1 (a) and  2 (b) displacement components at three times ( 1 <  2 <  3 ) in the laterally varying model.The incoming wavefront is initially polarized along  1 only, corresponding to a back-azimuth of 90°.

Figure 9 .
Figure 9.The (a) radial and (b) transverse-component XKS waveforms recorded at the receivers on the surface of the laterally varying model.The back-azimuth is 90° as in Fig. 8.An amplification factor of 3 is applied to the transverse components relative to the radial components.Note that the vertical boundary between the anisotropic domains occurs at range 950 km along the station profile.

Figure 10 .Figure 11 .
Figure 10.The (a) radial and (b) transverse-component XKS waveforms as functions of back-azimuth and for all stations along the profile (superimposed).An amplification factor of 3 is applied to the transverse components.

Figure 12 .
Figure 12.Splitting intensities as functions of back-azimuth for all stations along the profile (lower panel).Superimposed splitting intensities for all stations (upper panel).Circles denote analytical results according to , (eq.22), with effective 1-layer parameters obtained from eq. (23).

Figure 13 .
Figure 13.(a) Transverse-component waveforms for the 20 best models obtained after application of 1-layer splitting to minimize the energy of the transverse waveforms.(b) The corresponding 20 best 1-layer model parameters (colored symbols).Large open circles and squares correspond to parameters of the four anisotropic domains (see Fig. 7).

Figure 14 .
Figure 14.(a) Transverse-component waveforms for the 20 best models obtained after application of 2-layer splitting to minimize the energy of the transverse waveforms.(b) The corresponding 20 best 2-layer model parameters (identical colored symbols belong to the same vertically varying model).Large open circles and squares correspond to parameters of the four anisotropic domains (see Fig. 7).