Heterogeneous seismic anisotropy beneath Madeira and Canary archipelagos revealed by local and teleseismic shear wave splitting

observe a small average of S -wave splitting times of 0.16 ± 0.01 s, which signiﬁcantly increase with source depth beneath El Hierro ( > 20 km) and Tenerife ( > 38 km) up to 0.58 ± 0.01 and 0.47 ± 0.05 s. This suggests an inﬂuence of melt pocket orientation in magma reservoirs developed at uppermost-mantle depths. Likewise, anisotropy increases signiﬁcantly beneath the islands with shield stage volcanism (up to 9.81 ± 1.78 per cent at El Hierro, western Canaries, against values up to 1.76 ± 0.73 per cent at Lanzarote, eastern Canaries). On average, teleseismic SKS -wave splitting delay times are large (2.19 ± 0.05 s), indicating sublithospheric mantle ﬂow as the primary source for anisotropy in the region. In the Canaries, the western islands show signiﬁcantly smaller average SKS delay times (1.93 ± 0.07 s) than the eastern ones (2.25 ± 0.11 s), which could be explained by destructive interference above the mantle upwelling. Despite complex patterns of fast polarization directions throughout both regions, some azimuthal pattern across close stations can be observed and related to present-day mantle ﬂow and anisotropy frozen in the lithosphere since before 60 Ma. Additionally, we infer that the current presence of a mantle plume beneath the archipelagos leads to the associated complex, small-scale heterogeneous anisotropy observations

shear wave splitting of data from teleseismic events as seismic shear waves travelling in an anisotropic medium split in two orthogonal components with different speeds (Silver & Chan 1991).In case of multiple anisotropic layers, including measurements from local events allows to distinguish crustal from upper -antle contributions (Russo & Silver 1994).
Upper-mantle seismic anisotropy is thought to be mainly caused by dislocation creep-related lattice or crystallographic preferred orientation (LPO, CPO) of olivine due to mantle flow (e.g.Mainprice & Nicolas 1989;Karato & Wu 1993;Savage 1999;Skemer & Hansen 2016).The most abundant large-scale alignment of olivine under dry mantle conditions results in a crystallographic a-axis orientation in the mantle flow direction (A-type; Karato et al. 2008).Through the presence of water (Jung & Karato 2001) or change in pressure other alignments can occur (Mainprice et al. 2005).The asthenosphere can be viewed as a deformed low viscosity layer that accommodates the relative motion between the lithospheric plates and the underlying mantle flow.Current mantle movement continuously overprints anisotropy in the asthenosphere (e.g.Conrad & Behn 2010), making it a key probe of mantle flow (e.g.Kawakatsu & Utada 2017;Wang & Becker 2019).In particular, mantle plumes and plume-plate interactions affect the patterns of anisotropy into parabolic asthenospheric flow (e.g.Ribe & Christensen 1994;Walker et al. 2001Walker et al. , 2005;;Hammond et al. 2005;Collins et al. 2012;Ito et al. 2014), although laboratory experiments show that despite a simple geometric setup the splitting pattern can become more complex, partly facilitated by a tilted plume-head (Druken et al. 2013).Numerical simulations by Rümpker and Silver (2000) show that this pattern is also perturbed with the distance away from the central upwelling.On the other hand, in the highly viscous lithospheric mantle, the observed anisotropy is typically considered as 'locked-in', 'frozen' or fossil anisotropy and can thus be associated with past deformation events (e.g.Conrad & Behn 2010;Assumpc ¸ão et al. 2011;Vinnik et al. 2012).In the crust, anisotropy is mainly extrinsic due to the alignment on larger scale features (shape preferred orientation; SPO).These are caused by either intrusions, fractures or the fine ordered layering of multiple materials with different elastic properties but can be further influenced by the presence of fluids or temperature (e.g.Crampin & Booth 1985;Mainprice & Nicolas 1989).
Several shear wave splitting studies have been carried out in some hotspot regions such as, for example, in Hawaii (e.g.Walker et al. 2001;Collins et al. 2012), Iceland (e.g.Xue and Allen 2005), Réunion (e.g.Scholz et al. 2018), French Polynesia (e.g.Barruol et al. 2009), Galápagos (Fontaine et al. 2005), Eifel (Walker et al. 2005), East Africa (e.g., Walker et al. 2004), Cape Verde (Lodge and Helffrich, 2006) and the Seychelles (Hammond et al. 2005).However, many other hotspot regions have not been studied yet, such as the Canaries and Madeira region.Continental-based studies have shown that in nearby Iberia and north-western Morocco a mostly uniform pattern of fast shear wave polarization direction (FPD) is dominant (see Díaz et al. 2015, for references), mimicking the broad plate-driven mantle flow modelled by Conrad & Behn (2010).The exception is the area around the more complex Gibraltar Arc, where mantle flow is deflected around a slab, significantly diverting from global mantle flow models.
As part of the SIGHT project ('SeIsmic and Geochemical constraints on the Madeira HoTspot'), this study carries out the first detailed observation of seismic anisotropy beneath the Madeira and Canary Islands archipelagos to better constrain mid-plate upward mantle flow.We use teleseismic and local shear wave splitting measurements of data collected from 69 seismic stations located on all major islands of both regions.

S E T T I N G
The Madeira archipelago is located east of the Mid-Atlantic Ridge (MAR), around 700 km off the coast of Morocco (Fig. 1a).Madeira is the largest island (737 km 2 ), followed by Porto Santo Island (42 km 2 ), around 40 km to the northeast, and the Desertas Islands (14 km 2 ), roughly 17 km to the southeast.A bathymetric connection of shallow seafloor (around 200 m) can be observed between Madeira and the Desertas, whereas Porto Santo is separated from Madeira Island by a channel with seafloor depths reaching 2500 m.The lithosphere in the region has been dated to an age around 130-140 Ma (Verhoef et al. 1991) and is supposed to be >80 km thick, being characterized by high admittance values, or geoid to depth ratios (Cazenave et al. 1988).In Porto Santo, the youngest volcanism dates back to 10.2 Ma (Ferreira et al. 1988), whereas the oldest records of volcanic activity on Madeira have been dated to 7 Ma (Ramalho et al. 2015).This agrees with the absolute plate motion (APM) in no-net-rotation frame, which is 23.48 mm yr −1 in the 40.14 • direction around Madeira (Kreemer et al. 2014).Volcanism on the main island has been as recent as 6400 yr ago (Geldmacher et al. 2000).Considering this age, the hypothetical occurrence of volcanism in 1748 CE (Zbyszewski et al. 1975) and He isotope evidence for present-day mantle degassing (Amaral et al. 2017), Madeira cannot be considered volcanically extinct.
About 440 km south of Madeira, the Canary Islands intraplate archipelago consists of seven major islands.They are, from west to east, El Hierro, La Palma, La Gomera, Tenerife, Gran Canaria, Fuerteventura and Lanzarote (Fig. 1b).The trend of increasing island and volcanic ages follows the APM, which is 24.59 mm yr −1 in 43.29 • direction around the Canary Islands (Kreemer et al. 2014).The different stages of volcanism (Carracedo 1999;Geldmacher et al. 2005;Gottsmann et al. 2008) also cause a difference in general topography, being relatively flat among the older eastern islands and containing steeper, higher peaks at the younger western islands.Martinez-Arevalo et al. (2013), using Ps receiver functions, showed that the crust beneath the islands increases from about 11.5 to 12.5 km in the west (El Hierro, La Palma) to about 20-30 km in the east (Fuerteventura and Lanzarote).
The exact nature of the formation of the Madeira archipelago has been a longstanding debate, since the islands are located within the African plate away from any plate boundary.Similar to the Canary Islands, it has been thought to be the surface expression of a hotspot linked to a plume-like structure, the origin of which is still debated.The coexisting Madeira and Canary hotspot tracks have been reconstructed in a parallel manner to past locations close to the southwestern part of the Iberian Peninsula and north-western Africa, respectively (Hoernle et al. 1991;Mata et al. 1998;Geldmacher & Hoernle 2000;Geldmacher et al. 2000Geldmacher et al. , 2005)).A recent tomographic study (Civiero et al. 2021) has shown that the Madeira and Canary plumes are in distinct stages of evolution with the Madeira plume being only traced down to 300 km depth, in opposition to what is observed for the Canary plume (see also Civiero et al. 2018).Due to their close proximity, a genetic link between these mantle plumes has been proposed to varying degrees.Although varying geochemistry in Madeira and in the Canaries supports the concept of separate plumes, isotope ratios converge towards a common composition (e.g.Geldmacher et al. 2011).This suggests a common origin for a low S-wave velocity anomaly thought to exist in the upper mantle (Hoernle et al. 1995) or, alternatively, beneath the 660 km discontinuity and linked to the deepest levels of the mantle (Civiero et al. 2021;see also French & Romanowicz, 2015;Marignier et al. 2020).This is consistent with results from noble gas studies that also  S1, Supporting Information); stations that provided data but do not yield good results are shown in grey.Bathymetry and topography are taken from ETOPO1 (Amante & Eakins, 2009), age contours (in Ma) are given by dark grey lines (Müller et al. 2008).Temporary stations are shown in black.Black arrows indicate the average annual plate motion predicted at exactly those points, taken from GSRM v.2.1 (Kreemer et al. 2014).They are expected to be constant within uncertainties throughout both areas.In (c), the major plates are abbreviated to Af-African plate, EuAs-Eurasian plate, NAm-North American plate and SAm-South American plate.
suggest a lower mantle contribution for the Canaries (Day & Hilton 2011).Alternative models of the Canary Islands region include the existence of a propagating fracture linking the islands and the African Atlas Mountains, local rifting and compression-related tectonic uplift processes (e.g.Anguita and Hernán 2000).Moreover, interactions between a plume and small-scale edge-driven uppermantle convection have also been suggested (e.g.Geldmacher et al. 2005).

M E T H O D
For the analysis of seismic anisotropy, we use the approach by Silver & Chan (1991) to determine the two defining splitting parameters, the fast polarization direction, FPD (ϕ) and the time lag between fast and slow shear waves (δt), also known as delay time.This approach utilizes the characteristic elliptical particle motion that appears in cases where δt is smaller than the dominant seismic wavelength.The resulting ellipticity is then described using a covariance matrix.During the analysis a correction is introduced that removes the effects of the splitting by minimizing the second eigenvalue of the covariance matrix.A successful computation can easily be identified by a linearization, which reveals the polarization of the initial isotropic shear wave.To find a stable solution, we work with multiple analysis time windows around the onset via cluster analysis (Teanby et al. 2004).The absence of detectable splitting (called 'null' results) can be caused by the absence of anisotropy, but also if the backazimuth is equal to the FPD or FPD + 90 • .It can, therefore, be another important parameter in characterizing subsurface anisotropy.
In case of LPO of olivine the most abundant alignment is A-type fabric, commonly induced by strain within the upper mantle, which results in FPD parallel to the mantle flow.A change in water content or stress can change the type (Karato et al. 2008).The change in FPD is minimal close to vertical incidence angles but can deviate significantly above 20 • (Lynner et al. 2017).In case of SPO, the FPD is parallel to the alignment of cracks and fractures, which can be used to infer the direction of maximum horizontal stress.Results from teleseismic phases are likely to be dominated by LPO but show incidence angles close to the vertical.Results from local phases can have larger incidence angles, especially if events are shallow, in which case the azimuthal variation in SPO influence is large and can lead to SPO dominated results (e.g.Song & Kawakatsu 2012;Smith et al. 2017).
Following Schlaphorst et al. (2017), we define a range of criteria to estimate the quality individually by visual inspection.These are: a clear onset of the shear phase is visible in the data (Fig. 2a); high signal-to-noise ratio on the radial component (SNR ≥ 3; SNR ≥ 10 for null results); a significant amplitude reduction on the corrected transverse component, which hints towards a reduction of the second eigenvalue of the covariance matrix that describes the particle motion (Fig. 2b); a similar shape of the fast and slow shear (R) and transverse (T) components before and after ( * ) the correction.For a good splitting measurement, the energy on the transverse component should be minimized after the splitting correction.(c) Top row: fast (solid line) and slow (dashed line) shear waves before (left) and after correction (middle and right).The amplitudes are normalized in the left and middle diagrams, true amplitudes are shown in the right-hand panel.Bottom row: particle motion in the polarization plane within the defined window before (left) and after correction (right).Shear wave splitting tends to an elliptical particle motion, which should be linearized after the correction is applied.(d) Big panel: contour plot of the grid search over δt and ϕ showing the second eigenvalue of the covariance matrix.The best splitting parameters are indicated by the cross, the thicker black line indicates the 95 per cent confidence contour of the F test.Above the panel information about the splitting parameters (δt and ϕ), the source polarization (spol) and the backazimuth (baz) is shown.Note that the source polarization is close to the backazimuth (∼8 • and 11.9 • ).Small panels: best splitting solutions for the 100 time windows around SKS (top) and cluster analysis of the solutions (bottom).A good result will be stable over different time windows, creating plateaus in both splitting parameters, with the best solutions clustering in one location.In this case the plateau spans the entire range of windows and all solutions cluster in one location.Examples of a local S splitting and a null result can be found in Figs S1 and S2 (Supporting Information).
waves (Fig. 2c, top row); a notable change from elliptical to linearized particle motion due to the correction (Fig. 2c, bottom row); a clear minimum of the second eigenvalue amplitude with small well-constrained 95 per cent confidence ellipses on the error contour plot spanning the two splitting parameters (Fig. 2d, big panel); a stable suite of splitting parameters over changing time windows, creating a plateau (Fig. 2d, top small panel); a clearly identifiable cluster around the best solution with few minor secondary clusters (Fig. 2d, bottom small panel).
Solutions with large differences between polarization angle and backazimuth (>30 • ) are disregarded as they hint to potential interference with further anisotropic layers in the deep mantle (Hall et al. 2004).Using this approach, we strike a balance between ensuring minimal interference from potential D anisotropy and sufficient inclusion of path deviations due to lateral heterogeneity in the upper mantle.Since we do not have this constraint for local S phases, we impose a further quality control by excluding all results with uncertainties larger than ±0.1 s (2σ uncertainty).Based on a combination of minimal energy on the transverse component, resulting in linear particle motion, and a low-quality factor (Wüstefeld & Bokelmann 2007) measurements are categorized as 'null'.
In a further step we combine individual splitting teleseismic splitting measurements per station to increase result robustness.For parameters that are stable over the entire azimuthal range the combination can be calculated for the entire station.For parameters varying by backazimuth, the combination can be split into multiple  1, show the stations.Depth distribution of events is shown in the histogram, using a bin width of 2.5 km.The dashed box shows the enlarged view of El Hierro in (b).(c) Map showing the locations of all local events from the Madeira archipelago.Note that events are shown with a uniform circle diameter, since no magnitude information exists for this catalogue.(d) Azimuthal equidistant map showing the event locations of all teleseismic events used in this study (white stars for Madeira andd grey stars for Canary Islands).The dashed concentric circles represent the search distance of 85 • to 135 • to the study region (black triangle).Note that this is a combined plot for all stations using an average station position of 30.1 • N, 16 • W.
azimuthal sectors with stable parameters.Combing can be done by either stacking the entire error contour plot (Wolfe & Silver 1998;Restivo & Helffrich 1999) or by calculating a weighted average of the best solution of each individual splitting measurement (e.g.Kong et al. 2015).
Observing differences in splitting parameters between results from local S phases from events with different depths and teleseismic XKS phases (the most common for splitting studies being SKS, SKKS and PKS) helps to narrow down the depth ranges of main anisotropic layers.Local S rays will only cover the subsurface above the hypocentre, whereas teleseismic XKS rays travel through the entire mantle.Equally, a difference in δt between local S splitting results of events with different depths is a strong indicator of an anisotropic layer located between the two, although variation in apparent splitting parameters with azimuth especially at high incidence angles from shallower events with just one deeper layer is possible.Other constraints can come from the azimuthal variation of splitting measurements across single and/or multiple close stations, taking into account the Fresnel Zones of XKS waves and their potential overlap (Alsina & Snieder, 1995).From the local events we calculate the percentage of anisotropy, A, (Savage, 1999): with d being the length of ray path approximated by the hypocentre-receiver distance and v S being an average shear wave velocity.
in 2011/2012.In contrast, most stations on the Canary Islands, as well as two permanent stations in Madeira (PMOZ) and Porto Santo (PMPST) have been operating for a longer time period and most of them are still active.In total, the local networks provide 26 stations for the Madeira archipelago (of which 12 are short period) and 43 stations for the Canary Islands.
Teleseismic data are collected at distances between 85 • å and 135 • (Fig. 3c) to ensure that the phases are temporally isolated to be distinguishable from other incoming phases.The observations can be carried out on any XKS phase, but PKS phases have a generally lower SNR and we were only able to observe results from SKS phases.To increase the likelihood of sufficient SNR, we set a lower earthquake magnitude limit of M w 5.5.Before applying shear ave splitting corrections, the data are filtered using a zero-phase Butterworth bandpass filter with corner frequencies of 0.04 and 0.4 Hz, two poles and two passes.Even though this cannot suppress oceanic noise present around the islands it was found to show the best SNR values and highest overall number of successful measurements.
For both study regions most events arrive from a backazimuthal direction between 0 • (north) and 90 • (east).These areas of earthquake origin form predominantly strings along the Aleutian Arc through Japan, splitting into two branches along the Izu-Bonin-Mariana Arc and the Ryuku Arc to the Philippine trench, and connecting to the Java and Sunda trench (Fig. 3c).This results in approximately 8000 event-station pairs for Madeira and over 20 000 pairs for the Canary Islands.
Local data are chosen to be located close to the station and deep enough to allow for steep incidence angles, which prevents P-to-S reverberations that can be mistaken for shear wave splitting; the incidence angle value depends on the Poisson's ratio and, in this setting, can be approximated with an upper limit of 35 • from the vertical (Evans 1984;Savage 1999).For local data we set a minimum magnitude limit of 2.5.Due to the different frequency content of local waves and the generally low SNR we had to use different filter boundaries to the teleseismic events of 2 and 5 Hz.It has been shown in continental settings that frequency dependency of splitting has an influence on the results (e.g.Eakin & Long 2013), which affects the comparability of teleseismic and local phases.Although, the same is not necessarily found in oceanic environments (e.g.di Leo et al. 2012;Schlaphorst et al. 2017), the possibility cannot be ruled out.To test for frequency dependency, we processed local S phases in the teleseismic frequency limits.Where signals could be observed at those lower frequencies, no evidence of frequency dependent effects could be found.
Local seismicity around Madeira is sparse as revealed by the Instituto Português do Mar e da Atmosfera (see Data Availability) catalogue and also referred in Matos et al. (2015) and currently there is no extensive local seismicity catalogue.Local seismicity around the Canary Islands (National Geographic Institute, IGN, see Data Availability) that fulfills the search requirements is mostly concentrated around El Hierro (Fig. 3a), but smaller clusters of events exist around stations on the other islands (Fig. 3b).We find 584 suitable event-station combinations for the Canary archipelago and 255 for the Madeira archipelago, mostly located at depths between 10 and 50 km.

Teleseismic SKS
An example of a good teleseismic SKS splitting measurement obtained in this study can be seen in Fig. 2 (further examples are shown in Figs S1 and S2, Supporting Information).We have compared the best solution weighted averaging of SKS splitting results with full error surface weighted stacking and can confirm that in complex settings both methods tend to agree on the FPD, while stacking is significantly underestimating the delay time (see also Kong et al. 2015).Therefore, in this study we focus on the weighted averaged values using as weights the individual measurement uncertainty, but for completeness the stacked results can be seen in Fig. S3 (Supporting Information).All averaged splitting results are listed in Table 1 and shown in Figs 4 and 5. Local and teleseismic delay times and local A averaged values (eq. 1) over individual islands are listed in Table 2. Individual splitting measurements are listed in Table S2 and shown in Fig. S4 (Supporting Information).Of all available event-station pairs, 198 results around Madeira (∼2.5 per cent, of which 30 are null results) and 230 around the Canary Islands (∼1.1 per cent, of which 28 are null results) passed the quality control.The low return of good results is to be expected at stations located in ocean island settings due to high noise levels (e.g.Barruol & Ismail, 2001).
In both regions the delay times are generally larger than 1 s and often exceed 2 s, suggesting at least one major anisotropic layer beneath the crust.However, a crustal contribution cannot be ruled out.
In Madeira, changes in both splitting parameters occur on very small length scales (Fig. 4).Still, patterns spanning multiple stations in the archipelago can be observed.Most events arrive at the stations from north-eastern directions (NE events).Stations in the central part of the island have FPD alignment close to the APM direction for NNE events, whereas stations to the east and west of the island centre show predominantly perpendicular alignment.In general, most FPDs are subparallel to the APM but variations are complex at many individual stations across the region.However, some stations in central Madeira show uniform patterns subparallel to the APM (DOC-01, DOC-05, DOC-16; see Fig. 6 and Fig. S4, Supporting Information).Interestingly, splitting measurements from all south events show a uniform alignment close to the APM, whereas all SW events show either oblique or perpendicular alignment.In contrast, the station in Porto Santo, though having azimuthal variations as well, shows alignment almost perpendicular to the plate motion, although the average delay time is smaller.
Similar to Madeira, we can observe azimuthal variations within single stations or station clusters in the Canary Islands (Fig. 5), although unlike Madeira, many patterns seem to be consistent over entire islands, showing less variation in individual sectors, especially on El Hierro and Gran Canaria.We observe differences between the eastern, central and western islands but in comparison to the Madeira archipelago the variation is more pronounced over the extent of the Canary Islands.In general, the western islands show significantly smaller average delay times.Even though the patterns are azimuthally variable on all islands, they tend to resemble uniform patterns more closely on Gran Canaria and the eastern islands, either almost perpendicular to the APM on Gran Canaria and Fuerteventura, or predominantly subparallel on Lanzarote.On Tenerife and the western islands, the patterns are more varied, although a majority contains subparallel FPDs.
On islands where different stations show significantly different splitting results in the same sectors, clusters of stations can be observed.On La Palma, for N events, only the central stations show perpendicular FPDs.Similarly, on Tenerife for NE and E events, predominantly the central stations and stations to the northeast of the island show FPDs perpendicular to the APM.However, while on La Palma all SW events show similar perpendicular FPDs, all south events in Tenerife show subparallel to oblique FPDs.Barruol & Ismail (2001).f Behn et al. (2004).In general, we do not concentrate on individual stations due to the similarities of close stations.However, we can compare our results to previous measurements taken at station TBT on La Palma.Our results indicate two sectors with FPDs of 36.2 ± 6.0 • (NE) and −18.5 ± 3.1 • (SW) when considering similarly aligned stations on La Palma (Table 1).Averaging only TBT we find a value of 21.1 ± 4.1 • (NNE), which is in agreement with previous results (19 ± 12 • - Barruol & Ismail, 2001; 17 ± 5 • - Behn et al. 2004).However, the variation among different azimuthal sectors in our results suggests that the subsurface is more complex than a simple one-layer approximation (Fig. 6).
The distribution of null results shows a few predominant azimuthal directions in both archipelagos (Fig. S5, Supporting Infoormation).In Madeira there are two at around −15 • and +80 • and in the Canary Islands the largest group is at around 0 • .However, these clusters are partly caused by multiple stations recording null results from the same event in combination with the fact that we recorded relatively few null results with confidence.

Local S
Local S splitting results shown in Fig. 7 and listed in Table S2 (Supporting Information).Of all suitable event-station pairs, a comparably large number of 241 pairs (∼41.3 per cent) show results that passed the quality control at the Canary archipelago, whereas around Madeira we find 14 pairs (∼5.5 per cent).Except for one event, all of these are located shallower than 50 km.Where a comparison is possible, results from events with incidence angles below 20 • are similar to events with higher incidence angles, suggesting no change in olivine fabric types.However, it is possible that other types can result in similar apparent splitting parameters due to their azimuthal variation at high incidence angles.Apart from El Hierro, Tenerife and Madeira, the sparse coverage of results prevents us from putting individual measurements in context, thus preventing further detailed analysis.However, a few general observations can be made about the other islands.
FPDs on Madeira are predominantly oriented in east-west direction.In the Canary archipelago, the FPD on most islands follows uniform patterns of roughly +30 • (NE) on Gran Canaria and Fuerteventura, together with most of Tenerife and northeastern El Hierro.More details can be observed where denser b Number of events.We do not stack all individual measurements but rather average the results of each station, weighed by their number of stacked events, since small-scale heterogeneities can cause anisotropic effects to cancel out, resulting in arbitrarily small delay times.c Number of averaged events.d El Hierro, La Palma, La Gomera.e Tenerife, Gran Canaria.f Fuerteventura, Lanzarote.clusters of results emerge (Fig. 7).In Tenerife, a southeastern region shows multiple events, consistently having among the largest delay times in the area with values of over 0.5 s.Adjacent in the north is a cluster of multiple results with similar FPDs but among the smallest delay times in the area.In the south of the island, a cluster of results breaks the pattern of otherwise uniform FPD.
El Hierro, together with La Palma being the currently most volcanically active part of the Canary Islands, witnesses the largest number of earthquakes and, thus, provides the most comprehensive set of splitting results.A larger number of results show bigger delay times of up to 0.5 s.The FPDs vary across the island (Fig. 7b), however three major clusters can be identified: (1) located offshore to the western tip of the island with a lower number of events ( 16) mostly at around 30 km depth; (2) along the southern coast with more events (90) with a depth range of 10-40 km and (3) in the northeast of the island with the largest number of events (316) at around 20 km.Cluster 1 shows two almost equally large groups of FPDs in +60 • or −60 • direction and average delay times of around 0.25 s.Cluster 2 shows FPDs shifting gradually from west to east along the coast with values of around −15 • to +15 • , keeping an average delay time of around 0.25 s.Cluster 3 shows the strongest parallel alignment of predominant FPD at around +30 • and average delay times below 0.5 s but larger than those of the other two clusters.Results at the southern end of the cluster, closer to the centre of the island, show larger delay times of over 0.5 s and FPDs closer to 0 • .In combination, the clusters show a radial pattern with its centre on the island.All results are plotted on the midpoint between the event and station location, which due to the curvature of the ray path will be further away from the station than the midpoint of the path.For local waves this effect can be neglected.Rose diagrams show the summation of all FPDs for individual islands.Due to the abundance of results the map is separated into three major clusters on El Hierro (1-west; 2-south and 3-northeast).The rose diagrams are normalized to an angular bin size of 30 • .Note that their respective sizes are not linearly representative of overall number of events; all smaller rose plots have been enlarged to ensure visibility.

Delay time and percentage of anisotropy comparison
A closer look at the delay times reveals generally significantly smaller values for local S splitting observations; they are all below 0.6 s and about half of all results show values smaller than 0.15 s (Fig. 8).The average delay time on Madeira Island is 0.15 ± 0.03 s, similar to the result of the entire Canary archipelago of 0.16 ± 0.01 s.However, across the extent of the archipelago the results vary, with the largest delay times to be found on the western and central islands (0.17 ± 0.01 and 0.18 ± 0.02 s), decreasing to the east (0.08 ± 0.01 s).Our results show that, especially for local shear wave splitting, relative measurements uncertainties in delay time tend to be much smaller than relative uncertainties in FPD.Although the effect is still visible for teleseismic measurements, it is less pronounced.We find generally very low delay time uncertainties around ±0.01 s for over half of the results.Around El Hierro, delay times with values above 0.25 s can be observed only for events with sources beneath 20 km.Similarly, around Tenerife this can be observed at greater depths of around 38 km.This trend is absent around Gran Canaria but due to the low number of data points the pattern is more ambiguous.
The distribution of A (eq. 1), using an average value of v S = 4.0 ± 0.5 km s −1 , follows the trend observed in delay times (Table 2 and Fig. S6, Supporting Information).Highest values of up to 9.81 ± 1.78 per cent (with an average of 1.75 ± 0.07 per cent) can be found on El Hierro.The upper limit decreases for La Palma (4.42 ± 0.83 per cent) and Tenerife (4.83 ± 0.93 per cent), and decreases further for Gran Canaria (0.62 ± 0.14 per cent), Fuerteventura (1.65 ± 0.36 per cent) and Lanzarote (1.77 ± 0.73 per cent).On average, the central islands have a much lower anisotropy (0.67 ± 0.14 per cent) compared to the western (1.74 ± 0.07 per cent) and eastern islands (1.45 ± 0.12 per cent).Madeira shows an average anisotropy comparable to the western and eastern Canary Islands (1.25 ± 0.31 per cent).However, due to the lack of data we cannot make definite statements for any islands apart from El Hierro, Tenerife and Madeira.
In contrast to local events, averaged teleseismic observations show values larger than the global median of 1.0 s (updated from  Becker et al. 2012) and nearly half of them surpass 2.0 s (Fig. 8 and Table 1).On average the value is 2.19 ± 0.05 s, but the Madeira archipelago shows larger delay times (2.29 ± 0.07 s) than the Canary archipelago (2.09 ± 0.06 s).However, the large lateral extent of the Canary Islands facilitates a significant delay time difference, showing smaller values at the stations on the western islands (1.93 ± 0.07 s) compared to stations on the eastern islands (2.25 ± 0.11 s).

Multiple layers of anisotropy
The significant variation of apparent splitting measurements with backazimuth strongly suggests a complex subsurface with multiple (potentially inclined) layers of anisotropy (Fig. 9).Even though azimuthal coverage of results would be sufficient in both archipelagos when combining adjacent stations, modelling of apparent splitting results using either one, two or inclined layers of anisotropy, proves to be difficult due to the complexity of the results.A π /2-periodicity might be present but is generally masked by smaller scale variations.In addition, null measurements do not seem to separate clearly from splitting observations azimuthally in individual sectors (Fig. 6) and are, therefore, complicated to integrate into our models.It is possible that an SNR close to our lower limit can be responsible for null measurements.
For each island we tested different two-layer models, using a priori information from the present-day as well as 60 Ma APM, and rift zones.These tests show that in case of a dominating bottom layer, the apparent splitting parameters will stay close to the values of that layer for the majority of the azimuthal range, but depending on the characteristics of overlying anisotropic layers can vary rapidly over small azimuthal changes.Further layers closer to the surface, which are more likely to be susceptible to short-scale regional heterogeneities, can then superimpose constructive or destructive interferences, further complicating the results.This is a very likely scenario in both our regions.

D I S C U S S I O N
In this study, we used a combination of dense recordings of local and teleseismic shear waves to map seismic anisotropy in the Madeira and Canary archipelagos.These allow to distinguish the effects of shallower (uppermost mantle to crustal) from deeper (infra-to sublithospheric) anisotropies to the observed splitting seismic patterns.
In addition, we can introduce limits of minimum or maximum depth for anisotropic layers based on the overlap of Fresnel Zones (Alsina & Snieder 1995).For epicentral distances between 85 • and 135 • the Fresnel Zones have radii of approximately 13 km at 10 km depth and 70 km at 200 km depth for SKS waves with a dominant frequency of 0.1 Hz.Consequently, deeper Fresnel Zones will overlap for stations on the same island in both archipelagos (Fig. S8, Supporting Information).Therefore, azimuthal changes at the same or close stations will impose an upper limit, while variations in the same azimuthal range over multiple stations will impose a lower limit.However, local heterogeneities of anisotropic features close to the surface can have an effect that varies from one station Here, the fast shear wave orientation of the lower layer is based on the current plate motion, whereas the upper layer is based on the general orientation of faults and rift zones (e.g.Carracedo 1999).The delay times are chosen to match the average delay times of the individual islands, assuming a dominating mantle anisotropic layer.Further islands are shown in Fig. S7 (Supporting Information).
to the next or on a single station azimuthal range.The results of local events indicate that stronger anisotropic layers that would lead to SKS delay times that we observe in both archipelagos have to be located deeper.

Uppermost mantle to crustal anisotropies
Our observations suggest crustal and upper lithospheric mantle anisotropic contribution.This can be observed especially beneath the western, recently volcanically active Canary Islands, as evidenced by significant increase in shear wave splitting delay times with event depth depicted by local shallow events.Translating these delay times to percentage of anisotropy, A, a clear pattern emerges at the Canary Islands related to the structural difference between the central islands and the eastern and western ones (e.g.Gottsmann et al. 2008, and references therein).The central islands (Tenerife, Gran Canaria) show significantly lower values.El Hierro, which is also currently experiencing uplift fed by a recharging of a crustal magma chamber in the centre of the island (González et al. 2013), shows the highest A values (apart from Lanzarote, which is only based on one result, leading to a large uncertainty).However, due to the lack of data we cannot make definite statements for any islands apart from El Hierro and Tenerife.
For El Hierro, the increased delay time only observed for events with a source deeper than 20 km suggests the existence of an anisotropic layer at a depth of around 18-20 km.Uncertainties in this value stem from the non-uniqueness of the problem due to the sparseness of the observational constraints.This depth is a few kilometres below the Moho (∼11.5 to 16 km, Watts, 1994;Martinez-Arevalo et al. 2013).In a study observing different hotspot islands, Park & Rye (2019) proposed that crustal fractures allow seawater to infiltrate, prompting serpentinization of the uppermost mantle and inducing textural anisotropy through the formation of serpentine mesh networks.Alternatively, the observed anisotropy could be caused by the existence of a sub-Moho plumbing system, which has been considered important in creating strong anisotropy beneath volcanoes (e.g.Magee et al. 2018), yet the type of influence on shear wave splitting can be variable (e.g.Miller & Savage 2001).The fact that similar anisotropic layers are absent in the older Canary Islands renders the latter option more probable.This hypothesis for the 18-20 km depth anisotropy at El Hierro is endorsed by the results of a low-frequency microseismic sounding study (Gorbatikov et al. 2013) showing a magmatic reservoir at the same depth, which is likely to feed the smaller crustal reservoirs detected by González et al. (2013).
A similar if somewhat less pronounced feature can be observed on Tenerife.The bulk of the shallow anisotropy seems to be located at two depths as shown by two increases in delay times around those depths (Fig. 8).One can be found around 36-38 km and potentially a shallower one at a depth of 18-20 km.However, this is a hypothesis, since the data are sparse and a reduction in delay times can also be the result of differences in incidence angle or azimuth and the symmetry in anisotropy.The shallower one is located closely beneath the Moho, estimated to be between 15 and 18 km (Lodge et al. 2012, Martinez-Arevalo et al. 2013).Since smaller delay times can be observed below that depth as well, these features seem to be localized.Although we are aware that the non-uniqueness of the problem prevents a definitive interpretation, crustal underplating mechanism akin to the one suggested beneath El Hierro could explain an anisotropic layer at that depth.This suggests that future local shear wave splitting analyses may help detecting underplating in other settings.In general, however, sub-Moho earthquakes are a rare occurrence.
Petrological studies have shown the development of plumbing systems characterized by magma reservoirs located at mantle depths to be common at oceanic volcanoes with low magma supply rates, such as the Canary Islands (e.g.Longpré et al. 2008), Madeira (e.g.Klügel & Klein, 2006), or the Cape Verdes (e.g.Mata et al. 2017).For El Hierro, the anisotropic layer observed at circa 20 km depth is compatible with the existence of a magma reservoir at 19-26 km depth, inferred from petrobarometric data (Stroncik et al. 2009).At Tenerife, results point to magma reservoirs at approximate depths between 20 and 45 km, which puts them directly beneath the Moho and at the base of the long-term elastic lithosphere, respectively (Longpré et al. 2008).Lodge et al. (2012) noted that in the Canary and Cape Verde archipelagos the thickness of underplating processes tends to increase with island age.Our data show that anisotropic layers are better defined at El Hierro than in Tenerife, while on Gran Canaria, an island at a more advanced stage of its life cycle (Carracedo 1999;Geldmacher et al. 2005), no detectable increase in delay time with depth is observed, hinting towards an absence of such lower crustal/shallow mantle anisotropic layer.Given the active volcano setting of our study region, this suggests that at the youngest islands the anisotropy is caused by aligned melt pockets rather than stress changes imposed by magma injection overpressure (Miller & Savage, 2001), although the influence of micro cracks cannot be ruled out.The role of magmatic intrusion to the endogenous growth of ocean islands has been considered significant (e.g.Klügel et al. 2005).However, our data suggest that their influence on anisotropy is only evident while melt exists.The presence of melt beneath El Hierro is supported by the recent eruption of 2011-2012, as well as the observation of two large intrusive bodies beneath the island, reaching a depth of 35 km which seem to be followed deeper (Gorbatikov et al. 2013).
On El Hierro, the radial pattern of local FPDs (Fig. 7b) matches the general orientation of the triple-rift system described by Carracedo (1999), González et al. (2013) and Becerril et al. (2015), suggesting that these rift zones impose part of the seismic anisotropy observed.Furthermore, the entire island experiences uplift fed by a recharging of a crustal magma chamber in the centre of the island (González et al. 2013), which also affects the crustal stress field and, thus, the anisotropic patterns around El Hierro.

Intra-to sublithospheric anisotropies
Due to the difference in delay time between local S and teleseismic SKS splitting measurements on all islands of the Canary archipelago (Fig. 8), we conclude that the majority of anisotropy in the entire region is deep in the lithosphere or beneath it.This agrees with previous observations made by Behn et al. (2004) on 13 ocean islands surrounding Africa.Around Madeira, the lack of local S results prevents us from making definitive statements.However, due to the similarity in splitting parameters in both regions, including rather large delay times observed at most stations on Madeira that are difficult to explain with solely lithospheric contributions, we also see here strong evidence for a primary contribution from asthenospheric mantle flow.
On Madeira Island all south and southwestern events show uniform SKS splitting parameters across all stations, making a deep source such as a sublithospheric mantle flow likely.The area covered by the Fresnel Zones can be located south and southwest off the island.The northern events can be sorted into two groups of stations with one mostly located in the centre of the island and the other located to the eastern and western sides.The Fresnel Zones will start to overlap at depth just like in the previous example.Therefore, it is likely that an additional source of anisotropy at shallower depth (likely the shallow mantle) is present for one of these groups, which could be caused by disturbances by mantle upwellings (e.g.Civiero et al. 2018).However, an exact prediction is not possible because of the lack of knowledge of the shape of the upwelling structure due to resolution limitations of large-scale seismic tomography studies.The eastern events in contrast experience more variability, which can be caused by their ray paths traversing a much larger portion of the subsurface close to the island and the underlying hotspot (cf.Civiero et al. 2021), as evaluated by the Fresnel Zones and based on the general shape of east-west elongation of the island.
On La Palma, southwestern events have the same FPD, suggesting a deep anisotropic layer, whereas the northern events are split into different groups of distinctively different FPDs.Here, like in Madeira, it is likely that an additional anisotropic layer at a shallower depth has some contribution.
For both archipelagos, there is a notable difference in SKS splitting parameters between eastern and western islands, though much more pronounced along the larger lateral extent of the Canary Islands (Figs 4 and 5;and Table 2).In both regions, this coincides with an increase in island as well as plate age and lithospheric thickness towards the east (Fig. 1).Beneath the western islands, deeper anisotropic contributions caused by sublithospheric mantle flow are likely to be dominant, whereas the older eastern islands could experience interference of the present-day sublithospheric mantle flow with a stronger layer of fossil anisotropy in the lithosphere (e.g.Silver & Savage 1994;Barruol & Hoffmann 1999).The stronger component of fossil anisotropy in the lithosphere, distinct from the present-day mantle flow, can explain the apparent counter clockwise mismatch of around 45 • -90 • between predominant FPD around the older islands of Porto Santo (Madeira archipelago), Fuerteventura and Gran Canaria (Canary archipelago) to the present-day APM (+40 • ).Indeed, the direction of the African plate motion, at least from 90 to 60 Ma, is estimated at approximately −19 • with a speed surpassing present-day motion (Müller et al. 2019).It matches the FPDs observed in the eastern Canary Islands and in Porto Santo (Figs 4 and 5), suggesting the attachment at the base of the lithosphere of a significantly thick layer (Becker & Lebedev, 2019).In contrast, the plate motion between 60 Ma and the present has been significantly slower, likely hindering large-scale olivine crystal reorientation.Previous studies of hotspot regions such as, for example, Hawaii (e.g.Collins et al. 2012) and Réunion (Scholz et al. 2018) also reported FPDs consistent with the palaeo-orientation of spreading ridges.Although a classic twolayer model with strong anisotropic layers in the lithosphere and the asthenosphere could be approximated (e.g.Lanzarote, Fig. 9), azimuthal complexity suggests additional features, which we are not able to resolve with the data and station coverage, as well as the a priori information needed to establish a model due to the nonuniqueness of the approach.
In contrast, many sectors in Madeira Island and the western younger Canary Islands (El Hierro, La Palma, La Gomera and Tenerife) show an FPD close to the present-day mantle flow direction (Figs 4 and 5).Still, the delay times and FPDs have significant variation.However, no obvious correlation with backazimuth is observed, such as a π /2-periodicity (Fig. 9), which would hint towards two major anisotropic layers with horizontal symmetry axis (Silver & Savage, 1994;Savage 1999).Most likely the results show a more complex pattern with small-scale heterogeneities, which can be caused by diverted lateral mantle flow due to a strong vertical upwelling flow component, for example due to a plume beneath the islands.There are other explanations for the observed patterns, such as small-scale flow or a thin anisotropic layer influencing the FPD (Kaviani et al. 2013).Especially in regions of high complexity there is a higher chance of multiple anisotropic layers interfering destructively with each other.Therefore, significantly smaller delay times are a likely effect around the western Canary Islands, which are located above the active part of the hotspot (Taylor et al. 2020;Civiero et al. 2021).Due to the steep incidence angles of SKS of normally less than 10 • to the vertical, we do not expect to observe large contributions from vertically oriented structures (Barruol & Hoffmann, 1999), resulting in reduced delay times.Indeed, our results show significantly lower delay times around stations on El Hierro and La Palma.However, a plume is likely to consist in part of steeply inclined mass movement, with significant structural heterogeneities over small length scales (Schwarz et al. 2004;Weis et al. 2011), such as those created by multidirectional spreading when encountering the base of lithosphere.These are able to distort measurements into a complex pattern (Savage 1999).While large vertical motion of the plume, considered in isolation, will create null measurements, the measured splitting at the station is always a result of the combination of any anisotropic influence on its path.To account for every effect, a very good understanding of the position and structure of the plume would be required, which is beyond the scope of this study.Still, we note that the number of null measurements is not high in this study, which is a result of the strict quality control we had to impose due to the noisiness of the data.
A small-scale change of the infinite strain axis (ISA) can be observed around the Canary Islands, resulting from a plume-like feature in the model (Fig. 10).Focused ascending mantle material is a well-known feature to disrupt the main flow field, therefore affecting the anisotropy pattern (e.g.Ito et al. 2014).Some early studies relating plume dynamics and seismic splitting suggested that the shape of the pattern is strongly dependent on the ratio between plume volumetric flow rate and APM.It will become approximately parabolic when the plate motion is comparatively significant, for example around the Hawaii hotspot and can be distorted due to an additional amount of differently oriented fossil anisotropy (Ribe & Christensen 1994;Walker et al. 2001Walker et al. , 2005;;Ito et al. 2014).More recent, sophisticated geodynamical and mantle fabrics calculations showed that the predicted splitting pattern over intraplate plumes can be more complex, for example, forming nested U shapes (Ito et al. 2014).Our observations show further complexity and an apparent randomization of the flow field due to the mantle plume, which are hinted towards in laboratory experiments (Druken et al. 2013) and might be explained by future fully dynamical models of plume-plate interactions.The change to FPDs perpendicular to the general ISA in the Canaries region is confirmed by the small-scale FPD changes in our results.With our results we can enhance the global ISA with a more detailed picture of the pattern around our study area.While our results support complex ISA patterns in association with the plume, the observed highly scattered splitting is more complex than existing simple parabolic models.This suggests that future, more sophisticated models of plume-lithosphere interactions emulating the conditions of our study region are needed.It is possible that anisotropy can be located shallower than the 225 km displayed by the ISA (Fig. 10), but it must be beneath the hypocentres of local seismicity.Viewed in the larger regional context, the changes in FPD on small length scales stand in stark contrast to the broadly uniform pattern among the majority of land station observations in Iberia and Morocco (Fig. 10), with the exception of the area around the Gibraltar Arc.There, the subsurface structure is more complex, with mantle flow being deflected around the slab and significantly diverting from the global mantle flow modelled by Conrad & Behn (2010).We note that due to the station coverage limitation our results alone are not sufficient to demonstrate two separated mantle upwellings beneath the Madeira and Canary archipelagos but can be used in combination with results from geochemistry (Geldmacher et al. 2001;Civiero et al. 2018Civiero et al. , 2021)).It is also important to note that although signals from plume-related flow, as well as modern plate motions will combine with existing anisotropic patterns, the strength and complexity of the result would be different.A relatively slow-moving plate needs a much longer time to reorient large patterns, so the overwriting process could be slower than with a plume.The combination of broadscale orientation from plate motion and broadscale existing anisotropic signals will result in broadscale features, whereas the introduction of a plume would result in smaller-scale features that are combined with broadscale patterns, thus losing their broadscale appearance.
In contrast, around Madeira the global mantle flow model of Conrad & Behn (2010) does not show any disturbance in ISA (Fig. 10) because it does not consider the existence of a mantle plume.However, our averaged teleseismic SKS splitting FPDs divert significantly from that direction in a similar pattern to the one observed around the western part of the Canary Islands.This, in combination with the geochemical (e.g.Mata et al. 1998;Geldmacher et al. 2001) and seismic tomographic studies (Civiero et al. 2018(Civiero et al. , 2021)), supports the presence of a plume-like feature beneath Madeira, with similar processes to those occurring in the Canary Islands hotspot.Our new observations should help reassess mantle flow models of the region and update the global hotspot catalogue.

. C O N C L U S I O N S
In this study, we used a combination of local and teleseismic shear waves to map seismic anisotropy in the Madeira and Canary archipelagos.The results show complex behaviour that cannot be explained by only one or two uniform anisotropic layers.To better resolve the additional complexity a denser seismic network and larger number of events would be needed, which is beyond the scope of this study but could be possible with a longer station operation time and additional stations (e.g. from the UPFLOW OBS network; UPFLOW 2021).Still, we can attribute sublithospheric mantle flow as an important cause of anisotropy in the study region.Uppermost mantle anisotropies are evidenced by local seismic data revealing the existence of melt associated with crustal underplating at El Hierro (20 km depth) and Tenerife (18 and 38 km depth), at depths consistent with petrobarometric data indicating the presence of magma reservoirs.Furthermore, we observe strong indications of mantle flow perturbance, most likely from vertical movement confined to smaller upwellings within the regions, which is additional evidence for the existence of a plume beneath the Canary archipelago, notably below the westernmost islands.Similar observations around Madeira lead to the conclusion that similar mechanisms are active in that area, emphasizing the necessity of a reassessment of the mantle flow models for this region of the Eastern Atlantic.More generally, our work highlights the small-scale complexity of the patterns of mantle flow in hotspot regions and how anisotropy is a powerful tool to pinpoint the location of plume heads and to unravel elusive phenomena such as underplating.

A C K N O W L E D G M E N T S
The project was conceptualized by DS, GS and JM.Data acquisition was carried out by DS, GS, FK and TD.The data were processed and analysed by DS.The original draft was written by DS with reviewing and editing being carried out by all authors.
The authors would like to thank J-Michael Kendall, Ricardo Ramalho, Andréa Tommasi and the anonymous reviewers for their comments and suggestions.We also thank Manuele Faccenda for fruitful discussions.The data presented in this study were obtained through a project proposed and led by the National Laboratory of Energy and Geology (LNEG) to the Electric Company of Madeira (EEM) to assess the geothermal potential of the Madeira Island.This project benefitted from the framework of Deep Ocean Test Array (DOCTAR) project fieldwork during the data collection.Data for the permanent stations in Madeira were gathered with the help of the Instituto Português do Mar e da Atmosfera (IPMA).We are also grateful to Resurrección Antón and the Instituto Geográfico Nacional (IGN) for the Canary Islands data.This is a contribution to project SIGHT (ref.PTDC/CTA-GEF/30264/2017).The authors would like to acknowledge the financial support Fundac ¸ão para a Ciência e a Tecnologia (FCT), I.P./MCTES through national funds (PIDDAC)" -UIDB/50019/2020 -IDL.This study is also part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research

Figure 1 .
Figure 1.Map of the study regions, Madeira (a) and the Canary Islands (b), as well as their global location (c).All seismic three-component broad-band and short-period stations that provided data for this study are shown by triangles on the map (also listed in TableS1, Supporting Information); stations that provided data but do not yield good results are shown in grey.Bathymetry and topography are taken from ETOPO1(Amante & Eakins, 2009), age contours (in Ma) are given by dark grey lines(Müller et al. 2008).Temporary stations are shown in black.Black arrows indicate the average annual plate motion predicted at exactly those points, taken from GSRM v.2.1(Kreemer et al. 2014).They are expected to be constant within uncertainties throughout both areas.In (c), the major plates are abbreviated to Af-African plate, EuAs-Eurasian plate, NAm-North American plate and SAm-South American plate.

Figure 2 .
Figure 2. Example of good SKS-splitting measurement recorded at station DOC12 (Deserta Grande, see more details in Table S1, Supporting Information).(a) Filtered three components (E, N, Z).Solid vertical bars indicate the start (A) and end time (F) of the window used to isolate the SKS wave.Note that 100 different windows were tested and this one gave the best result.The dashed vertical bars indicate theoretical arrival times of different phases.(b) Radial (R) and transverse (T) components before and after ( * ) the correction.For a good splitting measurement, the energy on the transverse component should be minimized after the splitting correction.(c) Top row: fast (solid line) and slow (dashed line) shear waves before (left) and after correction (middle and right).The amplitudes are normalized in the left and middle diagrams, true amplitudes are shown in the right-hand panel.Bottom row: particle motion in the polarization plane within the defined window before (left) and after correction (right).Shear wave splitting tends to an elliptical particle motion, which should be linearized after the correction is applied.(d) Big panel: contour plot of the grid search over δt and ϕ showing the second eigenvalue of the covariance matrix.The best splitting parameters are indicated by the cross, the thicker black line indicates the 95 per cent confidence contour of the F test.Above the panel information about the splitting parameters (δt and ϕ), the source polarization (spol) and the backazimuth (baz) is shown.Note that the source polarization is close to the backazimuth (∼8 • and 11.9 • ).Small panels: best splitting solutions for the 100 time windows around SKS (top) and cluster analysis of the solutions (bottom).A good result will be stable over different time windows, creating plateaus in both splitting parameters, with the best solutions clustering in one location.In this case the plateau spans the entire range of windows and all solutions cluster in one location.Examples of a local S splitting and a null result can be found in Figs S1 and S2 (Supporting Information).

Figure 3 .
Figure 3. Locations of events used in this study.(a) Map showing the locations of all local events from the Canary Islands used in this study.Circle diameters are scaled by the magnitudes of the events.Triangles, like in Fig. 1, show the stations.Depth distribution of events is shown in the histogram, using a bin width of 2.5 km.The dashed box shows the enlarged view of El Hierro in (b).(c) Map showing the locations of all local events from the Madeira archipelago.Note that events are shown with a uniform circle diameter, since no magnitude information exists for this catalogue.(d) Azimuthal equidistant map showing the event locations of all teleseismic events used in this study (white stars for Madeira andd grey stars for Canary Islands).The dashed concentric circles represent the search distance of 85 • to 135 • to the study region (black triangle).Note that this is a combined plot for all stations using an average station position of 30.1 • N, 16 • W.

Figure 4 .
Figure 4. Averaged SKS-splitting measurements in the Madeira Archipelago.(a) Overview map.Circles represent individual stations with coloured azimuthal sectors showing splitting results in that azimuthal range.(b) The results are averaged over multiple stations for Madeira and the Desertas Islands, and (c) averaged over only one single station on Porto Santo.The colours of the sectors represent difference between FPD and present-day plate motion direction, taken from GSRM v.2.1 (Kreemer et al. 2014), also shown by the black arrow.The grey arrow indicates plate motion 60 Ma ago (Müller et al. 2019).(b) and (c) Averaged results of the different sector.Colours are the same as on the map.The sector boundaries are shown by grey lines.If stations produce significantly different patterns, multiple solutions are shown in a single sector.Arrows show plate motion.Details can be found in TableS1(Supporting Information).The time lag between fast and slow shear waves (δt) is represented by the length of the bar and the FPD (ϕ) is represented by the orientation of the bar.The number of individual results used to form the average is indicated by the thickness of the bars.For individual and averaged SKS splitting measurements including null measurements of these stations, see Figs S3 and S4 (Supporting Information).

Figure 5 .
Figure 5. Averaged SKS in the Canary Archipelago.()a Overview map.Islands that contain more than one station are indicated by a white circle, islands that only contain one station with results (La Gomera, Lanzarote) are indicated by a white triangle.Detailed maps of the other stations are shown in (b)-(f).(c) The total averaged result of station TBT (indicated by rectangular black outlines) on La Palma is shown to compare results from (i) Barruol & Ismail (2001), (ii) Behn et al. (2004) and (iii) this study (inset).

Figure 6 .
Figure 6.Individual splitting measurements for (a) DOC-16, Madeira, (b) PMPST, Porto Santo and (c) TBT, La Palma.Null results are shown by grey crosses.Results of other stations can be found in Fig. S4 (Supporting Information).

Figure 7 .
Figure 7. (a) Local S-splitting measurements in the Canary Islands, (b) also showing a detailed map of El Hierro, (c) as well as Madeira.Stations are shown by triangles; red triangles depict stations with S-splitting measurements.The time lag between fast and slow shear waves (δt) is represented by the length of the bar and the FPD (ϕ) is represented by the orientation of the bar.The colour of the bar indicates the event depth.All results are plotted on the midpoint between the event and station location, which due to the curvature of the ray path will be further away from the station than the midpoint of the path.For local waves this effect can be neglected.Rose diagrams show the summation of all FPDs for individual islands.Due to the abundance of results the map is separated into three major clusters on El Hierro (1-west; 2-south and 3-northeast).The rose diagrams are normalized to an angular bin size of 30 • .Note that their respective sizes are not linearly representative of overall number of events; all smaller rose plots have been enlarged to ensure visibility.

Figure 8 .
Figure 8. Local shear wave splitting delay times against event depth and in comparison with averaged teleseismic SKS splitting delay times.The shapes of the symbols correspond to the different islands on which the stations are located.For Canary Islands local S results, darker colours represent islands that are located to the east.For SKS results, the measurements are sorted by island from west to east with colours representing number of averaged events.Note that the depth information for the SKS-waves is arbitrary to facilitate visibility.Figures with results of further individual islands, see Fig. S6 (Supporting Information).

Figure 9 .
Figure 9. Apparent splitting parameters with azimuthal variation based on two-layer models for (a) and (b) Madeira and (c) and (d) Lanzarote, (a) and (c) showing fast shear wave orientation and (b) and (d) delay time between fast and slow shear wave.Horizontal bars indicate the weighted average of the individual sectors shown in Figs 4 and 5, coloured by the difference to the APM direction.The thickness of the bar indicates the number of individual results included in the sector.Individual results are shown by circles, using the same colours as the bars to which they belong.The curves show the analytic results at different periods.Here, the fast shear wave orientation of the lower layer is based on the current plate motion, whereas the upper layer is based on the general orientation of faults and rift zones (e.g.Carracedo 1999).The delay times are chosen to match the average delay times of the individual islands, assuming a dominating mantle anisotropic layer.Further islands are shown in Fig.S7(Supporting Information).

Figure 10 .
Figure 10.SKS shear wave splitting results from this study (red) in context of results from other studies (in white) in Morocco and Iberia, taken from the updated Wüstefeld et al. (2009) database (see references therein).The blue lines show the predicted ISA predictions derived from a global mantle flow model at 225 km depth (Conrad & Behn, 2010).

Table 1 :
Details of all averaged SKS splitting measurements in both regions (italic values are used to distinguish results from previous studies).
b DOC stations indicated by their numbers only.c Number of averaged events.d The azimuthal boundaries for CTAC are: −45 • -3 • ; 3 • -55 • .e

Table 2 :
Average delay times and percentage of anisotropy, islands ordered from west to east.Number of stations with good SKS/S spitting (number of stations with data but a lack of good results in brackets). a