Shear-wave attenuation anisotropy: a new constraint on mantle melt near the Main Ethiopian Rift

The behaviour of fluids in preferentially aligned fractures plays an important role in a range of dynamic processes within the Earth. In the near-surface, understanding systems of fluid-filled fractures is crucial for applications such as geothermal energy production, monitoring CO2 storage sites, and exploration for metalliferous sub-volcanic brines. Mantle melting is a key geodynamic process, exerting control over its composition and dynamic processes. Upper mantle melting weakens the lithosphere, facilitating rifting and other surface expressions of tectonic processes. Aligned fluid-filled fractures are an efficient mechanism for seismic velocity anisotropy, requiring very low volume fractions, but such rock physics models also predict significant shear-wave attenuation anisotropy. In comparison, the attenuation anisotropy expected for crys-talpreferredorietationmechanismsisnegligibleorwouldonlyoperateoutsideoftheseismicfrequencyband. Herewedemonstrateanewmethodformeasuringshear-waveattenuationanisotropy,applyittosynthetic examples,andmakethefirstmeasurementsofSKSattenuationanisotropyusingdatarecordedatthestation FURI,inEthiopia.AtFURIwemeasureattenuationanisotropywherethefastshear-wavehasbeenmoreat-tenuatedthantheslowshear-wave.Thiscanbeexplainedbythepresenceofalignedfluids,mostprobably melts,intheuppermantleusingaporoelasticsquirtflowmodel.Modellingofthisresultsuggeststhata1% meltfraction,hostedinalignedfracturesdippingca.40°thatstrikeperpendiculartotheMainEthiopianRift,is requiredtoexplaintheobservedattenuationanisotropy.ThisagreeswithpreviousSKSshear-wavesplitting analysiswhichsuggesteda1%meltfractionbeneathFURI.Theinterpretedfracturestrikeanddip,however, disagreeswithpreviousworkintheregionwhichinterpretssub-verticalmeltinclusionsalignedparallelto theMainEthiopianRiftwhichonlyproduceattenuationanisotropywheretheslowshear-waveismoreatten-uated.Theseresultsshowthatattenuationanisotropycouldbeausefultoolfordetectingmantlemelt


Introduction
The presence of fluids within a fractured host rock has important effects on its seismic and mechanical properties.In the crust, there are many systems where the presence of fluids is critical.These include melt-water processes.Observed low seismic velocity zones in the mantle transition zone (e.g., Schmandt et al., 2014;Liu et al., 2016b) and ultra-low velocity zones (ULVZs, e.g., Liu et al., 2016a;Li et al., 2022) in the lowermost mantle have been interpreted in terms of melt.Aligned melt pockets are a very efficient mechanism for generating seismic anisotropy (e.g., Kendall, 2000;Holtzman and Kendall, 2010).This makes it difficult to discriminate between melt, other shape-preferred orientation (such as dry cracks), and lattice-preferred orientation models of seismic anisotropy from the crust (e.g., Bacon et al., 2022) to the lowermost mantle (e.g., Asplet et al., 2022).
Rock physics models predict that aligned sets of fluidfilled fractures, or melt inclusions, produce an effective medium that exhibits both velocity and attenuation anisotropy (e.g., Hudson, 1980;Chapman, 2003;Jin et al., 2018).This result can be achieved by either treating cracks as scatterers (Figure 1a,b;Hudson, 1980) or through the poroelastic squirt flow of fluids in saturated (or partially saturated) meso-scale fractures (Figure 1d,e;Chapman, 2003;Galvin and Gurevich, 2009;Rubino and Holliger, 2012;Jin et al., 2018;Solazzi et al., 2021).The squirt flow model, in particular, predicts a strong dependence of attenuation anisotropy on the presence of fluids (such as melt) and fracture properties.
Whilst attenuation anisotropy can be observed for P waves (e.g., Liu et al., 2007;Ford et al., 2022) it is the attenuation of S-waves that interests us here.Both the crack scattering (Figure 1b) and squirt flow (Figure 1e) models predict an attenuation anisotropy which can be used to complement studies that measure velocity anisotropy using shear-wave splitting (e.g.Kendall et al., 2005;Verdon and Kendall, 2011;Al-Harrasi et al., 2011;Baird et al., 2013Baird et al., , 2015;;Bacon et al., 2022;Schlaphorst et al., 2022).Attenuation anisotropy is a highly sensitive tool for detecting fluids within the earth that are hosted within aligned fractures.For microseismic settings, where the mechanism of seismic anisotropy is known to be fluid-filled fractures, measurements of anisotropic attenuation in shear-waves have been used to help constrain fracture and fluid properties (Carter and Kendall, 2006;Usher et al., 2017).Attenuation anisotropy can be observed directly in experiments (e.g., Best et al., 2007;Zhubayev et al., 2016), albeit at higher frequencies.Numerical models also show that attenuation anisotropy is sensitive to fluid transport properties (Wenzlau et al., 2010).
Measurements of differential attenuation between different teleseismic shear-wave phases, typically S-ScS, have been previously used to measure isotropic Q S in the Earth's mantle (e.g., Lawrence and Wysession, 2006;Ford et al., 2012;Durand et al., 2013;Liu and Grand, 2018).This differential attenuation can be measured by either taking log-spectral ratios or by measuring instantaneous frequency relative to a reference seismogram (Matheney and Nowack, 1995).Here we employ an instantaneous frequency method, which has been shown to be more robust than spectral ratios for teleseismic shear-waves (Ford et al., 2012;Durand et al., 2013).By making measurements of differential attenuation between fast and slow split shear-waves it is pos-sible to measure attenuation anisotropy.As attenuation anisotropy is primarily predicted by effective medium models of fluid-filled fractures, these measurements are highly sensitive to the presence of fluids, such as melt, within the Earth.
We outline how an instantaneous frequency matching method can be applied to measure attenuation anisotropy using shear-wave splitting.Using synthetic shear-wave data we demonstrate the frequency domain effects of attenuation anisotropy and the implications this can have for measurements of shear-wave splitting.We explore the pitfalls of measuring attenuation anisotropy and demonstrate the efficacy of our instantaneous frequency-matching method.We then demonstrate the application of joint measurements of attenuation anisotropy and shear-wave splitting using SKS data recorded at FURI, Ethiopia.

Models of attenuation anisotropy
When a shear-wave propagates through an anisotropic medium, seismic birefringence -or shear-wave splitting -occurs.The fast and slow shear-waves are polarised along the fast velocity direction and an (assumed) orthogonal direction and propagate at different velocities through the medium.This introduces a time delay between the two and can decouple the two (quasi) shear-waves, although in the teleseismic case the time delay time, δt, is much less than the dominant period of the waveform.Assuming that the medium can be described by a single elastic tensor c ijkl the phase velocities and polarisation of each wave can be found by solving the Christoffel equation, where V is phase velocity, ρ is density, p k is polarisation unit vector and n j,l are propagation unit vectors.Solving this eigenproblem yields three positive, real eigenvalues corresponding to ρV P , ρV S1 , ρV S2 with corresponding eigenvectors describing the polarisation directions, which are mutually perpendicular (Mainprice, 2015).
If the medium is also attenuating, then both shearwaves experience a frequency-dependent loss in amplitude and dispersion.The isotropic attenuation of a shear-wave over its path length, l, can be described by the anelastic delay time t * which is given by where v S is the isotropic shear-wave velocity and 1/Q S is the isotropic shear-wave dissipation coefficient.It can be shown that an attenuating medium requires frequency-dependent velocities, or physical dispersion, where the intrinsic seismic velocity of waves propagating through a medium varies with frequency (Aki and Richards, 1980).If this physical dispersion is also anisotropic, then the seismic velocity anisotropy is frequency-dependent and it follows that attenuation is anisotropic also (Carter and Kendall, 2006).
In the case of an anisotropic attenuating medium, where the shear-wave dissipation coefficient 1/Q S varies with propagation direction, the fast and slow split shear-waves will experience different anelastic delay times.We define this difference in anelastic delay times as where S1 is the fast split shear-wave and S2 is the slow split shear-wave.Following this definition, a positive ∆t * represents the case where the slow shear-wave is more attenuated than the fast shear-wave and a negative ∆t * is where the fast shear-wave is more attenuated that the slow shear-wave.It is also worth noting that due to the definition of anelastic delay time (3) velocity anisotropy will produce a ∆t * even if there is isotropic attenuation (i.e., where This effect, however, due to the difference in travel times through the attenuating medium, is negligible compared to the ∆t * that can be predicted for anisotropic attenuation and will always produce ∆t * > 0 (Supplemental Figure 1).

Anisotropic attenuation due to fluidfilled fractures
We consider two main models of seismic anisotropy due to fluid-filled fractures which also allow for the modelling of attenuation anisotropy.These model attenuation due to scattering (Hudson, 1980) and due to poroelastic squirt flow of the hosted fluids (Chapman, 2003).Hudson (1980) employs an effective medium approach to model attenuation due to preferential scattering by the aligned fractures.For this reason, we refer to this model as crack scattering or simply scattering.The attenuation predicted by this model is anisotropic and frequency-dependent (e.g., Crampin, 1984).Crack scattering also predicts anisotropic attenuation for unsaturated (or dry) aligned cracks, although the attenuation profiles are sufficiently different to allow the dry and saturated cases to be distinguished (Crampin, 1984).The thin layering of material could also produce an effective medium with frequency-dependent anisotropy (Backus, 1962;Werner and Shapiro, 1999) and therefore attenuation anisotropy through a similar scattering mechanism.
There are, however, several limitations to this effective medium approach.It does not model the frequency dependence of the elastic constants, limiting the sensitivity to fracture size, and it neglects the effects of fluid exchange between fractures or between fractures and the host rock matrix.Work to extend the models to include such fluid interchange and equant porosity in the rock matrix show that this has a significant effect on the predicted seismic anisotropy (e.g., Thomsen, 1995;Hudson et al., 1996;Tod, 2001).To adequately model this system an approach that considers the poroelastic squirt flow of fluids held in a random collection of grain-scale microcracks and spherical pores along with aligned meso-scale fracture sets (i.e., fractures much larger than the grain scale) was developed (Chapman, 2003).In the poroelastic squirt flow model, the propagation of a seismic wave causes fluids to migrate between connected meso-scale fracture, micro-scale crack and pore spaces which results in frequency-dependent velocity and attenuation anisotropy.These poroelastic effects can also be modelled by treating the effect of pores and fractures as perturbations in an isotropic background medium (e.g., Jakobsen et al., 2003;Galvin andGurevich, 2009, 2015).More recent developments in squirt flow models allow for partially saturated media (e.g., Rubino and Holliger, 2012;Solazzi et al., 2021) and for multi-phase fluids such as water and supercritical CO 2 (Jin et al., 2018).In both cases, squirt flow predicts attenuation anisotropy but we shall only consider the fully saturated case here.It should be noted that this model, and the scattering model, assume perfectly aligned fractures which is unlikely to represent real-world fracture systems completely.The models are also limited to very low aspect ratios which ultimately derives from the low aspect ratio limit of Eshelby's theory (Eshelby, 1957), which results in very low volume fractions (ca.2×10 −5 ) of fluids required to be in aligned fracture to produce significant velocity and attenuation anisotropy.Recent numerical modelling of squirt flow dispersion models has shown that dispersion increases with fracture density and decreases with aspect ratio, with aspect ratios ≥ 0.1 showing very weak attenuation (Sun et al., 2020).
To calculate seismic velocity and attenuation anisotropy for both the crack scattering and squirt flow models we follow the approach of Crampin (1981).We can include attenuation in the definition of a medium's elastic tensor c ijkl by introducing imaginary parts c I ijkl of complex elastic constants, where the real components c R ijkl are the elastic constants.Solving the Christoffel equation for this complex elastic tensor now yields complex eigenvalues λ = λ R + iλ I , with the dissipation coefficient 1/Q given by the ratio of the imaginary and real parts (Crampin, 1984), For the crack scattering model, the imaginary components of the complex elastic tensor can be constructed using equations from Crampin (1984).Figure 1 shows the seismic velocity (Figure 1a) and attenuation (Figure 1b) profiles modelled as a function of propagation angle relative to the crack normal for a saturated, cracked solid for a frequency of 0.1 Hz.The isotropic solid has velocities v p = 6.5 km s −1 , v s = 3.6 km s −1 and a density of 2700 kg m −3 with fractures which are filled with a fluid with a P-wave velocity v p = 2.7 km s −1 , a crack radius of 5 km, a crack density of 0.05 and an aspect ratio of 1 × 10 −4 .The predicted attenuation anisotropy, ∆t * , as a function of propagation angle (Figure 1c) is calculated using (3) assuming a path length of 50 km through For the fluid-filled models we use an isotropic solid with velocities v P = 6.5 km s −1 and v S = 3.6 km s −1 which contains melt inclusions with v S = 2.7 km s −1 , ρ = 2700 kg m −3 .These parameters are chosen to be broadly consistent with previous effective medium modelling of melt-induced seismic anisotropy (Hammond et al., 2014).
the medium.This broadly represents teleseismic shearwaves propagating through the upper mantle.As ∆t * represents the difference in attenuation between the fast and slow shear-waves there is a discontinuity at θ = 60 • in the scattering model where the polarisation direction of S1 and S2 swap (Figure 1b,c).The importance of this is that crack scattering only predicts ∆t * > 0. It is also worth noting that the scattering model predicts non-physical negative 1/Q values for propagation angles around θ = 45 • due to approximations used to calculate the imaginary components of the elastic tenor.This result can also be seen in Crampin (1984), where the approximations are developed.
The complex elastic tensor for the squirt flow model is calculated following the method of Chapman (2003).As a numerical example, we calculate velocity (Figure 1d), attenuation (Figure 1e), and ∆t * (Figure 1f) as a function of propagation angle for a frequency of 0.1 Hz using the same isotropic solid and crack fill properties and aspect ratio as before.Additionally, we specify a total porosity, Φ = 0.05; a grain-sized microcrack density, c = 0.05; a meso-scale (i.e., larger than grain size) fracture density, f = 0.1; a fracture length, a f = 10 m; and an aspect ratio, r = 1 × 10 −4 .Fracture and microcrack density are related to the respective porosities (or volume fractions) Φ f and Φ c in the squirt flow model by and Φ c = 4 3 π c r (9) (Chapman, 2003).This yields a fracture porosity Φ f = 4.2 × 10 −5 and a microcrack porosity Φ c = 2.1 × 10 −5 , with the remaining porosity modelled as spherical pore spaces.An important assumption of the squirt flow model is that microcracks and pores interact with only one meso-scale fracture, which in turn requires a low fracture density to be valid.We use a mineral-scale relaxation time τ m = 2 × 10 −5 s and grain size ζ = 120 × 10 −6 m, which are taken from Chapman (2003)'s numerical example.From these numerical examples, we can see that the inclusion of poroelastic squirt flow effects has a significant effect on the predicted seismic velocities and attenuation.Furthermore, squirt flow is sensitive to fracture length, with only a small range of fracture lengths producing measurable ∆t * for a given frequency (Figure 2).This frequency range is determined by the characteristic fracture relaxation frequency ω f which is related to fracture length a f by where It follows that different frequencies will induce squirt flow in different fracture sizes (Supplementary Figure 2).In practice, the fractures will not have a uniform length and there will be a range of frequencies.In this modelling the frequency used (0.1 Hz) is assumed to be the dominant frequency of the seismic phases.
The squirt flow model also assumes that the fractures are perfectly aligned.One effect of this assumption of identically sized and perfectly aligned fractures is that squirt flow predicts no attenuation anisotropy when θ = 90°(i.e., propagating parallel to the aligned fractures), whilst crack scattering predicts the maximum ∆t * .The squirt flow model produces a characteristic change in the polarisation of S1 and S2: in the example shown here this occurs at θ = 45 • , but the exact angle where this occurs depends on the model parameters used.Unlike the crack scattering mechanism, this allows for both positive and negative ∆t * .Observing this change of sign in ∆t * (or consistently observing ∆t * < 0) is a clear indicator of squirt flow and, therefore, of the presence of aligned fluid-filled fractures.This has been previously observed in microseismic datasets (Carter and Kendall, 2006;Usher et al., 2017).In particular, squirt flow could reasonably explain the results of Carter and Kendall (2006), who observed some cases where ∆t * < 0 in microseismic data recorded at the Valhall Field, in the Norwegian sector of the North Sea.Fractures on the order of 0.6 m − 6 m would produce attenuation anisotropy for microseismic frequencies (Supplementary Figure 2).Due to the length scales of both squirt flow and crack scattering, we would not expect significant attenuation anisotropy to occur for crystal latticepreferred orientation mechanisms.The effects of velocity anisotropy, where S2 is more attenuated due to its larger travel time, are negligible (Supplementary Figure 1) and even if there are grain-scale fluid inclusions, such as grain boundary wetting, the squirt flow effects would occur well outside of the seismic frequency band.This, combined with the sensitivity of attenuation anisotropy to very low volume fractions of aligned fluid inclusions, makes measuring attenuation anisotropy a promising tool to detect fluids in the subsurface.
3 Instantaneous frequency as a measure of attenuation anisotropy

Instantaneous frequency
As we have shown, the crack scattering and squirt flow mechanisms both predict attenuation anisotropy which we could potentially measure in shear-wave splitting datasets.If the shear-waves S1 and S2 share the same source, geometrical spreading, and effective receiver transfer functions, then they should have equivalent frequency spectra if the intrinsic attenuation along the ray path is isotropic, barring the small difference caused by velocity anisotropy.Therefore, if we can measure a significant difference between the frequency content of each shear-wave, this might be attributed to attenuation anisotropy.
To measure the difference in attenuation between fast and slow shear-waves we apply the instantaneous frequency matching method of Matheney and Nowack (1995).Instantaneous frequency matching has been shown to be less sensitive to noise when measuring attenuation than taking spectral ratios (Matheney and Nowack, 1995;Engelhard, 1996).Instantaneous frequency matching also gives more robust estimates of isotropic mantle attenuation for teleseismic shear-wave phases than spectral ratios (Ford et al., 2012;Durand et al., 2013).This method also does not require the assumption of frequency-independent attenuation, which is useful for the case of fluid-filled fractures where frequency-dependent anisotropic attenuation is predicted even for seismic frequencies (e.g., Chapman, 2003;Jin et al., 2018).A similar approach can be taken by performing frequency shifts in the frequency domain (Quan and Harris, 1997), although we prefer the instantaneous frequency method as all measurements are kept in the time domain.Instantaneous frequency is a concept that arises from complex trace analysis (Gabor, 1946).A time-domain signal x(t), such as a seismic wavelet, can be described in terms of its instantaneous amplitude (or envelope), a(t), and instantaneous phase, which is equivalent to representing the signal by its complex Fourier spectrum (Engelhard, 1996).To construct the complex trace we apply a Hilbert transform to x(t) to give the orthogonal quadrature (or imaginary) trace with the complex trace then given by: From this complex trace, we then obtain the following expressions for instantaneous amplitude, and instantaneous phase The instantaneous frequency of our signal x(t) is given by the rate of change of the instantaneous phase with respect to time (Taner et al., 1979).This requires taking the derivative of an arctangent function, which results in where is a damping factor that can be added to reduce the large positive and negative amplitude spikes that can occur (Matheney and Nowack, 1995).As we did not observe large spikes in our instantaneous frequency traces, and are only interested in a single time window, we did not add a damping factor.The instantaneous frequency values are also weighted by the squared instantaneous amplitude.This gives a damped and weighted instantaneous frequency within a specified analysis window as When weighted by instantaneous amplitude the instantaneous frequency of a signal approaches the centre frequency, or spectral mean, of the signal's Fourier power spectra for a sufficiently large analysis window (Saha, 1987;Barnes, 1993).We use analysis windows picked for shear-wave splitting analysis, which isolate the phase of interest.

Instantaneous frequency matching of split shear-waves
The attenuation of a seismic phase is measured by matching the instantaneous frequency of the observed phase, f obs , to that of a reference phase, f ref .This is done by applying a frequency domain causal attenuation operator, where t * is the anelastic delay time (2) and ω r is the angular reference frequency (Muller, 1984), to the reference phase.Note that D(ω) affects both the amplitude and phase of the waveform, which has important effects when we relate attenuation anisotropy to shearwave splitting.It is also worth noting that this causal attenuation operator is different from the operator stated in Matheney and Nowack (1995), by a factor of −1 for the complex exponent term.Both expressions are valid, with the sign of the phase delay term being chosen to ensure that D(ω) produces a causal signal (Supplementary Figure 3).This depends on the choice of reference frequency and the sign convention of the fast fourier transform (FFT) implementation used.We choose to follow previous work (Ford  et al., 2012;Durand et al., 2013) in using (20).Following Muller (1984) the reference frequency is set to the Nyquist frequency, as this ensures ω < ω r and imposes a negative phase shift for all frequencies when using (20).Another common choice of reference frequency is 1 Hz (e.g., Futterman, 1962;Aki and Richards, 1980;Ford et al., 2012;Durand et al., 2013), which works provided that it is outside the frequency range of interest.One final important point to note, which may be slightly obfuscated by our choice of notation, is that this choice of attenuation operator implicitly assumes that Q is constant with frequency.This is a reasonably safe, and common, assumption to make for the seismic frequency band (e.g., Aki and Richards, 1980).However, this does mean that whilst there is no assumption of constant Q in the measurement of instantaneous frequency (Dasios et al., 2001;Ford et al., 2012), the common choice of D(ω) adds this assumption to the instantaneous frequency matching process.
Where a match in the instantaneous frequencies is achieved (i.e., ∆f = f ref − f obs = 0) the t * operator that is retrieved represents a differential attenuation between f ref and f obs .The physical meaning of the measured differential attenuation depends on the selection of f obs and f ref .For example, to measure lowermost mantle attenuation, the lower mantle transiting S phase can be used as a reference phase for ScS.The differential attenuation between the S and ScS phases can then be attributed to the divergence of the phases' ray paths in the lower mantle (Ford et al., 2012;Durand et al., 2013).
To measure attenuation anisotropy, instead of choosing a separate seismic phase as the reference phase we take advantage of shear-wave splitting and use one of the split shear waves as the reference phase.This gives the differential attenuation between S1 and S2, which we have previously described as ∆t * (equation 3).The sign of ∆t * indicates whether the fast (S1) or slow (S2) shear-wave has experienced more attenuation.
For this method to work, the fast polarisation direction must be correctly identified so that the fast and slow shear-waves can be separated.This is important as shear-wave splitting delay times are typically much smaller than the dominant period of the signal.This assumption is often made for teleseismic shear-waves (e.g., Silver and Chan, 1988;Chevrot, 2000).A consequence of this is that fast and slow shear-waves are not wholly split in time.This causes interference between the two shear-waves if they are viewed in the incorrect reference frame, which consequently affects the apparent frequency content of the two shear-waves.This makes the frequency content of each component dependent on the orientation of the reference frame.This is then further complicated by the phase shift introduced by the causal attenuation operator D(ω).We will expand on this further below, using example synthetic shear-waves.

Frequency domain effects of shear-wave component rotation and attenuation anisotropy
We can use synthetic data to explore the effects of component rotation and attenuation anisotropy on the frequency content of the apparent S1 and S2 phases, without the constraints attached to real observations of shear-wave splitting.All our synthetic examples are generated using a Gabor wavelet ) with a dominant, or carrier, frequency f 0 = 0.2 Hz and a time shift t 0 = 0 s.The parameters γ and ν control the shape of the wavelet.For small γ the wavelet has a delta-like impulse and for large γ it has an oscillatory character.The parameter ν describes the symmetry of the wavelet.For ν = 0, the wavelet is symmetric and when ν = −π 2 or π 2 it is antisymmetric (Červenỳet al., 1977).Here we follow Matheney and Nowack (1995) and use the parameters γ = 4.5 and ν = 2π/5.Synthetics are generated at a sample frequency of 20 Hz.The wavelet is then projected onto horizontal component seismograms from a desired initial source polarisation.Shear-wave splitting is applied to each synthetic by specifying the two desired shear-wave splitting parameters: the fast direction, φ f , and delay time, δt.Where attenuation anisotropy is applied the synthetic is rotated to the fast polarisation direction, to isolate the fast and slow shear-waves, and either the slow trace or the fast trace is attenuated to achieve a positive or negative ∆t * .
Using simple synthetic examples (Figure 3), the effect that attenuation anisotropy has on shear-wave splitting can be seen.Here we generate synthetics with a fast polarisation direction of 30 • , a lag time of 1.5 s and a source polarisation of 70 • (Figure 3a).The slow trace is attenuated by applying a causal attenuation operator (20) where t * = 1.0 s, introducing a differential attenuation (or attenuation anisotropy), ∆t * = 1.0 s (Figure 3b).A negative differential attenuation ∆t * = −1.0s can be introduced by instead attenuating the fast shearwave (Figure 3c).Visual inspection of these synthetics shows a loss of amplitude on the attenuated trace.However, we can also observe an additional time delay in the attenuated traces introduced by the phase terms of the causal attenuation operator (equation 20), which has significant implications for measurements of shearwave splitting.
The effect of component rotation on shear-wave frequency content can be further demonstrated using the synthetics from Figure 3a and 3b.The seismograms are rotated to the geographic reference frame (i.e., where a reference frame rotation φ r = 0 • returns the North and East components) and then rotated through reference frame angles in the range of −90 ≤ φ r ≤ 90.At each φ r the amplitude of the frequency spectra is calculated, along with the instantaneous frequency within a 10 s analysis window centred on the wavelets (Figure 4).In the case where ∆t * = 0 s the spectral amplitude of the fast (Figure 4a) and slow (Figure 4c) shearwaves vary with φ r .When the fast and slow shearwaves are correctly separated at φ r = 30 • or φ r = −60 • there is no difference in the respective instantaneous frequencies, which are both measured as 0.2 Hz.When ∆t * = 1 s is applied the frequency content of the fast shear-wave should be unchanged, which is the case at φ r = 30 • .The additional attenuation applied to the slow shear-wave reduces the effect of component rotation, but the effect is still strong enough to affect our instantaneous frequency matching method.These examples also show that instantaneous frequency retrieves the average amplitude-weighted frequency for each trace.
The synthetic shear-waves shown in Figure 3b, where φ f = 30, δt = 1.5, and ∆t * = 1, can also be used to demonstrate how instantaneous frequency matching can retrieve the applied attenuation anisotropy.Again the synthetic is initially rotated to the geographic reference frame and then rotated over the range −90 ≤ φ r ≤ 90.This simulates searching over the full range of reference frame rotations to test all potential fast shearwave polarisations.At each reference frame rotation the instantaneous frequency of the two horizontal components is measured (Figure 5a), along with the difference in instantaneous frequencies (Figure 5b).The second eigenvalue of the trace covariance matrix, λ 2 , is also calculated after correcting for the lag time δt = 1.5 s (Figure 5c).We calculate λ 2 as it is commonly used in shear-wave splitting analysis that employs eigenvalue minimisation (Silver and Chan, 1991;Wuestefeld et al., 2010;Walsh et al., 2013).Shear-wave splitting introduces a phase delay between the two orthogonally po-larised shear-waves, resulting in an elliptical particle motion.Eigenvalue minimisation methods use λ 2 to characterise this, or the seismic energy on the transverse component of a seismogram viewed in the radialtransverse reference frame.It therefore follows that if attenuation anisotropy introduces an additional phase delay term, this can add a source of systematic error to shear-wave splitting measurements.Only when the data is corrected for the applied ∆t * by attenuating the apparent fast shear-wave, which is the reference phase for a positive ∆t * , is the input fast polarisation direction able to be retrieved (Figure 5c).In this example, we know ∆t * = 1 s and can omit a search over a range of potential ∆t * values.
The instantaneous frequency of the apparent fast and slow shear-waves varies as a function of reference frame rotation φ r (Figure 4, 5a).For both the uncorrected (solid lines) and corrected (dashed) traces there are two points where the instantaneous frequencies match, which can be seen as minima in |∆f | (Figure 5b).In the uncorrected data, these points are separated by approximately 90 • and if ∆t * = 0 then one minima lies at the fast polarisation direction.When there is attenuation anisotropy these minima are not located at the true fast polarisation direction (solid line, Figure 5b).When the correction for ∆t * is applied these minima collapse towards one another, but do not necessarily converge to the same point.
Figure 5c shows the effect that attenuation anisotropy has on shear-wave splitting measurements, as characterised by λ 2 .If we do not correct for the applied attenuation anisotropy then the λ 2 minima can appear to be less pronounced and deflected from the true fast polarisation direction.In this example, this synthetic shearwave splitting has no clear λ 2 minimum when we correct for the imposed delay time of 1 s.The minimum λ 2 occurs at a fast polarisation direction of −78.24 • compared to the true fast polarisation direction of 30 • .When we correct for ∆t * this effect is entirely removed and we can retrieve the input shear-wave splitting parameters.This error in fast polarisation direction increases with ∆t * and may not be fully captured by standard methods of measurement uncertainty estimation such as, for example, using the F-test derived 95% confidence region of the measured λ 2 values (Silver and Chan, 1991;Walsh et al., 2013), as the frequency effects of attenuation anisotropy distort λ 2 with rotation angle (Figure 5c).The magnitude of this effect depends on the strength of attenuation anisotropy.

Grid searching over component rotation and attenuation anisotropy to match instantaneous frequency
These synthetic examples (Figure 5) highlight an important challenge in measuring attenuation anisotropy for shear-waves.The inherent rotational interference between the fast and slow shear-waves makes measuring ∆t * highly dependent on accurately identifying the correct fast polarisation direction.Meanwhile, the error that ∆t * can introduce into shear-wave splitting measurements means that we cannot treat measurements of the fast polarisation direction as independent.To successfully measure ∆t * we must, therefore, also identify the true fast polarisation direction.
One strategy to achieve this is to search over both the potential component (or reference frame) rotation angles φ r and differential attenuation ∆t * .To transform the instantaneous frequency matching process into a minimisation, simplifying the grid search, we adjust the objective function from ∆f , used in Matheney and Nowack (1995) In this form we have to fix the reference and observed traces to allow for automation of the grid search and set the apparent At each reference frame rotation, we then correct for the differential attenuation by attenuating the apparent fast shear-wave (blue) by t * = 1 s and repeat the measurements (dashed lines).The solid vertical line shows the applied (or true) fast polarisation direction, 30 • and the dashed vertical line shows the fast polarisation direction that would be recovered if the synthetics are not corrected for ∆t * .S1 phase as f ref and the apparent S2 phase as f obs , assuming that φ r is the fast polarisation direction.This assumption means that we are unable to immediately determine the sign of ∆t * as cases where ∆t * < 0 are reported at the 90 • from the true fast polarisation direction (i.e., the traces have been rotated such that S2 has become the reference phase).To find the correct sign for a ∆t * measurement we must return to input data, correct for ∆t * and then measure shear-wave splitting.If the measured fast polarisation agrees with φ r , within measurement uncertainty, this indicates a positive ∆t * .
If the difference between the fast polarisation and φ r is approximately 90 • , within measurement uncertainty, this indicates that φ r is the polarisation direction of the slow shear-wave which requires a negative ∆t * .
If we look at grid search results for individual shearwaves (Figure 6a), it becomes clear that we cannot uniquely constrain φ r and ∆t * for a single event using our grid search method.One property of the relationship between the instantaneous frequency of splitshear waves and component rotation that we can take advantage of to resolve this is that instantaneous frequency (as a function of component rotation) is also dependent on the source polarisation of the shear-waves.Performing a grid search over φ r and ∆t * for synthetics with example source polarisations of 45 • (Figure 6a), 130 • (Figure 6b) and 285 • (Figure 6c), we can see that whilst we are unable to retrieve the input parameters φ r = 30 • , ∆t * = 1 s in each case there is a different subset of the model space which minimises |∆f |.For each source polarisation, this subset includes the true model parameters.When the examples are summed, the model space which can minimise |∆f | is greatly reduced (Figure 6).In this simple, low noise example the minima of the sum returns the input φ r , ∆t * exactly.
Therefore, we can measure φ r and ∆t * if we have sufficient measurements of shear-waves with different source polarisations, where the assumption that all shear-waves sample the same attenuation anisotropy can be made.For this stacking method to work well, data with a good spread of source polarisations is desirable.For real data this does place constraints on where measurements can be made, as measuring shearwave splitting from sources with an even distribution of source polarisations that sample a single region of attenuation anisotropy could be challenging.

Synthetic examples
We demonstrate our |∆f | stacking method using synthetic shear-wave data.These examples show that our method can retrieve input shear-wave splitting and attenuation anisotropy parameters.As before, we use a Gabor wavelet and generate a set of 100 synthetic shear-waves.These synthetics are generated with a random source polarisation drawn from a continuous uniform distribution between 0 • and 360 • and with a dominant frequency drawn from f ∼ N (0.1, 0.02).Shearwave splitting, with a fast direction φ f = 30 • and delay time δt = 1.0 s, is applied to all synthetics.Attenuation anisotropy, with ∆t * = 1 s, is applied by attenuating the slow shear-wave.Random white noise with a noise fraction, or noise-to-signal ratio, of 0.075 is also added to the synthetics after rotating the components to the geographic reference frame.This represents a good signal-to-noise ratio, of ca. 13 : 1 for real data as this example is intended to represent the ideal case for attenuation anisotropy measurements.To mimic the preprocessing of real data the synthetics are bandpass filtered, using a two-pole two-pass Butterworth filter with corners of 0.01 Hz and 0.  (Restivo and Helffrich, 1999).Each |∆f | grid is weighted by 1/N , where N is the number of waveforms recorded in a 10 • source polarisation bin.The best-fitting φ r and ∆t * is found by taking the minima of the weighted stack (Figure 7).
To estimate the uncertainties in our measurements, we bootstrap our |∆f | stacking.The 100 |∆f | grids are bootstrap sampled, with replacement, 10,000 times.We repeat the source polarisation weighted stacking for each set of bootstrap samples.The resulting distribution of the minimum |∆f | for each bootstrap sample (Figure 8) can be used to define a 95% confidence region in the stacked |∆f |.An upper-tailed test, where any |∆f | that is below the 95% confidence threshold estimated from the bootstrapping (Figure 8) is con-sidered to reasonably explain our data, is used.This 95% confidence threshold can then be mapped back onto weighted |∆f | stack and estimate the uncertainties of φ r , ∆t * from the length and width of the confidence region (Figure 7), following a similar approach to shear-wave splitting studies (e.g., Wuestefeld et al., 2010;Walsh et al., 2013;Hudson et al., 2023).If the minimum of the weighted |∆f | stack sits outside of this confidence threshold, then this tells us that there is either data polluting the stacks that require removal, or that we are unable to confidently measure ∆t * for that station.
For the synthetic examples attenuation anisotropy parameters φ r = 31 ± 1 • , ∆t * = 0.95 ± 0.16 • are measured for the synthetics where ∆t * = 1 s was imposed (Figure 7a).In the case where ∆t * = −1 s was added, we instead measure φ r = −56 ± 1 • and ∆t * = 1.00 ± 0.08 s (Figure 7b).These results show that the source polarisation stacking method can correctly, and accurately, measure the attenuation anisotropy parameters φ r , ∆t * .It is worth noting that we are not able to exactly retrieve the input parameters as we are only correcting for the difference in frequency content between the fast and slow shear-waves and are not removing the effect of attenuation, which results in a permanent loss of amplitudes.The negative ∆t * example (Figure 7b) shows the expected result from imposing ∆t * > 0 in the grid search.The change in sign is instead mapped into the reference frame rotation, with the minimum |∆f | being approximately 90 • rotated from the fast polarisation direction.This has the effect of setting the slow shearwave as the assumed reference (less attenuated) phase and the fast-shear wave as the observed (more attenuated) phase.This allows for the measurement of both positive and negative ∆t * , which is important to enable us to distinguish between potential mechanisms of attenuation anisotropy.
In these synthetic examples, the sign of ∆t * is known.For real data and experiments, we do not necessar-ily have this a priori information.Determining the sign of ∆t * is very important to measuring attenuation anisotropy as it allows us to distinguish between crack scattering and squirt flow mechanisms (Figure 1c,f).Observing negative ∆t * is potentially a powerful diagnostic for the presence of subsurface fluids, as it cannot be explained by velocity anisotropy and requires attenuation anisotropy due to a more complex mechanism such as squirt flow.In turn, squirt flow requires very small volume fractions of fluids hosted by aligned fractures to generate a measurable ∆t * .To correctly find the sign of ∆t * the most convenient approach is to measure attenuation anisotropy (φ r and ∆t * ) and then use these results to remove the effect of attenuation anisotropy before measuring shear-wave splitting.The measured shear-wave splitting parameters, after correcting for attenuation anisotropy, will tell us the correct fast polarisation direction.If the measured fast polarisation agrees with φ r , within measurement uncertainty, this indicates a positive ∆t * .If the difference between the fast polarisation and φ r is approximately 90 • , within measurement uncertainty, this indicates that φ r is the polarisation direction of the slow shear-wave which requires a negative ∆t * .
This can be demonstrated by measuring shear-wave splitting for two synthetic datasets, where the positive ∆t * synthetics are generated with φ f = 30 • , δt = 1.0 s, ∆t * = 1 s and the negative ∆t * synthetics are generated using φ f = 30 • , δt = 1.0 s, ∆t * = −1 s.Here shear-wave splitting is measured before (Figure 9a,c) and after (Figure 9b,d) correcting for the previously measured attenuation anisotropy (Figure 7).Shearwave splitting is measured using eigenvalue minimisation as implemented in the analysis code SHEBA (Wuestefeld et al., 2010).The individual shear-wave splitting results are then stacked, with each result weighted by the signal-to-noise ratio and the number of measurements within a 10 • back azimuth bin (Restivo and Helffrich, 1999).
The results of our shear-wave splitting measurements highlight two key factors.Firstly, the subtle effects that attenuation anisotropy has on apparent shear-wave splitting are clear.In the case with a positive ∆t * , where the slow shear-wave is more attenuated, the additional phase shift caused by the attenuation anisotropy nearly doubles the delay time relative to the true value (Figure 9a,b).The opposite occurs for a negative ∆t * .When the fast shear-wave is more attenuated it is delayed by the phase term of the attenuation operator, which reduces the delay time.In this example, the effect is sufficiently strong to delay the 'fast' shear-wave such that it arrives after the 'slow' shear-wave, which causes the 90 • rotation in the apparent fast polarisation direction (Figure 9c).In both cases, after correcting for the measured attenuation anisotropy (Figure 7) we can retrieve the input shear-wave splitting parameters with significantly higher accuracy than if no correction had been applied (Figure 9b,d).
In the previous example the input splitting parameters are constant across the synthetics, but it is common to get some scatter in individual shear-wave splitting measurements.To test the effect of this on our source polarisation stacking method we repeat the previous experiment and instead randomly draw 100 samples for φ f and δt from φ f ∼ N (30, 5) and δt ∼ N (1.5, 0.15).This set of shear-wave splitting parameters are then used to generate two sets of synthetics as previously described, where we apply ∆t * = 1 s to one set and ∆t * = −1 s to the other.In these cases we are unable to perfectly retrieve the input attenuation anisotropy, measuring φ r = 30 ± 1 • , ∆t * = 0.95 ± 0.30 • and φ r = −62 ± 1 • , ∆t * = 0.90 ± 0.08 • (Supplementary Figure 4).However when we correct the synthetics using the bestfitting ∆t * and φ r we still significantly improve the accuracy of the shear-wave splitting measurements (Supplementary Figure 5).Again this highlights the effect that attenuation anisotropy can have on shear-wave splitting measurements and that after correcting for attenuation anisotropy, even though the corrections are not perfect, we are broadly able to retrieve the true shearwave splitting parameters (Supplementary Figure 5b,d) even with some scatter in the individual observations.This does, however, increase the uncertainty in the retrieved shear-wave splitting and it should be noted that in both cases it is not possible to exactly retrieve the input shear-wave splitting parameters.

Measuring shear-wave splitting and attenuation anisotropy for FURI, Ethiopia
To demonstrate the potential of ∆t * to detect melt or fluids in the subsurface we choose the station FURI, which is situated on the margin of the Main Ethiopian Rift (MER) close to Addis Ababa.FURI is operated as part of the Global Seismograph Network (Albuquerque Seismological Laboratory/USGS, 2014).We choose this locality as previous SKS shear-wave splitting studies have interpreted seismic anisotropy due to aligned melts beneath the MER (e.g., Ayele et al., 2004;Kendall et al., 2005;Bastow et al., 2010;Hammond et al., 2014).Melt has also been inferred by seismic tomography, using body waves (e.g., Bastow et al., 2008), Rayleigh waves (e.g., Chambers et al., 2022) and ambient noise (e.g., Chambers et al., 2019;Eshetu et al., 2021), receiver functions (e.g., Rychert et al., 2012) and magnetotelluric (Whaler and Hautot, 2006) studies.Aligned melt mechanisms should also produce a strong signal of attenuation anisotropy (Figure 1,2), making the MER and surrounding region a natural target to search for attenuation anisotropy.FURI is one of the few permanent stations in the region, with over 20 years of waveform data available, making it a good station to measure SKS shear-wave splitting.
SKS is an ideal phase to attempt attenuation anisotropy measurements using our source polarisation stacking method.As SKS travels through the outer core as a Pwave, it is only sensitive to anisotropy after it exits the core.Arriving near vertically beneath the station, we can assume that all SKS phases sample the same region of the upper mantle regardless of backazimuth.Furthermore, SKS is radially polarised when exiting the core due to the P-to-S conversion (Hall et al., 2004).Therefore the backazimuthal coverage at FURI (Figure 10) approximately maps onto the achievable source polarisation coverage.Whilst this source polarisation coverage is not ideal (Supplementary Figure 6), it is sufficent to ensure that the |∆f | stacking is stable.For other shear-wave phases, such as teleseismic or local S, it may not be possible to achieve sufficent source-polarisation  (c,d) show the result after we correct the synthetic data using measurements of φ r , ∆t * made using our source polarisation stacking method.The stacked λ 2 surfaces are normalised by the 95% confidence value, indicated by the bold contours, which is derived from an F-test (Silver and Chan, 1991;Restivo and Helffrich, 1999).
coverage without relying on raypaths that sample different regions.In that case, our source polarisation stacking method cannot be applied unless it is clear the S phases are sampling the same region of attenuation anisotropy.Removing the requirement for source polarisation stacking is, therefore, desirable and is an avenue for future research.
It is worth noting that two layers of anisotropy have been suggested across the Main Ethiopian Rift, with the upper layer interpreted as aligned melt pockets and the lower layer associated with density-driven mantle flow due to the African superplume (Hammond et al., 2014).As only the upper layer is likely to host aligned melt inclusions, we do not expect the two-layer problem to have a significant effect on our results.However, it is worth noting that the contribution from the lower layer will introduce additional frequency mixing.We would not expect this to mask the strong attenuation anisotropy predicted for aligned melt inclusions, but this complication may increase the uncertainty in ∆t * .
Data is collected for 584 earthquakes, which are at a sufficient epicentral distance (95 • to 120 • ) for SKS to be visible, recorded at FURI (Figure 10).Only earthquakes with a moment magnitude in the range of 5.5 ≤ M w ≤ 7.0 and a minimum depth of 50 km are used.All earthquake data were requested from the International Seismological Centre (2023) bulletin, with the dataset covering 21 years, from 1st January 2001 to 1st January 2022.Before analysis, all waveforms are corrected for instrument response and we detrend and demean the data using tools available in ObsPy (Beyreuther et al., 2010).
Shear-wave splitting is measured for all 584 SKS waveforms before measuring attenuation anisotropy.Whilst it is not essential to measure shear-wave splitting before attenuation anisotropy, and indeed we have already shown that attenuation anisotropy can affect shearwave splitting measurements (Figure 4, 9), it can be a useful first step in analysis and enables us to manually inspect the waveforms data quality before measuring attenuation anisotropy.The waveform data is filtered using a two-pass two-pole Butterworth filter, with corner frequencies of 0.01 Hz and 0.3 Hz.This enables a direct comparison of our results with previous SKS shearwave splitting station averages (Ayele et al., 2004).The filtered waveforms are visually inspected and analysis window start/end search ranges are picked for waveforms where a clear SKS phase can be picked.This manual inspection reduces the dataset to 73 waveforms where SKS can be clearly identified.Figure 11a shows an example SKS phase used.We then measure shearwave splitting using the shear-wave splitting analysis code SHEBA (Wuestefeld et al., 2010), which utilises the method of Silver and Chan (1991) as updated by Walsh et al. (2013).The optimum shear-wave splitting analysis window, which will also be utilised to measure attenuation anisotropy, is found using cluster analysis (Teanby et al., 2004).At this stage in shear-wave splitting analysis, one might seek to further reduce the dataset, by applying data quality thresholds based on Wuestefeld et al. ( 2010)'s shear-wave splitting quality parameter Q (which is not related to the attenuation quality factor) or by removing results which have large measurement errors in φ f or δt (e.g., Kendall et al., 2005).In this case, we do not want to reduce the size of our dataset as this may remove data that exhibits attenuation anisotropy.
For each SKS phase a |∆f | surface is measured by grid searching over −90 • ≤ φ r ≤ 90 • and 0 s ≤ ∆t * ≤ 4 s, in intervals of 1 • and 0.05 respectively.These measurements use the analysis windows previously defined, using Teanby et al. (2004)'s cluster analysis method, for the corresponding shear-wave splitting measurement.As outlined previously we stack our |∆f | measurements, weighted by source polarisation.Measurement uncertainties are determined by bootstrapping the stacking process as described for the synthetic examples (Supplementary Figure 7).The source polarisation of each waveform is estimated in the shear-wave splitting measurement process by SHEBA (Wuestefeld et al., 2010).From the stacked |∆f | the measured attenuation anisotropy is |∆t * | = 0.45 ± 0.20 s and φ r = −45 ± 3 s.As in the synthetic examples, we are not immediately able to determine the sign of ∆t * .To find the correct sign, each SKS phase must be corrected for the measured attenuation anisotropy.Then the shear-wave splitting of the corrected waveforms, with the effects of attenuation anisotropy removed, can be measured.The attenuation anisotropy corrections are applied by rotating the waveforms to φ r and attenuating the retrieved reference phase by the measured |∆t * |.An example corrected SKS phase is shown in Figure 11b.This has the effect of removing the phase shift introduced by at- tenuation anisotropy, although as this is a station averaged measurement the correction may not be perfect.There will be a permanent loss of amplitudes, but the difference in frequency content between the fast and slow shear-waves should be removed (Figure 11b) and this will not affect measurements of shear-wave splitting.
After correcting the SKS waveforms, we measure station averaged shear-wave splitting of φ f = 40 ± 5 • and δt = 1.60 ± 0.34 s.This result is also consistent with previous work, within measurement uncertainties, but the best-fitting delay time has increased by 0.45 s.As the difference between φ f and φ r is 90 • , within measurement uncertainty, we interpret that the fast shear-wave has been more attenuated than the slow shear-wave and that ∆t * < 0. This gives a final joint measurement of station averaged shear-wave splitting and attenuation anisotropy at FURI of φ f = 40 ± 5 • , δt = 1.60 ± 0.34 s and ∆t * = −0.45± 0.20 s.

Characterising fluid inclusions using velocity and attenuation anisotropy
In our examples, using both synthetic and real data, we have established that we can measure attenuation anisotropy in split shear-waves.Our observation of attenuation anisotropy in real SKS data for FURI, Ethiopia is an important result and corroborates previous work which has interpreted seismic anisotropy in terms of preferentially oriented melt inclusions both beneath FURI (Ayele et al., 2004) and potentially more broadly across the Main Ethiopian Rift (Kendall et al., 2005;Bastow et al., 2010).The additional measurement of attenuation anisotropy gives us further insight into this mech-  anism.The observation of ∆t * = −0.45s can only be explained by the poroelastic squirt flow of a fluid-filled medium given that alternate mechanisms, such as crack scattering or velocity anisotropy effects, always predict that the slow shear-wave should be more strongly attenuated and ∆t * > 0 (Figure 1e,f).Furthermore, our modelling of attenuation anisotropy due to crystal pre- ferred orientation of Olivine shows that, for reasonable mantle Q, the expected ∆t * is at least one order of magnitude smaller than what we observe (Figure 1c,f, Supplementary Figure 1).We would also expect other potential mechanisms for intrinsic attenuation anisotropy for crystal preferred orientation, such as grain boundary melt squirt, to operate at frequencies significantly above the seismic frequency band.This allows us to discount crystal preferred aligment mechanisms, making attenuation anisotropy a good indicator for the presence of aligned fluid inclusions.
With the observed ∆t * = −0.45± 0.20 s strongly suggesting the presence of aligned fluid inclusions, the natural next question is how can we characterise these inclusions.As we have already described, the squirt flow model requires a large set of parameters to characterise a fluid-filled fractured medium.One of the most im-portant parameters to have reasonable constraints on is mineral relaxation time, τ m , which is empirically derived and is proportional to the viscosity of the saturating fluid and inversely proportional to the permeability of the host rock (Chapman et al., 2003).Previous work inverting shear-wave splitting for fracture models using the squirt flow model has shown that the inversion is highly sensitive to the τ m used (Al-Harrasi et al., 2011).It has also been shown that varying τ m has a substantial effect on the expected frequency-dependent seismic velocity anisotropy (Baird et al., 2013).Greater constraints on plausible values for τ m in the upper mantle are required to enable detailed modelling of fracture characteristics.Any modelling of fracture properties is also dependent on the choice of grain size, ζ, and fracture length, a f .Together τ m , ζ and a f describe the fracture scale squirt-flow relaxation time, which is also expressed as the squirt-flow frequency ω f = 2π τ f and determines the frequency range of the fracture-dependent squirt flow effects.Whilst this trade-off makes it difficult to constrain the fracture or grain size, if some reasonable assumptions are made it is still possible to constrain potential fracture orientations.
To search for potential fracture orientations, given the lack of constraint on τ m we make some assumptions to simplify the problem.Outside of τ m , fracture length, and grain size, there are 9 other potential free parameters required to calculate a complex elastic tensor using Chapman (2003)'s squirt flow model.We fix these parameters to the values in Table 1, which leaves fracture strike, dip, and medium thickness as free parameters to search over.Seismic velocities and densities are chosen to be consistent with previous effective medium modelling of the region (Hammond et al., 2014).A total porosity, or melt fraction, of 1%, is chosen, along with a fracture density of 0.1, as previous work suggests SKS shear-wave splitting at FURI could be explained by a melt fraction ≤ 1% (Ayele et al., 2004).This represents a parsimonious choice of model parameters as we seek to explain our observations with a small melt fraction, where the implied fracture porosity (i.e., melt volume fraction hosted in the fractures) φ f = 4.2 × 10 −5 .If SKS is assumed to be vertically incident, then the fracture dip corresponds to the angle to fracture normal used in the earlier numerical examples (Figure 1), and the fracture strike is predominately controlled by the measured fast polarisation direction.This assumption also makes ray path length interchangeable with medium thickness.
We search for the best-fitting medium thickness, l, in the range 50 km ≤ l ≤ 150 km and fracture dip angle, θ, in the range 0 • ≤ θ ≤ 90 • by rotating the elastic tensor to θ and calculating the predicted delay time, δt and attenuation anisotropy, ∆t * .The misfit for these predicted parameters is calculated using a normalised least-squares approach.To reflect the lack of constraint on τ m , and therefore also τ f , in upper mantle conditions this exercise is repeated over a large range of τ m values, 10 × 10 −6 s ≤ τ m ≤ 10 × 10 −2 s, an assumed grain size of 1 mm and fracture lengths of 10 m, 100 m and 1000 m. Figure 14 shows the τ f required by the current choice of grain size, fracture length and τ m (Figure 14a) along with the misfit of the best-fitting model (Figure 14b), predicted ∆t * and δt (Figure 14c,d), and the best-fitting medium thickness (Figure 14e) and fracture dip (Figure 14f) for a given τ m .This modelling exercise can be repeated by fixing an assumed fracture length and varying the chosen grain size, which gives similar results (Supplementary Figure 8).Despite the lack of constraint on τ m one set of model parameters emerges: a medium thickness in the range ca.90 km − 120 km and fracture dip in the range ca.38 • − 48 • which can reasonably explain the observed delay times and δt * .It is worth noting that different modelled fracture lengths and grain sizes require a different range of τ m values to fit the results.Therefore with better constraints on τ m it would be possible to identify plausible fracture (or melt inclusion) lengths for a given grain size.This modelling also shows the value of measuring attenuation anisotropy.In addition to identifying the presence of aligned fluid-filled fractures, measurements of ∆t * add important constraints to fracture orientation.The uncertainty in the measurement of δt = 1.60±0.34s means that it can be reasonably explained by all τ m (Figure 14d) and the additional measurement of ∆t * adds an extra data point.This uncertainty largely maps into melt volume fraction, which has a strong effect on the seismic velocity anisotropy, which we have elected to fix at 1 %, and fracture density, which is required to be low and fixed to 0.1 The measured delay time can also be fitted by shallowly dipping or near-vertical fractures, with the addition of attenuation anisotropy, ∆t * = −0.45±0.20 s, requiring shallowly dipping fractures (Figure 1, Supplemental Figure 9).This relies on the assumption that the squirt flow model (Chapman, 2003) is valid in upper mantle conditions, where this model has not previously been tested.Poroelastic squirt flow requires that the melt is hosted in very low aspect ratio inclusions, due to the limitations of Eshelby's theory (Eshelby, 1957), and that the melt inclusions are near perflectly aligned.To our knowledge poroelastic squirt flow is the only model can explain the negative ∆t * observed.Furthermore, grain-scale melt squirt in the mantle has long been used to model isotropic velocity and attenuation (e.g., Mavko and Nur, 1975;Hammond and Humphreys, 2000).Adding melt squirt of aligned melt inclusions allows us to consider the contributions to velocity and attenuation anisotropy, but this requires that we can model low aspect ratio melt inclusions as fluid-filled fractures.Future work is needed to establish the theoretical attenuation anisotropy of larger aspect ratio melt inclusions, such as melt tubules.
The best-fitting fracture strike direction is found by setting the medium thickness to the thinnest plausible value from the previous modelling exercise, 90 km, and searching over fracture dip and strike angles, where we seek to fit ∆t * , δt and φ f again using a normalised leastsquare cost function (Figure 15).This layer thickness is broadly consistent with previous estimates of the thickness of anisotropy beneath FURI (Ayele et al., 2004), al-Figure 14 Results of modelling the fracture dip and medium thickness, or ray path length, which best explain the observed δt and ∆t * at FURI, Ethiopia using a squirt flow model (Chapman, 2003) though a thinner region of melt inclusions could be accommodated by increasing the melt fraction.We assume a fracture length of 100 m and a grain size of 1 mm and set the mineral scale relaxation time, τ m , to 9.55 × 10 −5 s.The best-fitting orientations give a fracture with a dip of 39 • and an NW-SE strike (Figure 15).This rift perpendicular fracture orientation complicates previous interpretations that seismic anisotropy across the MER is due to rift parallel, vertical melt inclusions in the uppermost mantle (e.g., Ayele et al., 2004;Kendall et al., 2005).It is worth noting that it is only the addition of ∆t * which requires shallowly dipping fractures.This shallow dipping fracture model can then only fit the observed fast polarisation direction, φ f = 40 • if the fractures have a NW or SE strike.Alternatively, melt could be accomodated in inclusions with aspect ratios above the limit set by Eshelby (1957).This scenario is not modelled here, and would require revisiting long standing assumptions of aligned melt mechanisms for seismic anisotropy in the region as low aspect ratio melt inclusions have been the prevailing interpetation (e.g., Ayele et al., 2004;Kendall et al., 2005;Bastow et al., 2010;Hammond et al., 2014).A final possiblity is that there is some other, as yet unknown, mechanism for attenuation anisotropy in the uppermost mantle which can explain our observations whilst allowing for near-vertical  (Chapman, 2003): results for fracture strike and dip that can explain the observed attenuation anisotropy and shear-wave splitting in vertically incident SKS phases at FURI.We model a medium with 1% total melt volume fraction, of which 0.1% is hosted in aligned fractures.Medium thickness is fixed to 90 km following our previous models (Figure 14e), which is consistent with previous estimates of the maximum thickness of anisotropy in the region (Ayele et al., 2004;Hammond et al., 2014).We use a normalised least squares cost function for ∆t, δt, and φ f .Our results favour shallowly dipping, 39 • , fractures which are oriented NW-SE, which is approximately perpendicular to the Main Ethiopian Rift.For details of other model parameters used see text.

melt inclusions.
With only one data point we cannot say if this negative attenuation anisotropy, and the requirement for shallowly dipping fractures, is localised to FURI or is more widespread across the MER.Comparison to a recent shear-wave velocity tomography model at a depth of 100 km does indeed show a low-velocity anomaly which extends perpendicular to the rift axis directly beneath FURI (Figure 10; Chambers et al., 2022).A linearly interpolated cross-section through the model (Figure 16) shows that the feature is up to 70 km thick and situated directly beneath FURI.This feature could represent a network of shallowly dipping aligned melt inclusions extending away from the MER beneath FURI, which is causing the observed attenuation anisotropy.This is slightly thinner than what our models find, but this could potentially be accommodated by a modest increase in the overall melt fraction or fracture density.The melt fraction and fracture density were fixed to 1 % and 0.1 to simplify the modelling done here, but could plausibly be increased.Further work, such as siting several additional stations further along the anomaly perpendicular to the MER, is required to more thoroughly test if there are shallowly dipping melt inclusions extending away from the MER and to better constrain the extent of melt present.A current limitation is that most deployments at the MER are temporary and therefore often do not record a sufficiently large sample size of SKS phases for our stacking approach to be robust.
This example serves to highlight the potential of at-tenuation anisotropy to enhance our understanding of melt or fluid-rich regions, even where we have a good understanding of seismic anisotropy in the region.At a minimum attenuation anisotropy is potentially a useful tool for identifying the presence of fluids in the subsurface, even at very low volume fractions.More extensive, dense, measurements of shear-wave splitting and attenuation anisotropy may, in the future, allow for strong constraints to be placed on important properties such as the volume fraction of melt present and the orientation of the melt inclusions.Currently our source polarisation stacking method can only be readily applied for SKS and other coretransiting shear-wave phases.In these cases we can make the assumption that all phases sample the sample region of the upper mantle beneath the station and, for some stations, achieved sufficient source polarisation coverage to measure attenuation anisotropy.For other shear-wave phases, such as local or teleseismic S, this is not the case.To achieve the requisite source polarisation coverage would require data that most likely samples different regions of anisotropy which makes taking a station average unsuitable.This poses a particular challenge to measuring shear-wave attenuation anisotropy in the near surface, where the potential for attenuation anisotropy to improve characterisations of fluid-filled fracture systems could prove powerful.Removing the requirement for source polarisation stacking is, therefore, desirable and is a promising avenue for future research.

Conclusion
Seismic attenuation anisotropy is a phenomenon which can be efficiently generated by models of fluid-filled fractures, particularly a squirt flow model.This attenuation anisotropy has a clear theoretical and observable effect on measurements of shear-wave splitting.The effect of attenuation anisotropy on the frequency content of split-shear waves can be measured using an adaptation of existing instantaneous frequency matching methods (Matheney and Nowack, 1995).Using synthetic shear-wave examples and SKS phases recorded at FURI, Ethiopia, we show these effects and that we can measure attenuation anisotropy and retrieve the underlying shear-wave splitting parameters.To explain the observed attenuation anisotropy, where the fast shearwaves appear more attenuated than the slow shear waves in SKS phases, a squirt flow model (Chapman, 2003) is required.Even allowing for a lack of constraints on the rock physics parameters it is clear that this requires shallowly dipping (ca.40 • ) melt inclusions which strike perpendicular to the Main Ethiopia Rift.Whilst the modelled strike and dip of the melt inclusions is contrary to expectations from previous work (e.g., Ayele et al., 2004;Kendall et al., 2005;Bastow et al., 2010;Hammond et al., 2014), there is some potential correlation with low shear-wave velocity anomalies seen in recent tomographic models that extend away from the rift.These results highlight the power of attenuation anisotropy measurements as a blunt tool to detect the presence of aligned melt inclusions within the Earth.2022).This cross-section is approximately perpendicular to the Main Ethiopia Rift and passes through FURI.The location and start/end points of the section (A-A') are shown in Figure 10.Black vertical lines indicate the approximate location of the Main Ethiopia Rift along the cross section.To reveal anomalies in the upper mantle, the colour scale is clipped at 3.9 km s −1 which masks crustal features.For details of crustal features which can be seen in the tomography, readers should refer to Chambers et al. (2022).
With further instrumentation and improvement of rock physics constraints, it may be possible to constrain the properties of fluid-filled fractures at a range of length scales within the Earth.

Figure 1
Figure1Seismic velocity, quality factor and attenuation anisotropy (expressed in terms of ∆t * ) predicted by rock physics models for cracked, fluid-filled media which considers crack scattering(a,b,c Hudson, 1981), and that model poroelastic squirt flow (d,e,f)Chapman, 2003).Also shown is ∆t * predicted solely by velocity anisotropy effects, computed using an elastic tensor for Olivine(Abramson et al., 1997) and Q = 50.Results for other isotropic values of Q are shown in Supplementary Figure1.Velocity and attenuation for P (blue), S1 (orange) and S2 (green) are calculated by solving the Christoffel equation (see text for details) assuming a 50 km thick medium and a dominant frequency of 0.1 Hz.For the fluid-filled models we use an isotropic solid with velocities v P = 6.5 km s −1 and v S = 3.6 km s −1 which contains melt inclusions with v S = 2.7 km s −1 , ρ = 2700 kg m −3 .These parameters are chosen to be broadly consistent with previous effective medium modelling of melt-induced seismic anisotropy(Hammond et al., 2014).

Figure 2
Figure 2 Anisotropic attenuation, ∆t * , as a function of fracture length as predicted by both squirt flow (dashed line Chapman, 2003) and crack scattering (solid line Hudson, 1980) models.∆t * is calculated for a propagation angle θ = 70 • relative to the crack normal and a frequency of 0.1 Hz.

Figure 3
Figure 3 Example of Gabor wavelet synthetics used.(a) shows a synthetic where shear-wave splitting, with fast direction φ = 30 • and delay time δt = 1.5 s.(b) shows the synthetics in panel (a) where differential attenuation ∆t * = 1.0 s has been applied by applying a causal attenuation operator to the slow shear-wave.(c) shows the synthetic from panel (a) where a differential attenuation ∆t * = −1.0s has been applied by attenuating the fast shear-wave.All synthetics in this Figure are generated with a source polarisation of 70 • , a dominant frequency of 0.2 Hz and a sample rate of 50 ms.

Figure 4
Figure 4 Amplitude spectra of synthetic shear-waves as a function of component reference frame rotation.At each reference frame rotation angle, we calculate the amplitude spectra and instantaneous frequency (black line) for the apparent fast and slow shear-waves.The left column shows the frequency content for the synthetics shown in Figure3a, where φ f = 30 • , δt = 1.5 s and ∆t * = 0 s.The right column shows the frequency content for the synthetics shown in 3b, where an attenuation anisotropy of ∆t * = 1 s has been applied to the synthetics shown on the right.This ∆t * is applied before rotating the components.

Figure 5
Figure 5 Example of the effects of component rotation and attenuation anisotropy on the frequency content and shear-wave splitting (parameterised by λ 2 for a synthetic shear-wave.Instantaneous frequency (a), the difference in instantaneous frequency (b), and the second eigenvalue of the trace covariance matrix (c) measured for the synthetic shear-wave shown in Figure 3b over the range of reference frame rotations of the horizontal components (solid lines).At each reference frame rotation, we then correct for the differential attenuation by attenuating the apparent fast shear-wave (blue) by t * = 1 s and repeat the measurements (dashed lines).The solid vertical line shows the applied (or true) fast polarisation direction, 30 • and the dashed vertical line shows the fast polarisation direction that would be recovered if the synthetics are not corrected for ∆t * .

Figure 6
Figure 6 Example |∆f | grid search results for individual synthetic waveforms generated at different source polarisations.Synthetics are generated with shear wave splitting parameters φ f = 30 • , δt = 1.0 s, attenuation anisotropy ∆t * = 1 s and source polarisations of 45 • (a), 130 • (b) and 285 • (c).These |∆f | surfaces can then be stacked (d), with the minima of the stack returning the input attenuation anisotropy ∆t * and fast polarisation direction φ f .

Figure 7
Figure 7 Source polarisation weighted, stacked |∆f | surfaces.Each panel shows the |∆f | stack measured for 100 Gabor wavelet synthetics generated with shear-wave splitting parameters φ f = 30 • , δt = 1.0 s and ∆t * = 1 s (a) or ∆t * = −1 s (b).|∆f | is calculated for each synthetic by grid searching over φ r and ∆t * .The delay time δt is not measured at this point in the workflow as it does not affect |∆f | provided that a suitable analysis window has been chosen.Each synthetic is generated with a random source polarisation and with a dominant frequency drawn from f ∼ N (0.1, 0.02).

Figure 8
Figure 8 Bootstrapped summary statistics for the parameters φ r , (a) and ∆t * (b) along with the minimum |∆f | of each bootstrapped stack (c) for synthetic |∆f | stacking example shown in Figure 7a.The initial set of 100 individual |∆f | measurement grids is resampled, with replacement, 10,000 times and we repeat the stacking for each sample.The red vertical line in panel (c) indicates the bootstrap estimated 95% confidence level in |∆f |.

Figure 9
Figure 9Results of synthetic shear-wave splitting measurement stacking, following the method ofRestivo and Helffrich (1999).We generate 100 synthetics with shear-wave splitting parameters φ f = 30 • and δt = 1.5 s.Attenuation anisotropy of ∆t * = 1 s (a,b) or ∆t * = −1 s (c,d) is applied.Panels (a,c) show the shear-wave splitting results if we do not correct for this attenuation anisotropy.Panels (c,d) show the result after we correct the synthetic data using measurements of φ r , ∆t * made using our source polarisation stacking method.The stacked λ 2 surfaces are normalised by the 95% confidence value, indicated by the bold contours, which is derived from an F-test(Silver and Chan, 1991;Restivo and Helffrich, 1999).

Figure 10
Figure 10 Map showing shear-wave velocity beneath the Main Ethiopia Rift at a depth of 100km, obtained from the joint inversion of ambient noise and teleseismic Rayleigh waves (Chambers et al., 2022).Thick black lines indicate border faults and red polygons indicate magmatic segments.The location of FURI is shown by the yellow triangle.Station averaged SKS shear-wave splitting, after correcting for attenuation anisotropy, indicated by the black bar plotted on FURI, where the length of the bar corresponds to the delay time, δt, and its orientation to the fast polarisation direction.Measured attenuation anisotropy is shown by the magenta bar and follows the same plotting convention as the shear-wave splitting result.The cross-section A-A' (white line) through the tomography model is shown in Figure16.The inset map shows the locations of the 584 events used in this study (grey circles).From these events, we can identify 73 that yield clear SKS picks which are used to measure shear-wave splitting and attenuation anisotropy, shown by the red circles.We only consider events with an epicentral distance between 95 • and 110 • , the dashed lines mark the distance from FURI (yellow triangle) in intervals of 30 • .Event locations are taken from the International Seismological Centre (2023) bulletin.

Figure 11
Figure 11 One of the SKS phases recorded at FURI, Ethiopia, which we use to measure attenuation anisotropy.Panel (a) shows the pre-processed SKS phase rotated to the station averaged fast polarisation direction, φ = 38 • , that we measure for FURI and time shifted by δt = 1.15 s.The fast shear-wave, S1, is shown in blue and the slow shear-wave, S2, is shown in orange.Note that S1 appears to have a slightly longer period than S2, which suggests it has been more attenuated.Dashed lines show the measured instantanous frequency for the chosen analysis window (black lines).Panel (b) is plotted in the same style, showing the SKS phase after we correct for the measured attenuation anisotropy.

Figure 12
Figure 12 Source polarisation weighted, stacked |∆f | surface for FURI, Ethiopia.This result is obtained by stacking 73 |∆f | surfaces measured for SKS waveforms recorded at FURI.We measure an attenuation anisotropy of φ r = −45 ± 3 • and ∆t * = 0.45 ± 0.20 s, indicated by the blue cross.The 95% confidence region in our solution is demarcated by the bold contour and coloured white.Our approach to measuring |∆f | means that we cannot initially determine the sign of ∆t * .Upon analysis of our corrected SKS shear-wave splitting results (Figure 13b) we can determine that ∆t * = −0.45± 0.20 s.
Figure 14Results of modelling the fracture dip and medium thickness, or ray path length, which best explain the observed δt and ∆t * at FURI, Ethiopia using a squirt flow model(Chapman, 2003).Due to the lack of constraint on the mineral-scale relaxation time, τ m , we search over a range of τ m values for an assumed grain size of 1 mm and fracture lengths of 10 m (blue), 100 m (orange) and 1000 m (green).Panel (a) shows the fracture-scale relaxation time, τ f , which is proportional to τ m (23).The normalized least-square misfit of the best-fitting model for each τ m is shown in (b), with the predicted ∆t * and δt shown in (c) and (d).The observed ∆t * = −0.45s and δt = 1.6 s are shown by the solid black lines in (c) and (d), with the measurement uncertainties indicated by the shaded region.Panels (e) and (f) show the medium thickness, assuming a single anisotropic layer, and fracture dip angle required.

Figure 15
Figure 15Forward modelling, assuming a poroelastic squirt-flow model(Chapman, 2003): results for fracture strike and dip that can explain the observed attenuation anisotropy and shear-wave splitting in vertically incident SKS phases at FURI.We model a medium with 1% total melt volume fraction, of which 0.1% is hosted in aligned fractures.Medium thickness is fixed to 90 km following our previous models (Figure14e), which is consistent with previous estimates of the maximum thickness of anisotropy in the region(Ayele et al., 2004;Hammond et al., 2014).We use a normalised least squares cost function for ∆t, δt, and φ f .Our results favour shallowly dipping, 39 • , fractures which are oriented NW-SE, which is approximately perpendicular to the Main Ethiopian Rift.For details of other model parameters used see text.

Figure 16
Figure 16Interpolated cross-section through the shear-wave velocity model ofChambers et al. (2022).This cross-section is approximately perpendicular to the Main Ethiopia Rift and passes through FURI.The location and start/end points of the section (A-A') are shown in Figure10.Black vertical lines indicate the approximate location of the Main Ethiopia Rift along the cross section.To reveal anomalies in the upper mantle, the colour scale is clipped at 3.9 km s −1 which masks crustal features.For details of crustal features which can be seen in the tomography, readers should refer toChambers et al. (2022).

Table 1
Parameters used in squirt flow modelling of attenuation anisotropy observed at FURI, Ethiopia.For details of microscale relaxation time, τ m , grain size, ζ, and fracture length, a f , used see text.