Evidence of a Shallow Magma Reservoir Beneath Askja Caldera, Iceland, From Body Wave Tomography

In August 2021, Askja caldera switched to reinflation following ∼ 40 years of continuous deflation that was first measured some 20 years after its last eruption in 1961. Various lines of evidence, including from geodetic modeling, suggest that both the deflation and reinflation events are related to a shallow magma body. To better understand the subsurface plumbing system, we derive P‐wave velocity (Vp), S‐wave velocity (Vs), and Vp/Vs models of the mid‐upper crust by leveraging a new local earthquake traveltime data set. A cylindrical low‐velocity zone, ∼ 3 km wide and extending to ∼ 8 km below sea level (bsl), was imaged beneath the caldera. Within it, two distinct lower velocity and higher Vp/Vs anomalies are illuminated, one centered at ∼ 0.5 km and the other at ∼ 6 km bsl. The shallower anomaly lies directly beneath the zone of uplift and is likely associated with the current reinflation event.


Introduction
Askja caldera is situated in the Northern Volcanic Zone (NVZ) of Iceland, as shown in Figure 1.After its most recent eruption in 1961, leveling measurements that began in 1966 within the Askja caldera indicate that the region was inflating during the period 1966-1973, after which a continuous signal of deflation was recorded between 1983 and 2021, albeit at a gradually declining rate (Sturkell, Sigmundsson, & Slunga, 2006).This observation has been explained in terms of the drainage, cooling, and contraction of a shallow magma body located ∼2 km bsl (de Zeeuw-van Dalfsen et al., 2012).In August 2021, a switch from deflation to inflation was detected within the Askja caldera through Interferometric Synthetic Aperture Radar (InSAR) measurements; this was corroborated by GPS measurements in the caldera, with peak surface uplift of over 35 cm in the first year of reinflation at a GPS station near the western edge of lake Öskjuvatn.
The existence of a shallow magma body beneath Askja caldera has been supported by studies from multiple disciplines, including but not limited to geodetic modeling (de Zeeuw-van Dalfsen et al., 2012); ambient noise tomography (Volk, 2021); observations of higher attenuation associated with seismic rays passing directly beneath the caldera compared to those that do not (Greenfield et al., 2016); local earthquake shear wave splitting analysis revealing scattered fast axis orientations of anisotropy beneath the caldera (Bacon et al., 2022); and a petrological study (Sigurdsson & Sparks, 1981).In the case of the ambient noise study, both horizontal and depth resolution were limited by the period range used, resulting in a somewhat nebulous low-velocity anomaly that has no clear spatial association with the current uplift.Shear wave velocity (Vs) is more sensitive to the presence of melt compared to primary wave velocity (Vp), although Vp/Vs is often regarded as the most definitive indicator (Hammond & Humphreys, 2000).Previous studies have used local earthquake tomography to image seismic velocity structure beneath Askja; for instance, a lower velocity and higher Vp/Vs region was found ∼6 km bsl, and interpreted as a magma storage region (Greenfield et al., 2016;Mitchell et al., 2013).However, imaging of the velocity structure at shallower depths was limited due to a lack of path coverage immediately beneath the caldera, which may explain why the shallow melt anomaly inferred by other studies has not been imaged by this technique.Fortunately, new seismic recordings between 2016 and 2022 that better sample the shallow crust beneath Askja caldera are now available, thus providing a timely opportunity to reappraise seismic properties of the volcanic plumbing system through local earthquake tomography.Here, we present the results of such a study, which assumes new significance in light of the recent reinflation episode.

Geology and Tectonics
Iceland is located on the divergent boundary associated with the mid-Atlantic Ridge, where the North American and Eurasia plates separate, and above an active upwelling mantle plume (Celli et al., 2021;French & Romanowicz, 2015;Morgan, 1971;Shen et al., 2002;Wolfe et al., 2002).This distinctive geological setting gives rise to both active tectonics and volcanism, predominantly concentrated along several rift zones.The abundant melt production has resulted in the formation of a remarkably thick oceanic-type crust, reaching up to ∼40 km thickness (Darbyshire et al., 1998), and the emergence of the spreading ridge above sea level.The Western Volcanic Zone (WVZ), Eastern Volcanic Zone (EVZ), and Northern Volcanic Zone (NVZ) (see Figure 1 inset) are regions of active extension and volcanism that accommodate the eastern ridge jump caused by the underlying plume (Einarsson, 1991).
Askja caldera is located within the NVZ and is characterized by a set of at least three nested calderas.A significant explosive eruption occurred in 1875, which produced the innermost caldera that accommodates the 11 km 2 Lake Öskjuvatn (Hartley & Thordarson, 2012).In the 20th century, eruptions occurred from vents close to Lake  (Einarsson & Saemundsson, 1987) of Iceland showing the location of fissure swarms (in lighter gray), the direction of plate spreading, and the location of the study region (orange box).The three volcanic zones are the Western Volcanic Zone (WVZ), the Eastern Volcanic Zone (EVZ), and the Northern Volcanic Zone (NVZ).The green lines denote transform zones that accommodate relative movement between the mid-Atlantic Ridge and the rift zones in Iceland.
Öskjuvatn (Trippanera et al., 2017).It has been proposed that such eruptions are triggered by the upward movement of tholeiitic basalt magma from deeper sources, which then replenish a shallower magma reservoir (Sigurdsson & Sparks, 1981).
A recent seismic anisotropy study that exploited local earthquakes (Bacon et al., 2022) detected the presence of shallow fracturing in the caldera and its surrounding regions, which is also indicated by the dense microseismicity (Greenfield et al., 2020;Winder, 2021).Furthermore, the orientation of the fast axis of anisotropy inside the caldera appears to be scattered in comparison to those outside the caldera, which tend to be rift parallel.Stress modeling (Bacon et al., 2022) supports this disruption being caused by a shallow contracting magma body.

Data and Method
We utilized arrival time data from 3,348 local earthquakes that occurred between 2006 and 2022, recorded by up to 92 seismic stations located in the neighborhood of Askja caldera (see Figure 1).The network consists of 3component broadband sensors, which are primarily Guralp 6TDs with a corner period of 30 s.A smaller number of Guralp 3ESPCDs (60 s) and 3Ts (120 s) are also used (Greenfield et al., 2020).The earthquake locations were determined using Nonlinloc (Lomax et al., 2000) with manually picked P-and S-wave arrival times that have picking uncertainties of less than 0.10 s for P-waves and 0.15 s for S-waves.The complete data set includes 53,474 P-arrival times and 47,231 S-arrival times.The initial data set, comprising 23,115 P-arrival times and 19,060 S-arrival times, was sourced from a previous study (Greenfield et al., 2016), and additional picks from recently collected data were incorporated.The majority of the new picks are distributed beneath the caldera in the upper crust, thereby providing additional path coverage to image the shallow plumbing system.For consistency, the new picks were made using a similar procedure to the picks from Greenfield et al. (2016), in which P-wave and S-wave arrival times determined by the detection and location software QuakeMigrate (Winder et al., 2020) are manually refined, with pick uncertainty estimates based on signal to noise ratio and clarity of phase onset, with emergent arrivals generally having larger errors.The dominant frequency of most earthquake records tends to lie between 8 and 12 Hz, and prefiltering between ∼0.5-20 Hz effectively removes the long-period noise.After the application of filtering, unrotated 3-component seismograms were used for picking.
To carry out our 3-D tomographic imaging, we utilized the FMTOMO package (Pilia et al., 2013;Rawlinson & Urvoy, 2006), which adopts an iterative non-linear framework to jointly constrain hypocenter location and 3-D Vp and Vs structure.The forward problem of travel time prediction is solved using the Fast Marching Method (de Kool et al., 2006), a robust grid-based eikonal solver.The data set we use contains a similar distribution and quantity of P-and S-wave picks, thus making the direct calculation of Vp/Vs from separate Vp and Vs models a relatively reliable undertaking.The initial 1-D velocity model (Figure S1 in Supporting Information S1) used in the inversion is obtained with VELEST (Kissling et al., 1994).This is done by first generating 100 random synthetic models by adding noise with a uniform distribution to the 1-D reference model of Greenfield et al. (2016), before feeding them into VELEST.From the results, we chose the average of 30 best-fitting VELEST output models as the starting model in the tomography (see Figures S1 and S2 in Supporting Information S1).Although the new model is broadly similar to the starting model, it generally exhibits lower velocity in the shallow crust, which is where much of the new data is concentrated.
Inversion for velocity structure and hypocentral location are done sequentially, with a regularized subspace scheme used to adjust velocity parameters (Kennett et al., 1988) and a fully non-linear grid search used for source location that takes advantage of the grid-based eikonal solution to achieve computational efficiency.Damping and smoothing are applied to explicitly regularize the velocity inversion process with appropriate parameter values determined by inspecting the tradeoff (Figure S3 in Supporting Information S1) between the variance and roughness of the model and the data misfit.The optimum values produce a model that is as smooth and least perturbed from the starting model as possible, while still adequately satisfying the data observations.
To complement the results from seismic tomography, deformation at Askja was measured through the utilization of InSAR on Sentinel-1 satellite images.Four tracks were subjected to analysis: two on ascending orbits (tracks 9 and 147) and two on descending orbits (tracks 45 and 111).Due to snow cover throughout the rest of the year, only summer images were taken into consideration.The times-series analysis of deformation was carried out as follows: First, the images were pre-processed to create a coherent stack using the InSAR radar Scientific Computing Environment (ISCE2) software (Rosen et al., 2012).Subsequently, an in-house implementation of the Small Baseline method (Berardino et al., 2002) was used to generate the time series.To mitigate topography effects, a 10-m Digital Elevation Model from Landmaelingar Íslands, the National Land Survey of Iceland, was employed.From the resulting time-series data, two key parameters were derived: the total deflation of ∼70 mm during the period from 2015 to 2020 and the total uplift of ∼650 mm that occurred between 2021 and 2023.

Validation via Synthetic Tests
A series of synthetic tests have been carried out to validate the recovered tomographic models.The first test (Figure S4 in Supporting Information S1) shows the result of inverting two separate checkerboard patterns defined by anomalies of different sizes with noise added and relocation turned on; this most closely replicates the inversion of our observational data set.The added noise is assumed to have a Gaussian distribution with a standard deviation equal to the average estimate of the pick errors: 0.018 s for P-wave arrival times and 0.057 s for S-wave arrival times.We also generate several end-member synthetic models (Figure S5 in Supporting Information S1), which provide insight into how the separate addition of noise and relocation contributes to the degradation of the checkerboard recovery.Given that the real noise level is unknown, and that the checkerboard pattern exerts a different influence on the relocation (and hence reconstruction) compared to the Askja model, we feel that it is useful to illustrate how our assumptions and choices affect our synthetic test results.As expected, the result with no noise and no relocation provides the best reconstruction, with the separate effects of noise and relocation (noting that it is fully non-linear) contributing roughly equally to the decrease in accuracy of the recovered model.Based on the checkerboard test results (Figure S4 in Supporting Information S1), the recovery is generally good at depths above ∼9 km bsl, with a tendency to smear at greater depths.The shortest scale length of structure that can be recovered is ∼4 km horizontally and ∼2 km vertically down to ∼9 km bsl.As expected, the recovery improves significantly with the larger checkerboard size of ∼8 km horizontally and ∼4 km vertically.
A spike test (Figure 2d) is also used to investigate whether nearby anomalies can be distinguished from each other, something that is very important if we are hoping to detect adjacent magma chambers.In the reconstruction shown in Figure 2d, both noise and relocation are implemented in the inversion.The spike test result indicates that it is possible to distinguish between low-velocity anomalies at ∼0.5 and 6 km bsl.There is little evidence of vertical smearing, suggesting that our ability to discriminate between low-velocity anomalies at different depths beneath the caldera is robust.It is worth noting that the amplitude of recovery is generally lower than the input amplitude used to define the synthetic checkerboard and spike models, which is to be expected when damping and smoothing regularization is used in the inversion.

Results and Discussions
Following the application of FMTOMO to the Askja P-and S-wave data set, the final models produce a variance reduction of 90% for P-arrival times and 88% for S-arrival times, resulting in a substantial improvement in data fit relative to the starting model.The corresponding reduction in the root-mean-square residuals is from 0.20 to 0.06 s for P-arrival times and from 0.34 to 0.12 s for S-arrival times.Information on the differences between the initial and final locations is shown in Figures S6-S8 in Supporting Information S1, which illustrate before and after hypocenter locations projected in depth (Figure S6 in Supporting Information S1) and in vertical section (Figure S7 in Supporting Information S1).The histogram of relocation distances (Figure S8 in Supporting Information S1) indicates that the vast majority of events were relocated by less than 2 km.
Although the model region extends well beyond the caldera, in this study we will focus primarily on structure in the mid-upper crust beneath the caldera.In the shallower crust, the volcanic plumbing system exhibits a range of wave-speed anomalies.Particularly noteworthy is a cylindrical anomaly that extends beneath the NW region of the crater lake-approximately the same location as the center of the recent reinflation episode-characterized by lower Vp, lower Vs, and higher Vp/Vs (Figure 2).The presence of two distinct low-velocity anomalies inside the cylindrical anomaly is evident when visualizing the results: one is located at a depth of ∼0.5 km bsl (the "I" in the Vs model in Figure 2), while the other is located at ∼6 km bsl (the "II" in the Vs model in Figure 2).Our findings bear some similarities with previous work.For instance, an ambient noise tomography study has also identified a lower Vs anomaly at around sea level (Volk, 2021), which is likely to be indicative of partial melt, though the lateral resolution is lower than our body wave study, and the depth extent is less well constrained.Similarly, a previous study (Greenfield et al., 2016) reported a lower Vp, lower Vs, and higher Vp/Vs anomaly at a depth of 6 km bsl, which is also likely to be associated with partial melt.The clear delineation of the shallow low-velocity and higher Vp/Vs anomaly can be explained by an accumulation of partial melt, which is supported by multiple complementary observations.First, the reduced levels of seismicity (Figure 2) within the low-velocity cylindrical anomaly hints at the presence of molten material, which may serve as a barrier against stress accumulation and failure.Second, the combination of reduced Vp, reduced Vs, and an elevated Vp/Vs is a characteristic of partial melting (Walck, 1988).Third, geodetic modeling suggests that there is a subterranean magma reservoir located ∼2 km bsl assuming an elastic half-space, although depth resolution is limited (de Zeeuw-van Dalfsen et al., 2012;Pagli et al., 2006).Finally, the observed increase in gravity points to an augmentation in subsurface mass likely caused by melt injection, which may have contributed to the inflation observed since August 2021 (Koymans et al., 2023).If we assume that anomalies I and II in Figure 2 are due to melt, then it is possible to estimate the melt fraction, subject to several caveats (Paulatto et al., 2022).Chief among them is that traveltime tomography tends to significantly underestimate anomaly amplitudes and smear them out in potentially multiple directions.This can lead to an underestimate of the melt fraction and an overestimate of the size of the melt region.Here we estimate the peak melt fraction by applying the Hashin-Shtrikman-Walpoe model that relates seismic velocity to melt fraction, as described by Lyakhovsky et al. (2021), which is suitable for basaltic systems.For an average Vp/Vs ratio of 1.73 at 2 km bsl, where Vp ∼5 km/s, an increase to ∼1.82 beneath the edifice corresponds to a melt fraction of about 12-15%, which is slightly larger than the ∼10% found by Greenfield et al. (2016) for the deeper anomaly.Spike tests (Figure 2d) have provided validation for the separation of two distinct low-velocity anomalies within the cylindrical low-velocity zone.The dynamic interaction between them is likely responsible for the observed inflation and deflation phenomena.The drainage, cooling, and contraction of magma results in deflation, while influx from depth causes inflation.
In light of the recent inflation that started in 2021, perhaps the most intriguing aspect of our results is that the shallow low-velocity anomaly appears to lie directly beneath the region of surface uplift, as defined by both InSAR and GPS measurements, rather than the crater lake or region to the north where the 1961 eruption occurred.Figure 3 highlights this correlation by superimposing the cumulative uplift since August 2021 on top of the shallow low-velocity and high Vp/Vs anomaly.Also shown is the cumulative deflation for 5 years preceding the switch to inflation, which is centered close to the inflation anomaly, but offset slightly to the northeast toward the approximate location of the 1961 eruption.The comparison between the tomography and InSAR appears to suggest that the shallow magmatic processes that are giving rise to inflation are not some lateral and potentially localized offshoot but are directly related to the primary zone of melt migration from depth.Recent geodetic modeling by Parks et al. (2024) suggests that the inflation source is a sill located some 2.5-3.1 km beneath the surface (so <2 km bsl), in a horizontal and vertical location that substantially overlaps with anomaly I. Since anomaly I pre-dated the inflation, it cannot simply be a smeared-out image of the recently formed sill, but instead is a larger body of partial melt that has likely experienced a fresh influx from depth.
The cylindrical low-velocity anomaly is surrounded by a zone of elevated Vp and Vs but reduced Vp/Vs (Figure 2), which is likely to be a signature of consolidated intrusive mafic rocks that have cooled (Brandsdóttir  et al., 1997;Gauntlett et al., 2023).Beyond the high-velocity anomalies that surround the caldera in the shallow crust, a significant portion of the region, particularly where there is a high density of seismicity, is marked by lower Vp, lower Vs, and slightly higher Vp/Vs, which gradually diminishes with increasing depth (see label "III" in the Vs cross-section in Figure 2b).This anomaly is distinctly different from anomalies I and II in Figure 2, which feature a much higher Vp/Vs owing to the large decrease in Vs, but only a moderate decrease in Vp.Based on microseismicity and anisotropy studies, this outer region exhibits significant tectonic earthquakes, tectonic fractures, and fissure swarms (Bacon et al., 2022;Greenfield et al., 2020;Winder, 2021).Anomaly III is mainly concentrated in the highly fractured zones, between Askja and Herðubreið in this case, which is similar to observations in the San Andreas fault and Eastern California Shear Zone (Berg et al., 2021;Fang et al., 2018;Lin et al., 2007).In particular, Vs is found to be sensitive to highly fractured rocks (Catchings et al., 2020), but this diminishes with depth as cracks begin to close, which is also supported by studies of seismic anisotropy (Bacon et al., 2022).
Our new tomographic results provide fresh insight into the magmatic plumbing system of Askja, and allow for comparisons with other rift-related volcanoes in Iceland, for instance, Krafla and Grímsvötn.The existence of a shallow magma body coupled with a cycle of deflation and inflation has also been suggested beneath Krafla and Grímsvötn (Sturkell, Einarsson, et al., 2006).Krafla is another large central volcano that lies in the NVZ, about 80 km north of Askja.Seismic tomography reveals a low-velocity anomaly ∼2 km bsl, interpreted to be a magma body, beneath the approximate center of Krafla (Schuler et al., 2015).Grímsvötn, the most active volcano in Iceland, lies beneath the Vatnajokull icecap.It features a low-velocity anomaly, which is surrounded by highvelocity anomalies in the shallow crust (Alfaro et al., 2007).Geodetic modeling suggests a magma body exists at ∼1 km bsl, and Grímsvötn is similar to Askja and Krafla in that the pressure change of these upper-crustal magma bodies drives the alternation between deflation and inflation (Hreinsdóttir et al., 2014).This limited comparison indicates that shallow melt accumulation beneath active rift volcanoes in Iceland appears common and that the behavior of these magma bodies can be directly related to surface deformation.In future work, a more comprehensive understanding can be achieved by including attenuation tomography, since Q is a very sensitive indicator of partial melt (Wang et al., 2017).Furthermore, it is crucial to acknowledge the limitations in resolution, particularly at depths exceeding 10 km, primarily due to the predominance of shallow seismicity.To overcome this obstacle, a joint inversion of local earthquake and teleseismic body wave data (Huang et al., 2015) has the potential to enhance the resolution through the crust and into the uppermost mantle, offering valuable insights into the shallow and deep plumbing system.

Conclusions
This paper takes advantage of a large new data set of local earthquakes beneath Askja caldera in Iceland to produce a more detailed picture of the shallow plumbing system than has previously been possible using seismic imaging.The main finding of the paper is that there are two distinct lower Vp, lower Vs, and higher Vp/Vs anomalies centered at depths of ∼0.5 km bsl and ∼6 km bsl beneath Askja caldera.Based on the lack of seismicity within these regions, they likely represent zones of partial melt or magma chambers that form part of the underlying magmatic system.In addition, the detection of the shallower anomaly is consistent with geodetic studies that have modeled melt influx at a similar depth.However, our results represent the first 3D imaging of this key part of the magmatic system.Critically, the shallower anomaly sits almost directly beneath the region of surface uplift that began in August 2021, suggesting a possible direct connection.Higher Vp and Vs anomalies that encircle the low-velocity anomaly and align with the layout of the caldera are interpreted as solidified intrusive rocks.Beyond the caldera, lower Vp and Vs and slightly higher Vp/Vs are observed in the shallow crust, matching the locations of fissures and fractures on the surface.

∼3
km wide and extending to ∼8 km bsl, has been imaged beneath Askja caldera • Within this zone, two distinct lower velocity and higher Vp/Vs anomalies are detected at ∼0.5 km and ∼6 km bsl • The shallower anomaly underlies an area of recent surface uplift, which is consistent with shallow melt accumulation Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.Map of the study region showing seismic stations as purple triangles and earthquakes (in the period 2006-2022) as dots colored by depth: shallow (above 3 km bsl-blue), intermediate (between 3 and 7 km bsl-green), and deep (below 7 km bsl-red).The dashed red line encompasses the Askja central volcano.The white box denotes the location of the imaged region shown in Figure 2. The golden box denotes the lateral extent of the full tomographic model.Inset: Tectonic map (Einarsson & Saemundsson, 1987) of Iceland showing the location of fissure swarms (in lighter gray), the direction of plate spreading, and the location of the study region (orange box).The three volcanic zones are the Western Volcanic Zone (WVZ), the Eastern Volcanic Zone (EVZ), and the Northern Volcanic Zone (NVZ).The green lines denote transform zones that accommodate relative movement between the mid-Atlantic Ridge and the rift zones in Iceland.

Figure 2 .
Figure 2. The top panels show a depth slice through the (a) Vp model, (b) Vs model, (c) Vp/Vs model, and (d) Vs spike test at 0.5 km bsl and the corresponding bottom panels show a west-east transect along the dashed line superimposed on the depth slices.The black dots are earthquakes that lie within 3 km of each slice.The synthetic test results, with circles marking the positions of input low-velocity anomalies of peak amplitude 6%, use a colormap to illustrate the recovery with noise added and relocation turned on.The topography is scaled up by a factor of three above the transects to highlight the location of the caldera.

Figure 3 .
Figure 3. Horizontal Vs slice at 0.5 km bsl with estimates of vertical deformation from InSAR -Left for subsidence, Right for uplift.The star marks the approximate location of the 1961 eruption.Note that Line-of-Sight velocities from four different tracks/orbits were used to decompose the signal to vertical deformation.