The expression of mantle seismic anisotropy in the global seismic w av eﬁeld

results conﬁrm that the seismic phases that are commonly used in splitting techniques are indeed strongly inﬂuenced by mantle anisotropy. Ho wever , we also identify less commonly used phases whose waveforms reﬂect the effects of anisotropy. For example, PS is strongly affected by splitting due to seismic anisotropy in the upper mantle. We show that PS can be used to ﬁll in gaps in global coverage in shear-wave splitting data sets (for example, beneath ocean basins). We ﬁnd that PcS is also a promising phase, and present a proof-of-concept example of PcS splitting analysis across the contiguous United States using an array processing approach. Because PcS is recorded at much shorter distances than ∗ KS phases, PcS splitting can therefore ﬁll in gaps in backazimuthal cov erage. Our wav eﬁeld differencing results further hint at additional potential novel methods to detect and characterize splitting due to mantle seismic anisotropy.

to detect and characterize it.Anisotropy in Earth's crust and upper mantle can be measured using receiver function analysis (e.g.Levin & Park 1997 ;Schulte-Pelkum et al. 2005 ;Nikulin et al. 2009 ;Haws et al. 2023 ) or surface wave tomography (e.g.Panning & Nolet 2008 ;Ferreira et al. 2010 ;Zhu et al. 2020 ).Waveform inversion techniques have also been developed to characterize radial anisotropy in the mantle (e.g.Kawai & Geller 2010 ;Suzuki et al. 2021 ).Probably the most commonly used method to study seismic anisotrop y in Earth' s mantle, and the one we focus on in this paper, invokes measurements of so-called shear-wave splitting (e.g.Long & Silver 2009 ).The distribution of shear wave energy away from its initial polarization direction, splitting the wave into two quasishear waves, is indicative of seismic anisotropy (e.g.Silver & Chan 1991 ).
Depending on which portion of the mantle is being studied, different shear phases (or combinations of phases) are typically used.Fig. 1 (a) summarizes most commonly used seismic phases for the analysis of mantle anisotropy using shear-wave splitting.To analyze seismic anisotropy in the upper mantle directly beneath the receiver, the most commonly used phase is SKS, sometimes supplemented with SKKS and PKS (e.g.Silver & Chan 1991 ;Chevrot 2000 ;Liu et al. 2014 ;Walpole et al. 2014 ;Graw & Hansen 2017 ;Lopes et al. 2020 ).For such an analysis of upper mantle anisotropy, it is often assumed that the influence of lowermost mantle anisotropy is ne gligible.An alternativ e and commonly applied approach to studying upper mantle anisotropy is to infer source-side anisotropy from splitting of teleseismic S waves using explicit receiver-side anisotropy corrections (e.g.Russo & Silver 1994 ;Lynner & Long 2013 ;Walpole et al. 2017 ;Eakin et al. 2018 ).Source-side direct S splitting can also be used to study transition zone and uppermost lower mantle anisotropy in places where deep earthquakes occur (e.g.Foley & Long 2011 ;Mohiuddin et al. 2015 ).These observational strategies are well-established and have been used to map upper mantle anisotropy across much of Earth's landmasses.Beneath the oceans, ho wever , shear-wave splitting constraints on upper mantle anisotropy are sparse (e.g.IRIS DMC 2012 ), due to the paucity of seismic receivers.
Lowermost mantle anisotropy is generally more challenging to measure than seismic anisotropy in the upper mantle.A major reason for this is that all seismic waves that may sample seismic anisotropy in the lowermost mantle also travel through the upper mantle, potentially accumulating an upper mantle splitting signal before they reach the station (Fig. 1 a).Different measurement strategies have been developed to distinguish between an upper and lowermost mantle contribution.Such techniques include the analysis of differential S-ScS splitting (e.g.Wookey et al. 2005 ;Nowacki et al. 2010 ;Creasy et al. 2017 ;Wolf et al. 2019 ), and differential splitting of SKS and SKKS (e.g.Niu & Perez 2004 ;Deng et al. 2017 ;Grund & Ritter 2018 ;Reiss et al. 2019 ;Tesoniero et al. 2020 ;Asplet et al. 2020 ).An alternative technique makes use of the long ray path of S diff through the lowermost mantle, inferring deep mantle anisotropy by measuring splitting of S diff (e.g.Vinnik et al. 1989Vinnik et al. , 1995 ; ;Cottaar & Romanowicz 2013 ;Wolf & Long 2023b ;Wolf et al. 2023bWolf et al. , 2024b ) ), typically comparing with SKS splitting to account for any upper mantle contributions.
The sensitivity of seismic body wave phases to mantle seismic anisotropy is often investigated on a phase-by-phase basis, sometimes b y appl ying splitting techniques to synthetic seismo grams that were calculated for known Earth input models (e.g.Nowacki & W ookey 2016 ;W olf et al. 2022b ).Other times, the sensitivity of individual seismic phases in an anisotropic Earth is investigated by calculating sensitivity kernels (e.g.Favier et al. 2004 ;Sieminski et al. 2007Sieminski et al. , 2008 ) ).While such approaches are very helpful for a detailed understanding of how and where individual phases are affected by seismic anisotropy, they only focus on a fraction of the seismic wavefield.In this study, we aim for a systematic search of the seismic wavefield in order to uncover nov el splitting strate gies.We implement a wav efield differencing approach that allows us to systematicall y anal yze the ef fects of anisotropy on the entire seismic wavefield, and to investigate which phases are most sensitive to anisotropy in which portions of the mantle.
While commonly used shear-wave splitting strategies using phases such as SKS, SKKS and ScS continue to yield valuable information on anisotropy at various depths in the mantle, they all have limitations, including those imposed by the distribution of seismic stations and earthquakes at the rele v ant distance ranges.Expanding the repertoire of seismic body wave phases that can be used for shear-wave splitting analysis is desirable, as this would allow for splitting data sets with better spatial and azimuthal coverage.In service of this goal, in this study we carry out a systematic investigation of how seismic anisotropy located at different depths in the Earth's mantle expresses itself in the global seismic wavefield.We analyze a large number of body wave phases (Fig. 1 ), many of which are not routinely used for shear-wave splitting measurements.Our results point towards potential new (or rarely used) techniques to map upper mantle anisotropy using shear-wave splitting, including those that rely on PS and PcS phases.We provide proof-of-concept examples of these strategies applied to real data and discuss how they might improve our ability to image anisotropy, deformation, and flow in the Earth's mantle.

Global w av efield simulations
We conduct global wavefield simulations using AxiSEM3D (Leng et al. 2016(Leng et al. , 2019 ) ), w hich is capab le of handling any 3-D input model and arbitrary seismic anisotropy (Tesoniero et al. 2020 ).In this work, we focus on axisymmetric simulations, for which AxiSEM3D is as efficient as its older relativ e AxiSEM (Nissen-Me yer et al. 2014 ).We compute global wave propagation simulations down to ∼5 s period using AxiSEM3D, following the methodology applied in previous work (e.g.Wolf et al. 2022a ).We al wa ys use a smoothed version of isotropic PREM (Dziewonski & Anderson 1981 , see below) as our background model and al wa ys consider (PREM-) attenuation as well as Earth's ellipticity.We intentionally keep our input models simple to generally assess where in the seismic wavefield seismic anisotropy manifests itself.Effects of laterally changing seismic anisotropy and 3-D input models have been addressed in previous work (e.g.Wolf et al. 2022a , b ).
The typical source-receiver configuration for our simulations is shown in Fig. 2 (a).We place our 100 m deep seismic source at the North Pole, simulating either a normal or strike-slip fault earthquake (Fig. 2 a).The very shallow focal depth helps to avoid surface reflections (depth phases) in our seismograms.We select two different focal mechanisms (Fig. 2 a); their details are only important inasmuch as they influence the initial polarization of the seismic wave.The stations are spaced on a re gular 2-de gree latitudelongitude grid, leading to a closer station spacing at the poles than at the equator (Fig. 2 a).This configuration ensures that stations are regularly spaced along the azimuthal direction.We implement seismic anisotropy in our models by replacing the smoothed PREM velocity input model by seismic anisotropy described by a horizontally transversely isotropic (HTI) elastic tensor created using MSAT (Walker & Wookey 2012 ).We consider five different models, each with anisotropy in a different depth range.The depth ranges of anisotropy that we investigate are 24-220 km (layer 1), 220-400 km (layer 2), 400-670 km (layer 3), 670-800 km (layer 4) and 2641-2891 km (layer 5), as shown in Fig. 2 (b).We ensure that the isotropic average of the anisotropy used in each of these layers agrees with our smoothed isotropic PREM model, for which velocities are constant within each layer.The anisotropic strength within our lowermost mantle layer (layer 5) is 2.75 per cent; we adjust the anisotropic strength in each of the other layers such that we would obtain the same splitting delay time ( ∼1 s) for each of them for a vertically incident wave.We implement two different elastic tensor arrangements (Figs 2 c and d), representing cr ystallog raphic preferred orientation of olivine: in one, the direction from which the elastic tensor is sampled changes with azimuth (arrangement 1; F ig. 2 c), w hile in the other it is al wa ys the same (arrangement 2; Fig. 2 d).We implement these two different elastic tensor arrangements to ensure that our results are not strongly affected by the specific direction that the seismic anisotropy is sampled from.Whenever we use the seismograms to calculate absAD and relAD (see below), we cannot output seismograms for all stations from our source-station setup (Fig. 2 a) at an appropriate sampling period (1 s) due to storage limitations.Instead, we output seismograms every 15 • azimuth and every 2 • distance.

Data processing
In order to understand how seismic anisotropy expresses itself, we compare the seismic wavefield computed for our isotropic input model with the wavefields from each of our anisotropic simulations.Specifically, we compute the displacement U i or U , the gradient of displacement D ij or D , the curl of displacement R i or R , stress S ij or S , and strain E ij or E , where i,j = 1,2,3 correspond to the radial (1), transverse (2) and vertical (3) direction.The difference between the isotropic and the anisotropic simulations will be indicated with a δsign in the following, for example δU 1 = U 1, iso − U 1, ani .From this dif ferential w a vefield, w e infer how much different seismic phases are influenced by mantle anisotropy.We can express this either in For elastic tensor arrangement 1, the direction from which the elastic tensor is sampled changes as a function of azimuth.(d) Elastic tensor arrangement 2, presented using the same plotting conventions as in panel (c).For elastic tensor arrangement 2, the direction from which the elastic tensor is sampled is the same for every azimuth.
terms of a scalar quantity that explicitly considers the amplitude of the incoming wave, or as a scalar quantity that is normalized to the amplitude.To do this, we derive the phase-specific absolute and relative normalized integrated apparent difference ( absAD and relAD ) from the delta field.For U 1 and the SS seismic phase absAD can be expressed as: where a denotes the azimuth index, d 1 is the shortest distance at which the phase arrives (but al wa ys > 20 • ), d 2 is the largest distance at which the phase arrives (but al wa ys ≤180 • ), d is the distance index, tSS is the phase arri v al time for the SS phase, and t the time index.Therefore, absAD quantifies how much a particular phase is influenced by mantle anisotropy by analyzing the 20 s after the phase onset, integrating over distance and azimuth and normalizing the v alue b y the distance range over which the phase occurs.We find a time interval length of 20 s especially useful for the more attenuated phases with a relati vel y long dominant period.Similarly, it can be calculated how much a particular phase is affected by mantle anisotropy relative to its isotropic amplitude.This quantity is where U iso denotes the displacement amplitude in the isotropic seismogram.absAD and relAD will also be calculated for For the interpretation of absAD and relAD , it is important to note that we do not explicitly consider phase interference in the calculations, but simply use the 20 s after the phase arri v als whether or not another seismic phase arrives in this time windo w.Moreover , the contribution to absAD is stronger at shorter distances as | δU i | will be larger due to geometric spreading.This is not the case for relAD .Therefore, phases for w hich a bsAD is large likely have large amplitudes for at least a certain distance range and are, at the same time, sensitive to seismic anisotropy.Phases for which relAD is large do not necessarily have high amplitudes but are strongly indicative of seismic anisotropy.

Exploring the full seismic w av efield
To visualize the difference in seismic wave propagation between isotropic and anisotropic simulations, we create movies showing its time evolution.Example snapshots from such movies are shown in Fig. 3 for seismic anisotropy in the upper mantle (upper panel) and the lowermost mantle (lower panel).We also provide Movies S1-S5 that show the time evolution of the differential wavefield for elastic tensor arrangement 1 and seismic anisotropy at different depths (layers 1-5).The ov erall wav efield difference is substantiall y more af fected b y upper mantle anisotropy due to its influence on surface waves.The influence of surface waves becomes less the deeper in the mantle the seismic anisotropy is placed.In the difference snapshots for lowermost mantle anisotropy, different body wave phases that are often used to characterize it can be very well distinguished (e.g.SKS, SKKS, ScS and S diff ).The patterns that are apparent for upper and lowermost mantle anisotropy on a wave front along a particular distance away from the source are a combination of the initial source polarization and the sampling of the elastic tensor from an azimuth-dependent direction (arrangement 1).
While these w avefield dif ference movies are informative, they mainl y visuall y emphasize the phase arri v als that are most affected by mantle seismic anisotropy at each point in time.To also focus on less affected phases, we investigate the differential wavefield at a particular azimuth, plotted on top of a trav eltime curv e (Fig. 4 for upper mantle and transition zone; Fig. 5 for the lower mantle).Fig. 4 demonstrates that, as expected, the deeper the anisotropy is placed in the upper mantle, transition zone, or uppermost lower mantle (la yers 1-4), the w eaker its influence on near surface reverberations and on surface waves.Phases that are strongly influenced by seismic anisotropy in layers 1-4 include, for example, SS, SSS and SSSS; in fact, these phases are affected so strongly that the difference plots show their minor arc siblings as features with an opposite moveout (upper right-hand corner).
These results for layers 1-4 can be compared to those for layer 5 (lowermost mantle anisotropy; Fig. 5 ).The magnitude of the wavefield differences depends on azimuth, but the absolute sensitivity of each seismic phase to lowermost mantle anisotropy is the same for all azimuths.Seismic phases that are primarily influenced by lowermost mantle anisotropy include the phases commonly used to infer its presence (e.g.ScS, S diff , SK(K)S), but also ScSScS and ScSScSScS.Moreover, at certain distances (P)PS seems to be influenced by lowermost mantle anisotropy, likely because it is sampling the lowermost mantle at these distances and merging with (P)PScS (analogous to S/ScS at ∼90-100 • ).The details of the differential w avefield (unsurprisingl y) depend on the focal mechanism and the elastic tensor arrangement (Figs S1 -S4 ).The overall patterns, however, are the same for all of these scenarios.

Influence of seismic anisotropy on individual seismic phases
Next, we calculate a bsAD (F ig. 1 ) for various seismic phases, which shows us how strongly absolute phase amplitudes are affected by mantle anisotropy.These results are shown for δU i in Fig. 6 and for D ij , R i , S ij and E ij in Figs S5 -S12 .As expected from the patterns in Figs 4 and 5 , similar phases are af fected b y seismic anisotropy in lay ers 1-4, w hereas the phases mainly influenced by lowermost mantle anisotropy (layer 5) are different.It is unsurprising that S , SS, SSS and SSSS are strongl y af fected b y upper mantle anisotropy since the y trav el through the upper mantle two or more times.We find that the transverse components ( U 2 ) of SK(K)S phases are particularl y af fected b y mantle anisotropy.This is expected since SK(K)S would be SV-polarized in an isotropic Earth due to the P-SV conversion at the CMB.Indeed, this insight underpins the popularity of SK(K)S splitting as a tool for measuring anisotropy.The fact that PKIKS and PKS seem largely unaffected by mantle anisotropy in this view has to do with the strike-slip focal mechanism used for the simulation (compare to Fig. S6 ); this is therefore a function of the radiation pattern, rather than an actual lack of sensitivity to anisotropy.
Our calculation of absAD for different phases can shed light on whether there are body wave phases that are strongly affected by mantle anisotropy that are not typically used in splitting studies.Such phases, which are strongl y af fected b y seismic anisotropy in layers 1-4 in all simulations (Figs 6 a-d and Figs S5 -S12 ) include PS and PPS.This is especially noteworthy because their polarization is determined by the P -SV conversion through of the surface underside reflection, which makes them much easier to analyze than, for example, SS.We also find that PcS is strongly influenced by upper mantle anisotropy in our simulations.This is partially due to the interference with S at some of the distance range, but in general, PcS splitting is analogous to * KS splitting, although PcS usually has smaller amplitudes (e.g.Liu & Grand 2018 ).For the lowermost mantle case, we find that anisotropy mostl y af fects SKS, SKKS and ScS phases (Fig. 6 e) and no unusual phases that are strongly affected by deep mantle anisotropy, but rarely used to infer its presence, are immediately apparent.
So far, we have analyzed how much seismic phases are affected by mantle anisotropy in terms of absolute amplitude.Results for relative amplitude or relAD are shown in Fig. 7 for δU i .Fig. 7 shows that for each layer, transverse components are particularly affected for phases that travel through the outer core and convert from P to SV at the CMB.This is logical because in the isotropic case transverse component energy is not present in the absence of phase interference.(Note that to avoid di viding b y zero, we add a water level divisor to our calculations).Seismic anisotropy in layers 1-4 tends to affect S, SS, SSS and SSSS increasingly strongly, since the higher multiples sample seismic anisotropy more often.For the lowermost mantle, prominent signals for PKS and ScSScS are visible for relAD which were less apparent for absAD in Fig. 5 (e).
Interestingly, the influence of seismic anisotropy in layers 2-4 on each individual seismic phase is similar (Figs 6 and 7 ).One reason for this is that for these three layers the influence of the anisotropy on surface waves and near surface reverberations is substantially lower than for layer 1.At the same time, layers 2-4 are sampled from a similar incidence angle for most phases, while the angle at which layer 5 is sampled can be very different.By construction, the strength of anisotropy sampled by a vertically incident seismic wave is the same in each layer (Section 2.1 ), and many phases sample layers 2-4 close to vertical for most distances.Therefore, the anisotropic signature of splitting due to anisotropy in each of those three layers that is visible in individual phases is comparable.

D I S T I N G U I S H I N G A N I S O T RO P Y A N D H E T E RO G E N E I T Y A N D T H E I M P O RTA N C E O F S H E A R -WAV E S P L I T T I N G
The results shown in Figs 3-7 demonstrate that anisotropy has a significant effect on the global seismic wa vefield.How ever, it also well known that isotropic widespread heterogeneity significantl y af fects the seismic w av efield.A ke y question, therefore, is how to distinguish the effects of anisotropy and heterogeneity.Here, we use a set of simulations that incorporate isotropic heterogeneity (as expressed in a global mantle tomography model) to illustrate the well-established principle that measurements of shear wave splitting can be used as a clear indicator of seismic anisotropy.
We first re vie w the well-established theoretical basis for shear wave splitting observations (Vinnik et al. 1989a ;Silver & Chan 1991 ).For simplicity, we assume an SV-polarized wave before it samples seismic anisotropy, such as any shear wave which has undergone a P-to-SV conversion (e.g.SKS).We can describe this wave as a harmonic wave with the angular frequency ω.Assuming that ωt < < 1 (with t being time), the radial component R ( t ) can be expressed as

R( t) cos ωt
(3) (Vinnik et al. 1989a ;Silver & Chan 1991 ).After traveling through the anisotropic material, the transverse component T ( t ) will have the shape of the radial component time deri v ati ve R ( t ) where δt is the time lag between the split waves, φ is the fast polarization direction (measured clockwise from the north direction) and α is the backazimuth (which is the same as the initial polarization) for core-refracted phases such as SKS.
A quantity that is also useful to introduce is the splitting intensity (Chevrot 2000 ), SI , defined as The splitting intensity is high if two conditions are met: First, the transverse component has a significant amplitude and , second , the radial time deri v ati ve has the shape of the transverse component.Eq. ( 4) illustrates that both is true in case of shear-wave splitting due to seismic anisotropy.We carry out a set of simulations that illustrates why shear-wave splitting is such a powerful tool for distinguishing the effects of seismic anisotropy and 3-D heterogeneity.We perform three types of simulations using the setup shown in Fig. 8 (a) with a 600 km deep source and otherwise the same parameters as described in Section 2.1 .First, we calculate seismic wave propagation using isotropic PREM as background model.Secondly, we also implement seismic anisotropy in layer 2 (Fig. 2 b); and , third , we replace SKS waveforms (Fig. 8 b) illustrate that while 3-D tomography shifts the SKS arri v al time quite significantly, it does not lead to effects that mimic splitting: No significant transverse component energy is observed for either isotropic scenario, which is why splitting intensities are negligible.For seismic anisotropy in layer 2, ho wever , the estimated splitting intensity is large (above 1).This is because the transverse component takes the shape of the radial component time deri v ati ve, causing an elliptical particle motion.The same is true for other * KS waves such as SKKS (Fig. 8 c) and PKS (Fig. 8 d).For seismic waves that are not SV polarized before sampling seismic anisotropy, the same argument can be made if the seismograms are simply rotated to the direction of initial polarization derived from the particle motion.This is shown for S (Fig. 8 e) and SS (Fig. 8 f) phases.For isotropic simulations, estimated splitting intensities are zero; in contrast, we observe splitting for the input model with seismic anisotropy.
These results show that, cruciall y, hetero geneity does not have the same effect as seismic anisotropy for different seismogram components and different seismic phases.Specificall y, the w aveform effects that seismic anisotropy has on different horizontal components make it distinguishable from isotropic effects through shear-wave splitting measurements.While time shifts caused by 3-D heterogeneity affect the overall seismic wavefield, these effects can be reliably distinguished from the effects of seismic anisotropy, as shown in Fig. 8 (a).Our simulations therefore illustrate the well-established principle of shear-wave splitting can reliably distinguish between anisotropic effects and isotropic heterogeneity (e.g.Vinnik et al. 1989a ;Silver & Chan 1991 ;Silver 1996 ;Romanowicz & Wenk 2017 ;Becker 2020 ).
We also calculate the differential displacement as in Figs 4 and 5 comparing simulations with isotropic PREM with a scenario in which the PREM-mantle is replaced b y 3-D tomo graphy (Fig. S13 ).This exercise is shown for completeness, although, strictly speaking, it is not directly comparable to our anisotropic simulations.
The reason is that the main aspects that we consider when comparing displacements for two isotropic models (one of which includes realistic 3-D structure) are differential phase arrival times and waveforms distortions due to isotropic effects.Such effects were deliberately avoided in our anisotropic setup (for which the differential displacement is primarily influenced by the redistribution of energy for the horizontal components; Fig. 8 ).Thus, the differential displacement for the 3-D isotropic model (Fig. S13 ) is mostly influenced by two factors: the amplitude of the seismic wave (which also plays a role for the anisotropic simulations), and the path length travelled by the wave, which is approximately proportional to the traveltime (which, by construction, plays no role for the anisotropic simulations, except as it relates to the number of times that an anisotropic layer is sampled).We account for the absolute amplitudes by calculating relAD for different seismic phases (Fig. S14 ).Unsurprisingly, the observed relAD values are much different than those for the anisotropic simulations (Fig. 7 ).While for the anisotropic simulations relAD is dominated by the redistribution of energy from one horizontal component to another, for the 3-D mantle model relAD values are much larger for radial and vertical components than for the transverse component.The reason for this is that P and SV waves are coupled while SH, which is recorded on the transverse component, travels independently (in isotropic media).Therefore, transv erse wav eforms are simpler for isotropic simulations and fewer differences are recorded.The results for the 3-D isotropic model comparison thus confirm the validity of our anisotropic wavefield differencing approach.
Ho wever , it is w orth mentioning that there are known cases for which isotropic structure can mimic splitting, specifically for some reflected and diffracted seismic phases.The phases analyzed in Section 3.2 that are most affected by this complication are ScS (and its multiples) as well as S diff .Detailed work on these phases has been conducted in past studies (Komatitsch et al. 2010 ;Figure 7. Same as Fig. 6 for rel AD U i .Borgeaud et al. 2016 ;Nowacki & Wookey 2016 ;Parisi et al. 2018 ;Wolf et al. 2023b ;Wolf & Long 2024 ).

Differential PS-SKS splitting: inferring upper mantle anisotropy close to the PS reflection point
We demonstrated in Section 3 that PS waves are strongly affected by seismic anisotropy in Earth's upper layers.For the purpose of splitting measurements, which are only sensitive to anisotropy structure along the path, PS is initially SV-polarized due to the P-to-S conversion upon reflection at the surface.It is worth pointing out that, because for the measurement of shear-wave splitting an S arri v al at the receiver is needed, the order of S and P is generally not reversible (e.g.we could not measure SP splitting).On the S-leg of the ray path, PS samples seismic anisotropy in the upper mantle close to the reflection point as well as beneath the seismic receiver, and at distances between 90 • and 115 • any potential influence of lowermost mantle anisotropy can be avoided.Beneath the station, SKS samples seismic anisotropy in a very similar way as PS, such that differential PS-SKS splitting would point towards upper mantle anisotropy close to the PS reflection point at the surface (analogous to SKS-SKKS differential splitting for the lowermost mantle).Su & Park ( 1994 ) invoked this argument and measured seismic anisotropy in a location beneath the southwestern P acific Ocean; howev er, this strategy has apparently not been used since.One possible reason for this is that there is no advantage to using the PS-SKS differential splitting technique when the alternative is to infer upper mantle anisotropy close to the potential reflection point by directly measuring SKS splitting at a seismic station located there.Ho wever , seismic stations are not regularly distributed across Earth's surface.Differential PS-SKS splitting is therefore potentially helpful to study anisotropy in regions that are sampled by PS bounce points but are not themselves well instrumented.
Here we present a proof-of concept example of PS-SKS splitting using real data.Before measuring PS and SKS splitting, we bandpass-filter our data, retaining periods between 6 and 25 s.To analyze seismic anisotropy, we use SplitRacer (Reiss & R ümpker 2017 ), a MATLAB-based graphical user interface.SplitRacer calculates the following shear-wave splitting parameters: the fast polarization ( φ), the time delay time ( δt ) and the splitting intensity ( SI ; Chevrot 2000 ) of the split wave (Section 4 ).To measure the first two quantities, SplitRacer uses the transverse energy minimization technique (Silver & Chan 1991 ) along with a cor rected er ror formulation by Walsh et al. ( 2013 ).A strength of SplitRacer is that time windows are picked automaticall y, thereb y ensuring that measurements are independent of a specific choice of time window.
We demonstrate a proof-of-concept example for PS-SKS splitting using station INK (Walpole et al. 2014 ) located in nor theaster n Canada, which exhibits null or nearly null SKS splitting over the entire backazimuthal range.Since INK is a null station (that is, a station that does not exhibit splitting due to upper mantle anisotropy beneath the recei ver), an y splitting contribution to PS waves must be caused by seismic anisotropy on the ray path far away from the station, most likely close to the PS reflection point at the surface.We analyze PS splitting parameters for events that occurred in the southeastern Pacific subduction zones between 10/1995 and 01/2023 at distances between 90 • and 115 • .Our results are shown in Fig. 9 .We obtain ∼100 robust splitting intensity and two robust ( φ, δt ) measurements.We identify a strongly anisotropic region in the upper mantle to the north of Fiji, in which fast polarization directions are oriented approximately south-north (Fig. 9 d).The PS waveforms and splitting diagnostic plots of the event (2009-01-19 03:35:18) are shown in Figs 9 (a)-(c); this event was used to infer the ( φ, δt ) values for this strongly anisotropic region.Fig. 9 (a) shows a clearly split PS phase, while the SKS phase for the same event is null.After correcting for the best-fitting PS splitting parameters (Fig. 9 c), the corrected particle motion is linear (Fig. 9 b), as expected in case of splitting due to seismic anisotropy.
Upper mantle anisotropy has been densely mapped beneath continents; ho wever , for seismic anisotropy beneath ocean basins very little shear-wave splitting data are available (IRIS DMC 2012 ).To the extent that previous splitting measurements have been reported for the oceans, they were often made at ocean island stations (e.g.Fontaine et al. 2007 ), which themselves represent an anomalous tectonic setting, potentially influencing the upper mantle anisotropy beneath.Moreover, sometimes splitting measurements have been obtained through measurements of splitting from direct S waves (Mohiuddin et al. 2015 ;Eakin et al. 2018 ).Alternati vel y, some splitting measurements have been made using ocean bottom seismometers (e.g.Zietlow et al. 2014 ;Lynner 2021 ), which are expensive to install and typically yield relatively noisy data.We therefore recommend that the differential PS-SKS splitting technique be used systematically to map upper mantle anisotropy beneath oceans.
Building on this idea, we recently published an application of the differential PS-SKS splitting strategy suggested here to map seismic anisotropy beneath the Pacific Ocean basin (Wolf & Long 2023a ), using a data set of ∼320 000 seismograms.This new work suggests that the PS-SKS differential splitting technique can also potentially be used for anisotropic tomography approaches (e.g.Mondal & Long 2020 ) to better resolve structure close to the PS bounce point.

PcS beam splitting: inferring upper mantle anisotropy near the r ecei ver
Our results for relAD (Fig. 5 ) suggest that PcS is significantly affected by splitting due to mantle seismic anisotropy.Ho wever , due to its usually low amplitudes, PcS is not commonly used for the purpose of splitting measurements, with a few exceptions (e.g.Murdie & Russo 1999 ;He & Long 2011 ).Here, we apply a recently established beamforming technique to increase PcS signal-to-noise ratios and to measure splitting from the resulting beams.
To beamform PcS phase arri v als, we follow the methodology of Frost et al. ( 2024 ), using data from an event that occurred on 08/24/2011 in the Peru-Brazil border region (Fig. 10 a).Wolf et al. ( 2023a ) demonstrated that shear-wave splitting measurements from beamformed SK(K)S data reflect an average of the single-station splitting from the seismograms used for the beam.Therefore, PcS beam splitting is equi v alent to a laterally averaged splitting contribution.Following Wolf et al. ( 2023a ) and Frost et al. ( 2024 ), w e construct subarra ys of betw een 10 and 20 stations across the contiguous United States, which have a size of approximately 3 • × 3 • .For each subarra y, w e estimate slowness and backazimuth of the incoming wave, and use this information to calculate radial and transverse component beams from velocity seismograms that were bandpass-filtered retaining periods between 4 and 50 s.As suitable given the size of our subarra ys, w e use a curved wave front approach instead of the typical plane-wave approximation (Rost & Thomas 2009 ).To enhance slowness-backazimuth estimates, we measure how similar individual records are to the calculated beam, a quantity known as the F -statistic (Selby 2008 ;Frost et al. 2013 ).The maximum amplitude of the F -trace is used to infer the best-fitting slowness and backazimuth values.To be able to fairly compare waveforms and amplitudes between different components, we then calculate the linearly stacked beam using these slowness and backazimuth values for the unfiltered data.
We measure PcS beam splitting using SplitRacer, analo gousl y to how we use it for PS and SKS splitting.We first apply a bandpass filter to our beamformed data, retaining periods between 4 and 25 s (Fig. 10 d).This approach leads to robust splitting intensity measurements for most of the subarra ys; how ever, w e cannot obtain any well-constrained ( φ, δt ) measurements because splitting is generally weak, below the detection level for the transverse energy minimization method at these periods.Ho wever , the beamformed waveforms show clear PcS signals and relati vel y low noise levels.We therefore next measure splitting at frequencies that are higher than usual for shear-wave splitting measurements, retaining periods between 1-10 s (Fig. 10 e).While the splitting intensity results are very similar to those obtained with the 4-25 s bandpass-filter, we successfully measure robust ( φ, δt ) splitting parameters for nine subarrays.The inclusion of higher frequency energy helps to resolve ( φ, δt ) because the time delay is larger compared to the dominant period of the signal.The measured ( φ, δt ) splitting parameters we obtain with this approach are similar to pre viousl y published * KS splitting measurements compiled in the IRIS shear-wave splitting database (IRIS DMC 2012 ).
Given this successful example, we suggest that PcS beam splitting measurements can be used in the future to fill in gaps in backazimuthal coverage that often exist for * KS measurements.PcS is useful because it is recorded at shorter distances than * KS (up to ∼63 • ) and therefore can largely increase the number of usable seismic ev ents.A re gion for which PcS beam splitting measurements might be particularly helpful to improve backazimuthal coverage is Japan, which also has a densely spaced seismic stations suitable for beamforming, and which suffers from poor backazimuthal coverage for SK(K)S (Long & van der Hilst 2005 ).Further, we suggest that splitting measurements from beamformed data can generally be made at higher frequencies, which can increase the number of robust ( φ, δt ) measurements in cases of weak splitting, as shown in our example (Fig. 10 e).The use of phases such as PcS that are recorded at a shorter distance than * KS also potentially has the effect that less high frequency energy is lost through attenuation along the ray path.

D I S C U S S I O N
The results we obtained in our w avefield dif ferencing approach show which parts of the seismic wavefield are sensitive to mantle seismic anisotropy.Ho wever , even if effects of mantle heterogeneity can be distinguished from seismic anisotropy, which is true for many cases (Fig. 8 ; e.g.Silver & Chan 1991 ;Long & Silver ;No wacki et al. 2010 ;Romano wicz & Wenk 2017 ) but not all (e.g.Komatitsch et al. 2010 ;Borgeaud et al. 2016 ;Parisi et al. 2018 ;Wolf & Long 2024 ), it is not al wa ys straightforward to determine where along the ray path seismic anisotropy is located.For example, as demonstrated in Section 3 , S -wave splitting can be strongly influenced by upper mantle anisotropy.This is why S is often used to characterize anisotropy near the earthquake source (e.g.Russo & Silver 1994 ;Foley & Long 2011 ;Walpole et al. 2017 ;Eakin et al. 2018 ).Ho wever , since there are limits to ho w reliable explicit corrections for receiver side anisotropy are Wolf et al. ( 2022b ), the source-side S-wave splitting technique should be applied with caution.
While S is relati vel y simple because it samples the upper mantle only twice (once on the downgoing leg near the source and once on the upgoing leg near the receiver), other phases can be more complicated.For example, the SS phase potentially samples upper mantle anisotropy four times: at the source side, twice close to the reflection point, and beneath the receiver.This renders the phase practically unusable for inferring the precise geometry of upper mantle anisotropy in any particular location, although this has been tried (e.g.Wolfe & Silver 1998 ).Theoretically, upper mantle anisotropy should be known along three of the four legs of the ray path through the upper mantle to infer seismic anisotropy for the unknown leg, which is unrealistic for almost all possible source-receiver configurations.Alternati vel y, the sampled reflection point anisotropy could be assumed to be the same along the up-and downgoing leg for SS.Even if this assumption was justified, it would likely be practically impossible to reliably infer the anisotropy, at least using explicit ray-theory based corrections (Wolf et al. 2022b ).This issue is applicable to any seismic phase that may be strongly influenced by mantle anisotropy and samples anisotropy in many separate locations.Therefore, while the splitting of phases such as SS, SSS, etc. may provide a general indication of the presence of seismic anisotropy along the ray path through analysis of their particle motions, they are not well suited to studying the precise geometry of anisotropy in any given region.
To infer lowermost mantle seismic anisotropy, we usually face a similar challenge.All phases that potentially sample lowermost mantle anisotropy also potentially sample seismic anisotropy in the upper mantle, which has to somehow be accounted for in studies of deep mantle anisotropy.Many of the phases that w e ha ve shown to be strongly influenced by lowermost mantle anisotropy (Section 3 ) are already commonly used to diagnose it.For example, SKS-SKKS differential splitting studies are common (e.g.Niu & Perez 2004 ;Reiss et al. 2019 ), S and ScS are often analyzed together to infer lowermost mantle anisotropy (Wookey et al. 2005 ;Nowacki et al. 2010 ;Wolf et al. 2019 ) and S diff waves are also commonly used to infer lowermost mantle anisotropy (e.g.Vinnik et al. ;Cottaar & Romanowicz 2013 ;Wolf et al. 2023b ).Our results do not clearly point to unusual phases that can be used to infer lowermost mantle anisotropy in future studies.While ScSScS unsurprisingly is strongl y influenced b y deep mantle anisotropy, it suf fers from the same technical limitations as phases such as SS, SSS, etc., as discussed pre viousl y.F ig. 5 shows that bey ond around 120 • , PS starts to be influenced by lowermost mantle anisotropy.However, it is unclear how this could be exploited in practice.In theory, a PS-PScS differential splitting technique could be applied in an analogous manner to S-ScS differential splitting.Ho wever , a complication is that PS and PScS do not generally have the same reflection point at the surface.Therefore, splitting close to the reflection point would be hard to characterize for PScS, making it impossible to distinguish it from a lowermost mantle contribution.

F U T U R E D I R E C T I O N S
Our wavefield differencing results point towards new strategies to measure and characterize seismic anisotropy using shear-wave splitting observations, two of which we demonstrated as proof-ofconcept examples.In order to identify the region of the mantle from which a splitting signal originates, the wave should (ideally) not sample seismic anisotropy in multiple different locations along the ray path.Additionally, to enable splitting analysis the initial polarization of the wave should be known.These two conditions are satisfied for the two novel splitting strategies we explore here: PS-SKS differential splitting and PcS beam splitting.SKS, PS and PcS all involve initially SV-polarized S waves due to the P -SV conversion either at the CMB or at the surface.
PcS phases do not usually have high amplitudes (Liu & Grand 2018 ), which is why it is helpful to use these waves in a beamforming framework.We suggest that it is often more helpful to measure splitting for seismic phases with lower amplitudes, for which it is straightforward to understand where along the ray path mantle anisotropy is sampled, than seismic phases that may be strongly influenced by mantle anisotropy, but in a complicated manner (for example, SS).Such low amplitude seismic phases can be enhanced using array techniques, making use of the fact that the splitting parameters measured from beams or stacks approximately agree with the average splitting parameters of the single-station seismograms contributing to them (Wolf et al. 2023a ).Such an approach has been applied to S3KS pre viousl y (Wolf et al. 2023a ) and to relati vel y high-frequency PcS in this work.Based on these findings, we suggest to apply beamforming routinely in studies that measure shear-wave splitting.For example, Fig. 5 shows that there are minor arc seismic phases affected by lowermost mantle anisotropy.Such phases, like minor arc SKKS or S3KS, can be used for measurements of mantle anisotropy more commonly than they currently are.Similarly, the measurement of differential PPS-SKS or PPS-SKKS splitting is concei v able, b y enhancing signals via beamforming if necessary.Such new splitting strategies involving minor arc phases and unusual source-receiver configurations will be helpful to infer upper mantle anisotropy in locations that suffer from poor coverage, as well as to study lowermost mantle anisotropy, which is often only detectable for specific ray path configurations.

C O N C L U S I O N
In this work, we have applied a wavefield differencing approach to analyze the differences betw een wa vefield for isotropic and anisotropic models, incorporating seismic anisotropy at different mantle depths.These wavefield differencing results demonstrate which seismic phases are most strongly influenced by mantle anisotropy.We show that some seismic phases are more suitable than others to infer splitting, even if they are influenced by mantle anisotropy to a similar degree.In particular, we suggest that the PS-SKS differential splitting technique can be commonly used to infer upper mantle seismic anisotropy beneath ocean basins, based on a proof-of-concept example using a station in Canada and earthquakes in the western Pacific.Additionally, seismic phases that are not usually used for splitting measurements because of their low amplitudes should be more routinely analyzed in areas of dense seismic array deployments using a beamforming approach.As a proof-of-concept example, we calculate high-frequency PcS beam splitting for one seismic event for stations across the United States.The w avefield dif ferencing results presented here may inform the design of future studies of mantle anisotropy using body waves.

Figure 1 .
Figure 1.Seismic phases used and this study and their tra veltimes.(a) Ra y path sketch for seismic phases commonly used to infer upper (SKS, SKKS and PKS) and lowermost (SKS, SKKS, S , ScS and S diff ) mantle anisotropy.(b) Seismic phases that usually have high amplitudes in seismograms investigated in this study.See legend for color key.(c) T ra veltime curve for all seismic phases presented in panel (b), using the same colors.The traveltime curve was calculated for a 100 m deep event, as used in the global wavefield simulations.

Figure 2 .
Figure 2. Setup for global wavefield simulations.(a) Source (yellow star) -station (black dots on map) configuration.Events either have a normal (left) or a strike-slip (right) geometry (top).(b) Depths at which we incorporate seismic anisotropy into our AxiSEM3D input models.Layer 1: Upper mantle; layer 2: low er upper mantle; la yer 3: transition zone; la yer 4: upper low er mantle; la yer 5: low ermost mantle.(c) Elastic tensor arrangement 1.The source is represented as a yellow star on the map, which uses a north pole centered projection, and the horizontally transversely isotropic elastic tensors are shown as stereoplots.For elastic tensor arrangement 1, the direction from which the elastic tensor is sampled changes as a function of azimuth.(d) Elastic tensor arrangement 2, presented using the same plotting conventions as in panel (c).For elastic tensor arrangement 2, the direction from which the elastic tensor is sampled is the same for every azimuth.

Figure 3 .
Figure 3.Time snapshots of the differential displacement wavefield (see color bar) for vertical, radial and transverse components from 1000 s (left) to 2000 s (right) after event origin time, for elastic tensor arrangement 1.The top panel shows δU i for seismic anisotropy in the lower upper mantle and the bottom panel for the lowermost mantle.Different seismic phases are marked on the snapshots.

Figure 4 .
Figure 4. Differential displacement wavefield (color bar) as a function of distance for an azimuth of 0 • , calculated using a strike-slip source and elastic tensor arrangement 1 (Fig. 2 c).(a) The dif ferential w avefield for seismic anisotropy in layer 1 (upper mantle) is plotted underneath the trav eltime curv e shown in Fig. 1 (c) with the phase traveltimes shown as black lines.The wavefield difference is presented for radial (left), transverse (middle) and vertical (right) components.(b) Same as panel (a) for seismic anisotropy in layer 2 (lower upper mantle).(c) Same as panel (a) for seismic anisotropy in layer 3 (transition zone).(d) Same as panel (a) for seismic anisotropy in layer 4 (upper lower mantle).White vertical stripes are due to the plotting convention are not to be interpreted.the PREM mantle by the 3-D tomography model S40RTS without incorporating seismic anisotropy.The results of these simulations are shown in Fig. 8 .SKS waveforms (Fig.8b) illustrate that while 3-D tomography shifts the SKS arri v al time quite significantly, it does not lead to effects that mimic splitting: No significant transverse component energy is observed for either isotropic scenario, which is why splitting intensities are negligible.For seismic anisotropy in layer 2, ho wever , the estimated splitting intensity is large (above 1).This is because the transverse component takes the shape of the radial component time deri v ati ve, causing an elliptical particle motion.The same is true for other * KS waves such as SKKS (Fig.8 c) and

Figure 6 .
Figure 6.abs AD U i (see text) for radial ( U 1 , blue), transverse ( U 2 , g reen) and ver tical ( U 3 , red) components for most seismic phases shown in Fig. 1 (b).(a) ab s AD U i , nor malized to the largest amplitude of either U 1 , U 2 or U 3 , for seismic anisotropy in layer 1 (upper mantle).b) Same as panel (a) for seismic anisotropy in layer 2 (lower upper mantle).(c) Same as panel (a) for seismic anisotropy in layer 3 (transition zone).(d) Same as panel (a) for seismic anisotropy in layer 4 (upper lower mantle).(e) Same as panel (a) for seismic anisotropy in layer 5 (lowermost mantle).For the calculation of abs AD U i and rel AD U i we sum over all calculated azimuths as laid out in eqs ( 1 ) and ( 2 ).

Figure 8 .
Figure 8. Representation of the differential effects that seismic anisotropy and 3-D velocity heterogeneity have on seismic waveforms.Results are from simulations that use (1) background model isotropic PREM (Dziewonski & Anderson 1981 , black), (2) isotropic PREM and anisotropic layer 2 (blue) and (3) replacing the PREM-mantle by the tomography model S40RTS (Ritsema et al. 2011 ; red).(a) Source (yellow star) -receiver (black circles) configuration for phases S, SS, SKS, SKKS and PKS at randomly chosen distances at which they are recorded.(b) Top: SKS waveforms for all three types of simulations.Significant transverse component energy only arrives for the anisotropic simulations.3-D heterogeneity shifts the arrival time compared to isotropic PREM.Bottom: Corresponding particle motions and SI values (bottom right in particle motion plots).SI values are 0 for isotropic simulations.(c) Particle motions and SI values for the SKKS phase.(d) Same as (c) for PKS.(e) Same as (c) for S. Note that particle motions are rotated in a coordinate system determined by the initial polarization of the wave.(f) Same as (e) for SS.

Figure 9 .
Figure 9. Differential PS-SKS splitting results for null station INK (Walpole et al. 2014 ).(a) PS (top) and SKS (bottom) radial and transverse component waveforms for an event that occurred on 01/19/2009.The PREM-predicted phase arri v al time is shown with a green line and the start/end of the automatically selected measurement windows are shown with orange lines.SKS transverse energy is at the amplitude level of the noise; therefore, SKS splitting is null.(b) PS particle motions before (left) and after (right) correcting for the best fitting splitting parameters.(c) Best-fitting splitting parameters in the φ-δt plane.The 95 per cent confidence interval is shown in black, red crossing lines show the best fitting combination of ( φ, δt ).(d) Map view of the source (orange stars)receiver (red triangle) configuration for the differential PS-SKS splitting analysis.Ray paths are shown as gray lines and colored circles indicate absolute PS splitting intensities (see legend).Black lines represent ( φ, δt ) splitting parameters.Bottom right: Zoom-in for the region of strong upper mantle anisotropy.

Figure 10 .
Figure 10.Summary of PcS beam splitting measurements.(a) The event used for the beamforming, which occurred on 08/24/2011, is represented as a yellow star (moment magnitude: 7.0; depth: 149 km).Subarray central stations across US are shown as black circles.(b) Radial (R) and transverse (T) component PcS beam waveforms for an example subarray with central station 537A.The PREM-predicted phase arri v al time is shown by a green line and the start/end of the automatically selected measurement windows is presented by orange lines.(c) Best fitting splitting parameters ( φ , δt ) for the waveforms shown in panel (b).φ denotes the fast polarization direction measured clockwise from the backazimuthal direction.The 95 per cent confidence interval is represented in black, with contour lines showing different transverse energy component levels.(d) Splitting measurements after applying a 4-25 s bandpass filter.Colors represent the splitting intensity (see legend), and are plotted a the central station location of each subarray.(e) Similar to panel (d) for a bandpass filter between 1 and 10 s.Nine well-constrained ( φ, δt ) measurements (red sticks) are obtained.