On the Origin of Seismic Anisotropy in the Shallow Crust of the Northern Volcanic Zone, Iceland

The Icelandic crust is a product of its unique tectonic setting, where the interaction of an ascending mantle plume and the Mid‐Atlantic Ridge has caused elevated mantle melting, with the melt accreted and cooled in the crust to form an oceanic plateau. We investigate the strength and orientation of seismic anisotropy in the upper crust of the Northern Volcanic Zone using local earthquake shear‐wave splitting, with a view to understanding how the contemporary stress field may influence sub‐wavelength structure and processes. This is achieved using a data set comprising > $ > $ 50,000 earthquakes located in the top 10 km of the crust, recorded by up to 70 stations over a 9 year period. We find that anisotropy is largely confined to the top 3–4 km of the crust, with an average delay time of 0.10 ± 0.05 s, and an average orientation of the fast axis of anisotropy of N014°E ± 27°, which is perpendicular to the spreading direction of the Eurasian and North American plates (N106°E). These results are consistent with the presence of rift‐parallel cracks that gradually close with depth, the preferential opening of which is controlled by the regional stress field. Lateral variations in the strength of shear wave anisotropy (SWA) reveal that regions with the highest concentrations of earthquakes have the highest SWA values (∼10%), which reflects the presence of significant brittle deformation. Disruption of the orientation of the fast axis of anisotropy around Askja volcano can be related to local stress changes caused by underlying magmatic processes.

. (a) Overview of Iceland with the major glaciers outlined. The orange bands delineate the en-échelon fissure swarms that characterize the on-land expression of the northern Mid-Atlantic Ridge. The study region shown in panel (b) is outlined in red. (b) Shaded digital elevation map for the region around Askja volcano. Red triangles are seismic stations operated by the University of Cambridge used in this study. The purple square is the Icelandic Meteorological Office station, MKO. The entire earthquake catalog of Greenfield et al. (2018) is shown as gray dots, with those colored by hypocentral depth representing earthquakes used in this study. Two fissure swarms are highlighted: Askja's (purple) and Kverkfjöll's (green). The dashed line delineates the region associated with the Askja central volcano. The arrows show the regional direction of plate spreading, N106°E (c) An east-west section showing the earthquake catalog locations. (d) and (e) are circular histograms of surface features (fractures, fissures, faults) mapped by Hjartardóttir et al. (2016) for the Askja and Kverkfjöll fissure swarms, respectively, with average strikes of 18.4° and 23.9° shown by the black bars. n is the number of features in each sample.
Askja is a large, active central volcano located at the southern end of the NVZ (see Figure 1). A complex, nested sequence of at least three caldera-spanning 20 km-constitutes the main volcanic edifice, which is composed primarily of hyaloclastite and pillow lavas erupted during the last glacial maximum. The last eruption of Askja was in 1961, when a 2 km-long fissure opened up, with lava breaching the eastern side of the main caldera wall. Surface mapping around Askja has revealed a complex pattern of both caldera-concentric and rift-parallel features, including faults and surface fractures, which deviate in orientation from those observed elsewhere in the associated rift zone (Graettinger et al., 2019;Hjartardóttir et al., 2016).
Deformation around Askja has been monitored since the last eruption in 1961, at first with a tilt line within the caldera Tryggvason, 1989), but more recently using satellite-based GPS and InSAR measurements (de Zeeuw-van Dalfsen et al., 2012;Pagli et al., 2006;Sturkell et al., 2006). The long-term trend since 1961 is one of deflation, albeit at a decaying rate. Forward modeling based on geodetic observations suggest a shallow (3.5 km) Mogi-type source beneath the Askja caldera can explain the observed deflation, though most studies have assumed an isotropic, elastic half-space, which may be inappropriate around Askja (Drouin et al., 2017;Heimisson & Segall, 2020). For instance, rheological models based on a visco-elastic ridge appear to be key in the interpretation of geodetic data (Pedersen et al., 2009). Nonetheless, the derivation of such models from geodetic data provides insight into the contemporary stress state of the crust, wherein strain is gradually accumulated between rifting episodes. This complements the available seismic data set.
Large systems of fissures and faults are widespread across the rift segments associated with Askja and another central volcano, Kverkfjöll, situated to the southeast (see Figure 1). Broadly speaking, the fissure swarms associated with each central volcano are markedly seismically quiet in between rifting episodes. This part of the NVZ, however, stands out in that intense seismicity-which has been ongoing since records began in 1974 (Einarsson, 1991)-is observed between the fissure swarms of the Askja and the Kverkfjöll volcanic systems, extending northeast toward Herðubreið (a tuya formed during the last glacial period- Figure 1b). This seismicity has been explored in detail (e.g., Green et al., 2014;Greenfield et al., 2018;Soosalu et al., 2010), most recently by Winder (2021), who suggests that the anomalous levels of seismic activity may be a result of the combined effects of fluids migrating from depth, groundwater flow away from the hydrothermally active regions around Askja, and, consequently, the deposition of unusually high concentrations of clays, weakening the faults. These faults form a network of conjugate strike-slip faults (see Figure S1 in Supporting Information S1) that are bisected by the strike of the rift, suggesting some relation to plate spreading as well as a degree of interplay between the faults and pre-existing structures that are controlled by the stress field (Green et al., 2014;Winder, 2021). This seismicity tends to occur in swarms (where the earthquakes are clustered in both space and time), located primarily above the well-mapped brittle-ductile transition at around 6-8 km depth (Soosalu et al., 2010). There is also significant seismic activity in the Öskjuvatn caldera, which lies within Askja, associated with the migration of geothermally heated fluids, as well as a number of deep clusters of earthquakes thought to be associated with the migration of melt within the trans-crustal melt storage system .
Seismic anisotropy, the directional dependence of seismic wave-speed, has been observed in the crust across a range of environments (Boness & Zoback, 2006;Gao et al., 2019;Illsley-Kemp et al., 2018;Johnson et al., 2011;Kaviris et al., 2020). The nature of anisotropy can be broadly classified as either effective (i.e., a long-wavelength, bulk property of an otherwise heterogeneous medium) or intrinsic, arising from the anisotropic elastic structure at the crystal lattice level. Effective anisotropy is typically invoked as the primary mechanism by which seismic anisotropy is generated in the shallow, brittle crust. Here, mechanisms are typically related to either stress, through preferential closure of micro-cracks, known as Extensive Dilatancy Anisotropy (EDA; Crampin, 1994;Zatsepin & Crampin, 1997), oriented melt pockets (OMP; Bastow et al., 2010;Holtzman et al., 2003;Keir et al., 2005), structure, such as repeating isotropic layers (Backus, 1962) or structurally controlled anisotropy, related to damage zones around faults (Boness & Zoback, 2006;Li & Peng, 2017;Pastori et al., 2019;Zinke & Zoback, 2000). We seek here to determine the mechanism-or mechanisms-responsible for generating seismic anisotropy in the crust around Askja in order to understand better the present-day state and anisotropic structure of the nascent crust formed at a mid-ocean ridge, as well as how large-scale stress fields (e.g., due to plate spreading) can be disrupted by local perturbations, such as the inflation or deflation of shallow magma storage systems or magmatic intrusions. Spatially mapping the seismic anisotropy in the crust around Askja, and other volcanoes, provides important context in which to study how the crust in volcanic environments responds to such transient perturbations to the local stress field, which-if well understood-may lead to real-time insights into the state of active volcanic systems (Johnson et al., 2011).
Shear-wave splitting is one of the most unambiguous indicators of seismic anisotropy. When a linearly polarized shear wave impinges on an anisotropic medium, it is partitioned into two quasi-S waves, which propagate with different wavespeeds. The polarizations of these two waves, commonly referred to as the 'fast' (hereafter denoted ϕ) and 'slow' axes, are controlled by the symmetry and orientation of the anisotropic elastic tensor. A time lag, δt, accrues between the polarized waves as they propagate through the region, with the final integrated value proportional to both the path length and the strength of anisotropy. Significant work has been done to establish methods that can distinguish between structural and stress-induced anisotropy (Boness & Zoback, 2006;Johnson et al., 2011;Li & Peng, 2017;Zinke & Zoback, 2000). Being able to distinguish these causes is critical for the potential application of time-series analysis to shear-wave splitting observations as a means of monitoring the evolution of the stress field in volcanic environments in response to seasonal signals, long-term temporal signals (such as deflation and inflation), and stress transients resulting from volcanic processes such as caldera collapses and dike intrusions (e.g., Johnson & Poland, 2013).
Here, we perform local earthquake shear-wave splitting analysis in the neighborhood of Askja volcano, in order to relate observed anisotropy to the underlying processes responsible for the accretion of new crust at a mid-ocean ridge and the development of associated volcanic systems. The results provide a new perspective on a region that is already well studied using complementary geophysical methods (de Zeeuw-van Dalfsen et al., 2012;Drouin et al., 2017;Greenfield et al., 2016;Sturkell et al., 2006).

Data
We use continuous seismic data recorded by a network of 3-component seismometers operated by the University of Cambridge since 2008, with additional data from one instrument operated by the Icelandic Meteorological Office (MKO, denoted by the purple square in Figure 1). Over time, the network has consisted of between 30 and 70 broadband instruments, primarily Güralp 6TDs (30 s corner frequency). All data used in this study were recorded using Güralp 6TDs. For the shear-wave splitting analysis, we use the earthquake catalog of Greenfield et al. (2018) which spans 2009-2015, updated (using the same methodology outlined in their paper) to include data recorded between 2015 and 2018 . These earthquakes were detected and located using the automatic Coalescence Microseismic Mapping algorithm (CMM; Drew et al., 2013). The details of pre-processing applied to the data to generate this catalog is available in Greenfield et al. (2018). The CMM algorithm produces automatic arrival time picks for P-and S-phases that were used, along with some manually picked phase arrivals, to relocate the events using the probabilistic, non-linear earthquake location method implemented in NonLinLoc (Lomax et al., 2000). The final catalog consists of 58,143 individual earthquakes spanning a local magnitude range of −0.6-4.0, with a magnitude of completeness of ∼1.
The majority of earthquakes (52,141, or 89.7%) occur in the brittle, upper 7 km of the crust, generated primarily by a network of cross-cutting conjugate strike-slip faults oriented approximately N-S and NE-SW, located to the northeast of Askja volcano and to the south of Herðubreið. The remaining shallow seismicity is related to geothermal processes at Askja volcano. The depths of these shallow events are well-distributed throughout the brittle crust. The final 10.3% (6,002) of events in the catalog occur in pockets at depths between 7 and 25 km in the typically aseismic lower, ductile crust and are thought to be associated with magmatic processes Martens & White, 2013). We limit our analyses to splitting observed from earthquakes originating in the upper 10 km of crust in order to focus specifically on anisotropy in the shallow crust. Finally, we exclude any events that occurred between August 2014 and January 2015 in order to remove the possible effect of stress transients related to the 2014-15 eruption of Bárðarbunga. Figure 1 illustrates the spatial distribution of earthquakes and seismic stations between 2009 and 2018 that have been used in this study.

Shear-Wave Splitting
We measure the shear-wave splitting parameters (ϕ, δt) using the Multiple Filter Automatic Splitting Technique package (MFAST version 2.2; Savage et al., 2010;Teanby et al., 2004), which uses the eigenvalue minimization algorithm of Silver and Chan (1991). Figure 2 illustrates the output from MFAST for a good quality event; further examples can be found in the supplementary information (see Figures S2-S5 in Supporting Information S1). Unlike other methods, this does not require any knowledge of the initial polarization, which is often difficult to assess a priori in local earthquake datasets, though at the cost of being more prone to cycle skipping. A grid search over ϕ and δt is used to find the pair of values that best remove the observed splitting, determined by measuring the linearity of particle motion on the horizontal components within a window around the S-phase arrival. This is further automated by trialing multiple windows around the S-phase arrival and applying cluster analysis to the ensemble of results in order to identify stable solutions. This is a particularly effective means of identifying cycle skipping. Errors for individual measurements are calculated by conducting an F-test and finding the 95% confidence interval on the optimal (ϕ, δt) (Walsh et al., 2013). Each measurement is automatically graded based on the distribution of clusters and the tightness of the misfit contours from the grid search (Savage et al., 2010). MFAST also trials a suite of filters over the S-phase pick in order to determine the filter that most effectively boosts the signal-to-noise ratio. Table S1 in the Supporting Information S1 provides an overview of the final suite of MFAST parameters used in the analysis stage.
We limit our analyses to the subset of measurements that satisfy the following criteria: a signal-to-noise ratio (as defined in Savage et al., 2010) greater than 4; clusters graded "ACl" (a measure of the number of clusters identified and how tight they are); errors in ϕ < 10° in order to mitigate erroneous observations resulting from cycle skipping; values of δt < 0.48 s, equal to 0.8 times the maximum delay time of the search; and errors in δt < 0.05 s as an additional filter against 'null' measurements and poorly constrained results. A null measurement can occur when there is no anisotropy in the plane of the shear wave particle motion, or when the source polarization of the shear wave is along the fast or slow axes of the anisotropic medium. Source polarizations are determined from the corrected horizontal particle motion. Measurements of ϕ within 20° of the source polarization are considered too ambiguous in that they cannot be definitively distinguished from nulls, and are subsequently excluded from Figure 2. An example of a good splitting measurement. (a) shows the raw data for the East (green), North (orange), and Vertical (blue) components. (b) shows a zoom in around the S phase arrival rotated onto the nominal 'radial' (P) and 'transverse' (T) axes before and after correction for splitting. (c) and (d) show the phase arrivals rotated onto the 'fast' and 'slow' axes before and after correction, with (e) and (f) showing the corresponding particle motion. There is a clear linearization of the particle motion of the horizontal components and removal of energy from the transverse component. Panels (g)-(i) show the results of the multiple window trials and the cluster analysis. Finally, (j) shows the resultant grid of the minimized eigenvalue. The blue cross denotes the optimal (ϕ, δt) pair. further analysis. After applying these criteria, we are left with 69,119 measurements of shear-wave splitting. We further remove measurements for which the angle of incidence of the shear wave at the surface falls outside the shear-wave window (Nuttli, 1961). This window, defined by sin −1 (V s /V p ), is the angle to the vertical at which there will be non-negligible interactions with the free surface that would alter the phase and amplitude information on the horizontal components (Crampin, 1984). A V p /V s ratio of 1.78 corresponds to a shear-wave window of ∼34° from the vertical. However, volcanic environments typically exhibit very low velocities in the topmost layers (Lesage et al., 2018), which will cause significant deflection of the ray toward the vertical. Therefore, we limit our analysis to event-station pairs with a straight-line angle-of-incidence at the surface of <50°. Finally, all remaining splitting measurements were visually inspected to filter out any poor results, with over 90% of the measurements being retained; the primary cause of a poor result was either cycle skipping or a poor initial S-phase pick. The results of this manual filtering stage are presented in Figure S6 of Supporting Information S1. This leaves 9,974 high-quality measurements of shear-wave splitting.

Shear Wave Anisotropy
The delay time, δt, is an integrated measure of the strength of anisotropy along the raypath, making it an unsuitable metric for direct comparison between different event-station pairs. This is overcome by converting the observed delay times to shear wave anisotropy (SWA: Thomas & Kendall, 2002), which is a measure of the strength of anisotropy as defined by the fractional perturbation, a, from the average shear wave speed, : where t slow and t fast are the slow and fast traveltimes, respectively. SWA is also a function of both the path length, d, and the velocity structure along the ray, therefore representing a more appropriate metric to compare between individual observations. We assume straight-line raypaths and use an optimal 1-D velocity model determined by inverting microseismic arrival times (Mitchell et al., 2013). Nowacki et al. (2018) demonstrated that the errors introduced by the straight-line raypath assumption are negligible for shallow events, for which the raypaths do not deviate far from a straight line, with up to around 1% overestimation in SWA for the deepest events. Additionally, they show that the uncertainty in SWA arising from inaccuracy in the velocity model is estimated to be less than 1% from bootstrap modeling. The geological setting of Afar and the corresponding seismic data set are markedly similar to that of Iceland in a number of ways. While the erupted products at Aluto, in Afar, are more silicic in composition, the velocity structures in both regions exhibit low velocities in the shallow crust that cause rays to turn sharply toward the vertical. Furthermore, the seismicity used in Aluto spans a similar depth range as the data set presented here, with the majority of events occurring above ∼10 km depth. Given these similarities, we think it likely that this uncertainty analysis remains appropriate for the Iceland data set.

Delay Times
From the entire catalog of shear-wave splitting measurements, we recover an average delay time of δt = 0.10 ± 0.05 s. This value is consistent with similar datasets elsewhere, e.g. ∼0.2 s around Soufrière Hills volcano, Montserrat (Baird et al., 2015); 0.1-0.2 s in the Western Volcanic Zone, Iceland (Menke et al., 1994); and 0.11 ± 0.06 s around Aluto volcano, Ethiopia (Nowacki et al., 2018). We find the distribution of delay time observations to be sufficiently normal to justify the extraction of a regional 1-D depth profile as the central tendency of the data via the application of a rolling arithmetic mean (Figure 3). We use a 1.5 km rolling window, spaced every 0.75 km, which is chosen to reflect the uncertainty in hypocentral depth for shallow events . We observe a constant delay time at depths 3 km. Between 3 km depth and the surface, there is a suggestion that the delay time begins to trend toward 0, which is consistent with a finite-thickness anisotropic layer in the very shallow crust, a common observation across volcanic environments (Johnson et al., 2011;Menke et al., 1994;Nowacki et al., 2018). This does not preclude structural control on anisotropy, but it is a key requirement for stress-induced anisotropy due to the preferential closure of micro-cracks. In oceanic-type crust, most pore space has been closed by lithostatic pressure at around 4-5 km below the surface (Christensen, 1984). The relationship between crustal porosity and depth can be expressed as the exponential function (e.g., Athy, 1930;Audet & McConnell, 1992): where c is a constant (∼6.15), Φ1 is the surface porosity, P(r) is the lithostatic overburden pressure (=ρ(r)gd, where ρ is the density, g is the acceleration due to gravity, and d is the depth), and P c is the characteristic closing pressure of the material (Han et al., 2014). We perform a simple fit of a similar exponential function to the depth profile, shown in Figure 3, which suggests that the 1-D behavior of the delay time is consistent with the presence of crustal cracks that gradually close with increasing depth.

Fast Axis Orientation
We observe an average orientation of ϕ∼N014°E ± 27° for the fast axis of anisotropy, though we recommend caution in drawing too much from the exact value of, and the uncertainty on, this measure as the circular statistics used are only appropriate if the observations are drawn from a unimodal distribution. Small variations in fast polarizations across the region, such as those expected in response to, for example, a rotation in the stress field, may be contributing to the large spread in observed ϕ values. The average orientation correlates well with the normal to the plate-spreading direction, as shown in Figure 4, as well as the mapped surface structures which likely reflect the long-term regional stress field (Hjartardóttir & Einarsson, 2012;Hjartardóttir et al., 2016). This is consistent with observations made at other spreading centers, such as the northern Main Ethiopian Rift . Exactly how the orientation of the fast axis of anisotropy varies across the region is investigated further in Section 3.2.2.

Shear Wave Anisotropy
Measurements of delay time are converted to SWA using Equation 1, as described in Section 2.3. We constrain the shallow anisotropic layer to be entirely above ∼3 km b.s.l. (i.e., a 4 km thick layer), inferred from the constant delay time below this depth observed in the 1-D profile (Figure 3). Assuming that the mechanism generating seismic anisotropy is aligned fractures in the shallow crust, this value is consistent with measures of fracture density from other independent measures, such as radial anisotropy constrained by surface waves extracted from ambient noise (Volk, 2021), response of velocity changes, dv/v, to seasonal changes in load (Donaldson et al., 2019), and general profiles of pore space as a function of depth in oceanic crust (e.g., Carlson & Herrick, 1990). While there is an element of bias in assigning the splitting observation to a single point in space, we follow precedent and use the mid-point of the raypath ( Figure 5) between the source and receiver before re-gridding the data. For near-vertical raypaths, as is the case for the majority of our data set due to the shear-wave window constraint, this introduces negligible systematic error in the pattern of lateral variations. The application of a symmetric 2-D Gaussian spatial filter to the re-gridded observations further reduces the impact of this assumption on the observed lateral patterns. Here, we present the results for a grid with 0.5 × 0.5 km 2 cells and a minimum observation count of 10, and 2-D Gaussian spatial filter with a half-width of 1 km ( Figure 6). The key features of the lateral variation in anisotropy strength are robust to perturbations to both the grid parameters and the smoothing radius. We trialed cell sizes varying from 0.25 × 0.25 km 2 to 1 × 1 km 2 , minimum number of observations per bin between 3 and 15, and a smoothing radius of 1-3 km, and found that the results did not vary significantly (see Figure S7 in Supporting Information S1). We acknowledge that the process of re-gridding the data in this way means that some azimuthal information is lost, but we deem it acceptable for the purpose of identifying trends in the strength of anisotropy across the rift segment. We measure an average anisotropic strength of ∼5%, with values ranging between 2% and 12%, which spans the appropriate range expected for mechanisms proposed for elastic anisotropy of the crust.

Fast Axis Orientation
We re-grid the observations of ϕ by grouping them laterally by the midpoint along the event-receiver raypath ( Figure 5), with the results presented in Figure 7. We use an adaptive quad-tree gridding method, which allows us to increase the detail (down to a minimum cell size of 2 × 2 km 2 ) where we have a higher density of observations. The minimum cell size used is the same order as the uncertainties in the epicentral locations for the earthquakes in the catalog. Starting from a single cell spanning the entire study region, this process recursively subdivides a cell into four sub-cells if the number of observations in the cell exceeds 200. Any cells with fewer than 50 observations are omitted from the final grid. Again, we trialed a number of values for these three parameters, the results of which are shown in Figure S8 of Supporting Information S1. Within each cell, the resultant vector is evaluated from which both the average orientation and the mean resultant length, ̂ , is determined. ̂ is a measure of dispersion analogous to the variance (in the opposite sense)-values close to 0 imply near uniform dispersion, whereas values close to 1 suggest that the orientations are tightly bunched around a particular orientation (e.g., Davis & Sampson, 1973). This allows us to observe the lateral trends in the orientation of anisotropy, without constraining the source of anisotropy to be in the vicinity of the source or the receiver.

Anisotropy Orientation and Strength
Our analysis of shear-wave splitting from earthquakes in the brittle, upper 10 km of crust around Askja has constrained the primary source of anisotropy to be in the top 3-4 km of crust, with the dominant orientation of the fast axis of anisotropy correlating strongly with the strike of the rift (Figure 4). Together, these two observations provide compelling evidence for extensional stress related to plate spreading as the underlying mechanism generating the observed seismic anisotropy. This is consistent with other studies of local shear-wave splitting in similar environments, such as the East African Rift (Illsley-Kemp et al., 2017;Keir et al., 2005;Nowacki et al., 2018), the volcanic zones of New Zealand (Illsley-Kemp et al., 2019;Johnson et al., 2011), and the Corinth Rift (Kaviris et al., 2017). The average delay time of shear-wave splitting observations (δt = 0.10 ± 0.05 s) is also consistent with these studies.
Although we attribute our observations of shear wave anisotropy to fractures or cracks in the shallow crust, there are other causes of anisotropy that may be a factor. For instance, aligned melt pockets could produce a signature of effective anisotropy with ridge-parallel orientation of the fast axis, as has been suggested in the upper mantle and lower crust of the Main Ethiopian Rift (Hammond et al., 2014;Kendall et al., 2005). However, there is no evidence for the presence of melt in large volumes in the shallow crust beneath the NVZ, outside of the central volcanic systems, and ambient noise studies that constrain azimuthal variations of radial anisotropy are not consistent with such a mechanism (Volk, 2021). Furthermore, it is difficult to propose a physically coherent reason as to why melt pockets would be focused in the very shallow crust, yet be absent at greater depth, which would need to be the case to explain the trend shown in Figure 3. Another possibility is Lattice Preferred Orientation (LPO) anisotropy associated with deformation, lava flows, or depositional processes. Recent measurements of radial anisotropy from ambient seismic noise  support the presence of LPO in the crust resulting from internal deformation or flow, but this appears to be largely restricted to depths below about 15 km, and therefore is unlikely to influence our results. Lava flows can align minerals such as plagioclase and clinopyroxene (Boiron et al., 2013), but this tends to occur at very short scale lengths horizontally and in depth, and consequently are also unlikely to substantially contribute to our pattern of anisotropy over what is a relatively large study area.
When interpreting the map of SWA (Figure 6), we recommend that a greater importance be placed on the relative values, as opposed to the absolute values, which can be 'tuned' by varying the thickness chosen for the anisotropic layer. We primarily see elevated values of SWA in regions with elevated rates of seismicity, which is consistent with the idea that stress is the primary control on the mechanism generating anisotropy. There is a region of elevated SWA to the south of Herðubreið (Figure 6), which corresponds with a region of elevated seismic activity on the network of cross-cutting conjugate strike-slip faults, which have been extensively studied (Green et al., 2014;Winder, 2021). Consequently, we can infer that this section of crust is heavily fractured and highly stressed, two conditions under which one would expect to see a higher anisotropic signal. This may, however, also be an artefact of the assumption that the anisotropic layer has a uniform thickness across the region. For example, elevated volumes of fluids within the crust may hold pore spaces open at greater depths. There is indeed some suggestion from the distribution of seismicity that the thickness of the brittle layer within the crust does exhibit some variability over the region (Soosalu et al., 2010), the impact of which could be explored in the future. The relatively low values of SWA to the northeast of Herðubreið correspond to a region of elevated V p /V s observed in a tomographic study of the region (Greenfield et al., 2016), which was interpreted to be a sign of elevated fluid content. This is consistent with the suggestion from Nowacki et al. (2018) that a higher V p /V s may indicate that there are more fluids present, which in turn causes lower effective anisotropy, and may also explain the relatively low SWA values around the Askja geothermal field, on the eastern edge of the Öskjuvatn caldera. However, we should note that elevated V p /V s need not necessarily imply lower anisotropy; for instance, Wang et al. (2012) made laboratory observations of cracked samples and carried out effective media modeling, which suggested that the presence of high V p /V s ratios may also be indicative of significant crack-induced anisotropy. Consequently, it may prove beneficial to explore models for the effective elastic stiffness of a medium hosting cracks with varying aspect ratios and fluids and compare these to the anisotropic response (if any) of the crust to stress transients, such as the 2014 Bárðarbunga-Holuhraun dike intrusion, which triggered seismicity on these faults (Winder, 2021).
The spatial trends in the orientation of the fast axis of anisotropy are broadly consistent with both the observed surface features from geological mapping and the plate-spreading direction (Figure 4). This is consistent with findings from other rift environments (Illsley-Kemp et al., 2017;Menke et al., 1994;Nowacki et al., 2018), where the fast axis of anisotropy was found to be aligned to the present-day minimum compressive stress i.e. rift parallel. In these studies, the source of anisotropy is also attributed to aligned cracks in the top 3-4 km of the crust. Such crack alignment in the very shallow crust is also present in other tectonic environments, including fold and thrust belts. For example, de Lorenzo and Trabace (2011) investigate local earthquake shear-wave splitting using data recorded in the central Appenines, and attribute anisotropy in the top 4-5 km of the crust to fault-parallel fluid-filled crack systems.
As Figure 7 illustrates, the orientations of the fast axis of anisotropy are not uniformly rift parallel; for instance, in the very south they have a stronger easterly component compared to those in the north. It is likely that the regional stress field in the southern region is overprinted by the ongoing deformation that is taking place around Askja, as well as potentially being affected by the loading of the crust due to the presence of Vatnajökull to the south. Subsidence of the main caldera has been ongoing since 1983(de Zeeuw-van Dalfsen et al., 2012, possibly due to the cooling and contraction of an underlying magma body, although recent micro-gravity increases may be due to magma flow into a shallow magma chamber (de Zeeuw-van Dalfsen et al., 2013). Such local stress changes and associated deformation may be responsible for scattered horizontal velocity vectors measured by GPS stations in the vicinity of Askja (Árnadóttir et al., 2009;Drouin et al., 2017); consequently, the disruption to the pattern of anisotropy around Askja is perhaps not surprising. In the presence of a deflating source, the orientation of the maximum horizontal stresses within the crust tends to be concentric to the center of deflation (e.g., Johnson & Poland, 2013). In contrast, the maximum horizontal stresses induced by an inflating source tend to be radial to the center of inflation. In the next section, we use stress modeling to investigate the relationship between the stresses associated with the deflation of Askja and plate spreading, and the measured orientations of the fast axis of anisotropy across the region.

Stress Modeling
Numerous studies have concluded that the orientation of anisotropy in the crust is generally controlled by the regional stress field and/or the alignment of structures, such as fissures and faults (Illsley-Kemp et al., 2018;Johnson et al., 2011;Savage et al., 2010). Distinguishing between stress-induced and structural anisotropy in the Northern Volcanic Zone is made somewhat more complex by the fact that the regional stress field is also the primary control on the orientation of structural features. It is observed, however, that the system of faults between the Askja and Kverkfjöll rift segments (responsible for a large proportion of the tectonic seismicity in the region) is composed of conjugate strike-slip faults oblique to the strike of the plate margin (see Figure 1). This suggests that we can rule out fabric resulting from the damage zones around faults as a mechanism generating (significant) Figure 7. Lateral variations in observed fast axis orientations, ϕ. The observations have been assigned to the midpoints between source and receiver, then re-gridded using a quadtree method. The resultant grid is plotted using faint black lines. Within each cell, the bar represents the average fast orientation, colored by the 'resultant vector' which is a measure of dispersion/coherence of the orientation data. Darker colors indicate stronger coherence. The same map, but with circular histograms instead of the circular mean within each cell, is shown in Figure S9 of Supporting Information S1. anisotropy, based on the regional averages. This observation may have implications as to the nature of faults that can influence the surrounding anisotropic field and the spatial extent over which they are capable of doing so. The strike-slip faults in question do not exhibit typical mainshock-aftershock behavior, instead relieving strain accumulated in the brittle crust via swarms of earthquakes that migrate along the fault surface. Since records began, no earthquakes with a local magnitude greater than 4 have been observed on these faults, despite a number of them being sufficiently long to do so.
We explored the role of stress in the generation of anisotropy by modeling the regional stress field around Askja using the Coulomb v3.3 software package (Toda et al., 2011). Whereas the ductile lower crust (deeper than around 6-8 km, based on the depth extent of the seismic catalog (Soosalu et al., 2010)) is able to deform by continuous creep under the extensional stresses, accretion and extension of the brittle upper crust is episodic in nature. Over time, elastic strain accumulates in the brittle crust, before being released over short, intense periods of diking and extensional faulting, as seen during the 2014-15 eruption of Bárðarbunga (e.g., Ágústsdóttir et al., 2019;Sigmundsson et al., 2014;Woods et al., 2019). We model this process using a buried dislocation, which has previously been used to model plate boundary deformation in the rift zones of Iceland LaFemina et al., 2005). This model assumes that spreading below the brittle-ductile boundary is constant and equal to the full-spreading rate, represented by an opening Okada dislocation (Okada, 1992) extending from the locking depth to infinite depth. The stress singularity at the upper edge of the buried dislocation is eliminated from the model by tapering the dislocation such that the opening gradient goes to 0 at the topmost edge (Heimisson & Segall, 2020). The spreading boundary is taken to pass through Askja, striking along the rift segment at N015°E. A small component of spreading is assigned to the Kverkfjöll rift segment, though it is debatable whether any active spreading is occuring in this region (Drouin et al., 2017). However, this inclusion does not significantly impact the modeled stress field. The ongoing deflation beneath Askja is incorporated using the best-fitting (analytical) solution from forward modeling of GPS data (Drouin et al., 2017). This results in a point Mogi source at 3.5 km depth beneath the Askja caldera (see Figure 8), with a volumetric change of 0.0013 km 3 /year. While both models are highly simplified, neglecting visco-elasticity in particular, they are sufficient to capture, to first order, the tectonic stress state of the crust. The input files for this modeling are available in the Supporting Information S1.
Using the method of Lund and Townend (2007), we extract the maximum horizontal stress vectors (S Hmax ) from the final model at a depth of 0 km b.s.l., where we expect the impact of the stress field to have the most significant effect on the opening/closure of cracks. We observe a strong correlation between the orientations of fast directions and S Hmax across the region, including a similar rotation moving from south to north. This provides a strong link between the stress field and the anisotropy, as would be expected for the EDA mechanism. The differences, particularly at the southern end of the region, are potentially due to the component of strain imparted by the presence of the Vatnajökull ice cap, which is not included in the modeling. Interpolating the strain field directly from the available GPS data may prove valuable in assessing how much of the observed rotation is due to the unmodeled components. Around Askja, the modeled strain field shows a similar level of scatter to what is observed in Figure 7, though there is no particular coherency in alignment. This is likely to be due to the limited spatial resolution of the splitting measurements, coupled with the simplifying assumptions made in the stress modeling. Careful analysis of the temporal changes in the anisotropic signal in response to stress transients, such as the 2014 Bárðarbunga-Holuhraun dike intrusion, may provide more supporting evidence for the EDA mechanism dominating the generation of anisotropy in the upper crust in the Northern Volcanic Zone.

Conclusions
We have presented shear-wave splitting results from the Northern Volcanic Zone, Iceland, based on a large data set of local earthquakes that span a period of 9 years. The dense, stable network has allowed us to image the anisotropic properties of the Icelandic crust with a high spatial resolution. These observations have allowed us to investigate the likely mechanisms generating this anisotropy, whether controlled by the stress state, or structural features in the crust. The main findings of this study include (a) based on earthquakes that occur between the surface and 10 km depth, anisotropy is largely restricted to the top 3-4 km of the crust; (b) delay time variations in the shallow anisotropic layer are consistent with the presence of cracks that gradually close with depth; (c) SWA is strongest in regions of elevated seismicity, particularly in the zone between the Askja and Kverkfjöll rift segments, which appears to be heavily fractured; (d) the dominant orientation of the fast axis of anisotropy is almost perpendicular to the spreading direction, which indicates that regional stress is the dominant control on anisotropy; and (e) in the neighborhood of Askja, the orientation of the fast axis of anisotropy becomes scattered, which is consistent with stress modeling results that use a Mogi source located 3.5 km beneath the main caldera. Future work will focus on the very deep earthquakes beneath the Northern Volcanic Zone, and the constraints they may be able to supply on anisotropy in the lower crust, which has previously been imaged by ambient noise tomography.

Data Availability Statement
All waveform data, except that recorded at the IMO station MKO, are available for download from the IRIS DMC, under the following network codes: 4F (2007-2011, https://doi.org/10.7914/SN/4F_2007), Z7 (2010-2015, https://doi.org/10.7914/SN/Z7_2010), and 8K (2016-2022, https://doi.org/10.7914/SN/8K_2016, which is expected to become fully available by 2023). Waveform data from MKO are only available under request to the Icelandic Meteorological Office (https://en.vedur.is/about-imo/contact/). The shear-wave splitting measurements and Coulomb stress modeling input file can be downloaded from https://doi.org/10.5281/zenodo.5007022. Data analysis was carried out using Python (3.9), with the following packages proving particularly useful: Ob-sPy (1.2.2, Beyreuther et al., 2010); SciPy (1.7.1, Virtanen et al., 2020); NumPy (1.21.3, Harris et al., 2020); Pandas (1.3.4, pandas development team, 2021); pyproj (3.2.1, Snow et al., 2021). Data visualizations were performed using Matplotlib (3.4.3, Hunter, 2007) and Generic Mapping Tools (6.2, Wessel et al., 2019). All code required to reproduce the analysis and visualizations presented in this study can be downloaded from https://doi.org/10.5281/zenodo.5636924.  (loans 857, 914, 968, 980, 1,022 and 1,115). The authors thank the many people who participated in seismic data collection fieldwork in Iceland over the years, in particular Bryndís Brandsdóttir and Sveinbjörn Steinþórsson, without whom this work would not be possible. The authors also thank the Iceland Meteorological Office (IMO) for kindly providing the additional data from MKO. The authors especially thank Tim Greenfield and Tom Winder for their considerable efforts in producing the earthquake catalog that made this study possible and helpful discussions along the way. Department of Earth Sciences, Cambridge contribution ESC.6035.