2021–2023 Unrest and Geodetic Observations at Askja Volcano, Iceland

Unrest began in July 2021 at Askja volcano in the Northern Volcanic Zone (NVZ) of Iceland. Its most recent eruption, in 1961, was predominantly effusive and produced ∼0.1 km3 lava field. The last plinian eruption at Askja occurred in 1875. Geodetic measurements between 1983 and 2021 detail subsidence of Askja, decaying in an exponential manner. At the end of July 2021, inflation was detected at Askja volcano, from GNSS observations and Sentinel‐1 interferograms. The inflationary episode can be divided into two periods from the onset of inflation until September 2023. An initial period until 20 September 2021 when geodetic models suggest transfer of magma (or magmatic fluids) from within the shallowest part of the magmatic system (comprising an inflating and deflating source), potentially involving silicic magma. A following period when one source of pressure increase at shallow depth can explain the observations.


Introduction
Askja volcanic system, in NVZ in Iceland, comprises both a central volcano and a fissure swarm covering ∼190 km × 20 km (Sigmundsson et al., 2020).The central volcano includes a series of nested calderas.Eruptions have been both basaltic effusive and silicic explosive, with the former much more common.Four eruptions occur on average per century (Hartley et al., 2016).The last eruption in 1961 was predominantly effusive and produced a 0.1 km 3 lava field (Blasizzo et al., 2022;Thorarinsson & Sigvaldason, 1962).The last plinian eruption occurred in 1875 (Sigurdsson & Sparks, 1981;Sparks et al., 1981).This major event and continued subsidence following the end of the eruption formed the most recent caldera, now filled with lake Öskjuvatn (up to 220 m deep, Gylfadóttir et al., 2017) (Figure 1).
Measurements of ground movements at Askja began in 1966.Yearly leveling measurements were undertaken between 1966and 1972and from 1983until present. Observations from 1968 to 1970 display a relatively flat trend, changing to inflation from 1970 to 1972.Since 1983 until July 2021 persistent subsidence occurred at Askja.The yearly height difference between the end points of a 1.2 km-long profile decayed in an exponential manner from 12.6 mm during 1983-84 to 1.7 mm during 2020-July 2021 (Sturkell et al., 2006(Sturkell et al., , 2023;;Tryggvason, 1989).Suggested explanations for the long-term subsidence include magma cooling and contraction, or magma withdrawal.Many models based on observed ground deformation and micro-gravity changes have been published: a single contracting shallow point source centered within the caldera (Sturkell & Sigmundsson, 2000;Tryggvason, 1989), a point source and dike model (Rymer & Tryggvason, 1993), an ellipsoidal shallow source (Pagli et al., 2006), and a model with two contracting point sources of various depth as well as horizontal plate spreading (Sturkell et al., 2006).Numerical models including rheological variations in the upper crust suggest an important role of the tectonic setting of Askja in creating observed long-term subsidence (de Zeeuw- van Dalfsen et al., 2012;Pedersen et al., 2009).A positive micro-gravity fluctuation between 2007 and 2009, without changes to the subsidence pattern has been explained by compressibility and magma inflow into a pre-existing magma reservoir (de Zeeuw-van Dalfsen et al., 2013).
The deformation at Askja changed at end of July 2021 (Koymans et al., 2023).Inflation was detected on continuous GNSS (cGNSS) station OLAC (in the Askja caldera) and on Sentinel-1 interferograms (Figures 1 and  2), coinciding with a minor increase in seismic activity.Initially, OLAC was the only existing cGNSS station inside of the caldera (installed in 2015), with two additional stations located within ∼5 km (DYNG and MOFC; see Figure 1).In fall 2021, additional continuous stations TANC, JONC, and KASC were installed.This study focuses on this recent period of inflation, providing an overview of the deformation observations to date, the best-fit analytical geodetic models that can describe these displacements and a discussion of likely processes involved and future eruptive scenarios.

Interferometric Synthetic Aperture Radar Processing and Decomposition
Four tracks of the Sentinel-1 SAR satellites regularly cover Askja at 6-or 12-day intervals: two ascending orbits (T118, T147) and two descending (T111, T45).Interferograms were processed using the InSAR radar Scientific Computing Environment (ISCE2) software (Rosen et al., 2012) and Icelandic Meteorological Office (IMO) inhouse software for Small Baseline (SBAS, Berardino et al., 2002) time-series analysis.The resulting four lineof-sight (LOS) displacement time-series were combined and decomposed to give an approximate estimate of the east component of horizontal displacement (near-East) and the vertical (near-Up) displacements (Drouin & Sigmundsson, 2019).Only summertime images were used because of yearly snow cover between end of September and end of June.The 10 m resolution IslandsDEMv1 DEM from Landmaelingar Islands (LMI, 2022) was used to remove topographic effects.Images were multilooked by 6 in range and 2 in azimuth to obtain approximately 30 m × 30 m pixels.Two separate time periods were considered: time period 1 (T1) spanning 26 July-20 September 2021 and time period 2 (T2) from 20 September 2021-10 September 2023 (Figure 2).To correct for additional processes occurring during T2 (plate spreading, glacial isostatic adjustment and the deep deflation signal occurring since 1983-see Section 1 above) we estimated LOS velocities from July 2015 to September 2020 (from the three Sentinel-1 tracks available during this period) and subtracted this from the corresponding LOS deformation grids spanning T2 prior to modeling.We do not apply a correction to T1 data as these are negligible for such a short time interval.

GNSS Processing
A combination of continuous and campaign GNSS measurements were used to map ground deformation.Five field campaigns were conducted around the beginning and after the onset of the Askja inflation; 27-31 July, 15-24 September 2021, 18-27 August 2022, 27-30 June 2023, and 17-28 August 2023, building on annual GNSS campaigns in the region.Each site was occupied for a minimum of 36-hr during each campaign.Data were analyzed using the GAMIT/GLOBK software, version 10.7 (Herring et al., 2018), using over 100 global reference stations to constrain daily site positions in the IGb14 reference frame.During the analysis we scaled position uncertainties using site-specific, elevation dependent, noise values estimated during the processing.In addition to station coordinates, the processing solved for satellite orbits and earth rotation parameters, atmospheric zenith delay every 2 hours, and three atmospheric gradients per day.We corrected for ocean loading using the FES2004 model (Lyard et al., 2006) and used the IGS14 azimuth and elevation dependent absolute phase center model for all antennas.We estimated a velocity solution for July-September 2021 and September 2021-September 2023.In addition, the background deformation was estimated using data from June 2015-January 2021 and subtracted from the T2 displacements prior to modeling.

Geodetic Modeling
Geodetic inversions utilized a modified version of the GBIS software (Bagnardi & Hooper, 2018), assuming deformation sources within a uniform elastic half-space.Input data were the campaign and continuous GNSS measurements and LOS displacements derived from SBAS time-series analysis of interferograms from the four Sentinel-1 tracks.Interferograms were decimated using a higher sampling within 10 km of the main inflation signal (coordinates used: longitude 16.75, latitude 65.03) and a sparser sampling beyond this.Simple source geometries were employed to determine the optimal source type/s and location.Model geometries explored included a point source (Mogi, 1958), a spheroid (McTigue, 1987), a prolate spheroid (Yang et al., 1988) and a rectangular shaped dislocation (Okada, 1985).The geodetic inversions were run separately for T1 and T2, as deformation is different within these time periods (Figure 2).

Previous Subsidence at Askja
Long-term subsidence within the caldera was continuing until summer 2021 (Figures S1 and S2 in Supporting Information S1) albeit at a reduced rate.Decomposition of interferograms from three separate Sentinel-1 tracks shows the location of the deflation signal from 2015 to 2020 (Figure S3 in Supporting Information S1).The maximum subsidence is located in the center of the main Askja caldera.This is in agreement with previous geodetic studies which located the center of the source in a similar location, at depths between 2.5 and 3.8 km (Pagli et al., 2006;Rymer & Tryggvason, 1993;Sturkell & Sigmundsson, 2000;Sturkell et al., 2006).

New Period of Uplift
Our geodetic observations show clearly an inflation signal at Askja volcano, beginning near end of July 2021, centered on the western to north-western edge of Öskjuvatn (Figure 2; Figures S4-S6 in Supporting Information S1).During the initial T1 period (26 July-20 September 2021), the rate of inflation was higher (>70 cm/yr) and the inflation signal was centered very close to cGNSS station OLAC (Figure 2a; Figure S4  The cGNSS station OLAC, had an average uplift rate during T1 of 709 mm/yr, and 282 mm/yr in T2.Transects through the near-Up displacement rate maps clearly show that during T1 a subsidence signal is evident on the eastern side of Öskjuvatn, but not discernible during the latter period T2 (Figures 2a-2d).
Geodetic inversions were run for the T1 and T2 periods.A variety of geometries were tested for both periods, with all initial model parameters free.For both time periods the best-fit single source model (model with the lowest misfit/minimum residuals) was an inflating sill-type source, modeled using an Okada dislocation.For T1, an additional deflating source is required to fit the subsidence signal in the east and to reduce the misfit.The weighted residual sum of squares (WRSS) (e.g., Cervelli et al., 2001;Parks et al., 2020) is reduced from 0.393 for a single source to 0.368 for the dual source model (Figures 2a and 3a; Figures S7 and S8 in Supporting Information S1).During T1, the inflating sill-type source (Okada dislocation) is horizontal with a median length of 5.9 km, width of 3.7 km and strike of 119°.It is located at a median depth of 2.8 km (2.6-3.1 km, 95% confidence interval (CI)) and corresponds to a volume change of (13.6 ± 2.7) × 10 6 m 3 .The deflating source has a medium dip of 4°, length of 5.4 km, width of 0.7 km and strike of 108°.It is located at a medium depth of 3.3 km (3.1-3.6 km, 95% CI) and corresponds to a volume change of ( 10.6 ± 2.5) × 10 6 m 3 (Figure 3a; Figure S9 and Table S1 in Supporting Information S1).
For T2 the best-fit single source was also determined to be a sill-type source with a medium dip of 12°, length of 4.5 km, width of 2.9 km and strike of 126°.It is located at a median depth of 2.7 km (2.5-2.8 km, 95% CI) with a rate of volume change of (13 ± 1.7) × 10 6 m 3 /yr (Figure 3b; Figure S10 and Table S2 in Supporting Information S1).The two shallow sources for T1 and T2 reside at a similar depth and geographical location.Considering the near constant rate of inflation between September 2021 and September 2023 (as indicated by both GNSS and InSAR), the total estimated volume change during T2 is (25.6 ± 3.4) × 10 6 m 3 and for T1 and T2 combined it is (39 ± 6) × 10 6 m 3 .

Discussion
The magma plumbing system beneath Askja has been described as highly stratified from petrological studies, comprising a rhyolitic upper part and ferrobasaltic lower part (Sigurdsson & Sparks, 1981).Sparks et al. (1981) proposed that an ascent of basaltic magma (during rifting) mixed with pre-existing rhyolitic magma triggering the 1875 plinian eruption.Recent petrological studies have suggested Askja 20th century basalts have mixed with evolved melts in the crust (Hartley & Thordarson, 2013).Earthquake tomography by Greenfield et al. (2016), resolved a magma storage region beneath the Askja caldera centered at depths of between 5 and 7 km bsl and also a low Vp/Vs ring around Askja in the shallow crust ∼2 km bsl (including the eastern side of Öskjuvatn Lake), which may correspond to a series of silicic magma bodies.
The total amount of subsidence at Askja between 1983 and 2021, has been estimated at 2.4 m, corresponding to a contraction volume of ∼0.1 km 3 (Sturkell et al., 2023).Previous studies suggested multiple processes were likely responsible for the observed long-term subsidence, including cooling and contraction of a pre-existing magma body, magma withdrawal and viscoelastic relaxation.However, these models struggle to fully explain the magnitude of the observed deflation signal.Another volcano in the NVZ, Krafla, had a rifting episode between 1975 and 1984.There inflation continued for about 5 years after the episode, but in 1989 this changed to subsidence, initially at a rate of ∼5 cm/yr.The onset of this subsidence at Krafla was attributed to contraction and ductile flow of material away from the spreading axis (Sigmundsson et al., 1997).Although extension along the plate boundary and viscoelastic relaxation may play an important role in recent subsidence at Krafla, the same processes may only account for a small fraction of the subsidence observed at Askja (Lanzi et al., 2023).Cooling and contraction of an existing magma body has been proposed as contributing to the deflation at Askja, however crystallization alone is unlikely to account for all the subsidence, especially when considering the low abundance of phenocrysts within the Askja 20th century lavas (<1 vol.%) (Hartley & Thordarson, 2013).Thus, additional processes are required-possibly magma outflow followed by a continued poro-viscoelastic response (e.g., Liao et al., 2021).
For T1 the volume change for the shallow inflating source is (13.6 ± 2.7) × 10 6 m 3 and for the deflating source ( 10.6 ± 2.5) × 10 6 m 3 .The median depth of the inflating source is 2.8 km and for the deflating source is 3.3 km.However, the depth of the sources and the distance between them may be greater, considering that these are analytical models assuming an elastic, homogeneous half-space, and initial finite element models (including the effects of topography and heterogeneity) indicate the source depths could be up to ∼1 km deeper (O'Hara, 2023).The sources extend to the eastern edge of Öskjuvatn caldera (Figures 1 and 3)-a known geothermally and seismically active area for decades (Einarsson & Brandsdóttir, 2021).The location of these sources is in proximity to a possible felsic magma storage region beneath Askja (Greenfield et al., 2016).Silicic melts beneath Icelandic central volcanoes are typically formed by magma bodies heating and remelting the surrounding basaltic crust (Martin, 2006;Muehlenbachs et al., 1974;Oskarsson et al., 1982;Sigmarsson et al., 1991), and the formation of silicic melts is favorable along the edges of magma bodies where the thermal gradient is greater and temperature conditions are conducive (e.g., along the edges of calderas, Jónasson, 1994).
Our modeling indicates the initial inflation may have been triggered by the upward migration of magma from an existing source located at a depth of around 3.1-3.6km to a shallower level.This would be consistent with the analysis of recent gravity observations, which, while the uncertainties are large, can be interpreted in terms of no increase in mass beneath the Askja caldera between summer 2021 and 2022 (Koymans et al., 2023).Since no additional shallow crustal intrusions/inflation events have been observed since 1983, it is feasible that the melt involved was evolved.It is possible that crystallization contributed to lowering of density and buoyant uprise.
For the T2 period, the inflation source is located at a similar depth (median depth of 2.7 km) with a volume change of (25.6 ± 3.4) × 10 6 m 3 (Figure 3), and a similar location.No deflation source is detected during T2.This may indicate recharge of melt from depth since September 2021.It is feasible that the continued inflation observed during T2 may also reflect additional processes related to volatile exsolution following the shallowing of magma during T1.However, the difference in depth of the two T1 sources is less than 1 km, so only minor volatile exsolution would be associated with magma transfer from the deeper to the shallower.Furthermore, no significant changes in gas emissions were detected from either the fumarole or water samples collected at both Lake Víti and Öskjuvatn during the summers of 2022 and 2023 (M.Pfeffer, personal communication, September 2023).It is also possible that at least part of the deformation during this T2 period may be attributed to poro-viscoelastic effects in response to initial magma migration related to the upper sill-type intrusion, if a liquid magma lens injected is surrounded by crystal mush.In that case, poroelastic diffusion and/or viscoelastic relaxation of the mush could cause continued inflation after the intrusion of magma stops (Liao et al., 2021;Sepulveda et al., 2023).
The schematic displayed as Figure 4 illustrates the approximate location of the sources identified in this study and that of the primary magma storage region and low Vp/Vs anomaly identified by Greenfield et al., 2016.Seismicity as recorded by the national seismic network of Iceland, within the entire study area since January 2016 until September 2023, and in subareas, is displayed in Figures S11-S16 in Supporting Information S1.The depths of the earthquakes are centered around 3-5 km.An initial increase in activity occurred between August and November 2021.In the known geothermal area on the eastern side of Öskjuvatn, the initial increase (possibly related to upward migration of pre-existing melt during T1) was followed by a drop in activity to background levels.Earthquake activity under Öskjuvatn, and around the main Askja caldera (excluding the geothermal area east of Öskjuvatn), shows an initial increase in activity coinciding with the onset of unrest, followed by a short hiatus in activity (at the start of November 2021) then an additional increase in activity from January 2022 onwards.The suggested variable seismic response of the caldera region is consistent with the findings of Winder et al. (2023), based on a dense local seismic network, who also observe an initial increase in earthquake rate during August 2021 and strong spatial and temporal variations.Additional analysis of the seismicity (including precise relocations) is required to better understand the temporal and spatial variation and the involved causal processes (e.g., stress triggering, magma migration etc.).
Possible scenarios for future eruption triggering include: (a) Additional magma inflow into a shallow source at 2.5-3.1 km beneath the caldera-increased pressurization leading to an eruption either beneath the lake or within the main caldera.In this scenario the eruption may be either effusive or explosive.(b) Dyke propagation out of caldera.This would likely result in an effusive basaltic fissure eruption which may occur anywhere within the Askja fissure swarm.(c) Interaction/intersection of basaltic magma with a pre-existing silicic magma body, as historic plinian eruptions indicate they were triggered by such interaction.Prior to the 1875 eruption, dikes initially propagated outside the caldera to the north of Askja resulting in basaltic fissure eruptions (Sigurdsson & Sparks, 1978).This effusive activity was followed by a rhyolitic eruption.A similar sequence of events appears to have taken place 11,000 years ago at Askja (Sigurgeirsson, 2016).At the time of writing (January 2024) the inflation is continuing but at a reduced rate since September 2023 (as observed by GNSS).Continued monitoring and additional observations and modeling are required to understand the process responsible for this change.

Conclusions
Since 1983 Askja volcano has been subsiding, albeit at a declining rate until the end of July 2021, when inflation began.We modeled the deformation assuming sources in a homogenous elastic half-space.For the period 26 July-20 September 2021, two separate sources were identified.An inflating source at a median depth of 2.8 km and with a volume change of (13.6 ± 2.7) × 10 6 m 3 and a deflating source at 3.3 km with a volume change of ( 10.6 ± 2.5) × 10 6 m 3 .The sources extend to the eastern edge of Öskjuvatn-a known geothermally and seismically active area, likely related to a relatively "open" part of the caldera faults in the Askja region where geothermal and/or magmatic fluids, including volcanic gases, can migrate.For the period 20 September 2021-10 September 2023 an inflation source is located at similar depth (median depth of 2.7 km) as in the earlier period, with a volume change of (25.6 ± 3.4) × 10 6 m 3 .No shallow deflation source is detected during this time.The modeling results can be interpreted in terms of the beginning of unrest being associated with magma transfer within the uppermost part of the Askja magmatic system.It is possible that silicic magma is involved.Absence of a shallow deflation source during T2 indicates magma recharge from depth and/or a triggered poro-viscoelastic response.The possibility that both silicic and basaltic magma bodies may be residing within proximity to each other beneath the main caldera should be considered when evaluating future eruptive scenarios.

Figure 1 .
Figure 1.Geological map of Askja showing the location of calderas and lava flows.Outlines of lava flows are from Sigvaldason et al. (1992).The cartographic data is from the National Land Survey of Iceland.Locations of GNSS stations are indicated.Outlines of fissure swarms on inset map are from Hjartardóttir and Einarsson (2021).Sill outlines in (c) are derived from this study.

Figure 2 .
Figure 2. InSAR and GNSS displacements.(a) Near-Up displacement rates for time period T1 (26 July-20 September 2021) and (b) T2 (20 September 2021-10 September 2023).(c, d) Transects through near-Up displacement rate maps for T1 (purple lines) and T2 (red lines).Location of traverses are displayed in (a, b).(e) Timeseries of vertical displacements at GNSS stations OLAC and MASK.Black triangle in (a, b) shows location of OLAC and white triangle MASK.
. During T2 (after 20 September 2021), the inflation rate is significantly reduced (∼35 cm/yr) and the center of the signal migrated further to the NE (Figure2b; FiguresS5 and S6in Supporting Information S1).

Figure 4 .
Figure 4. Schematic of Askja volcano magmatic plumbing system.Approximate location of inflation source during T1 is shown by the yellow disk and deflation source magenta disk.Inflation source in T2 is shown by the orange disk.Location of primary magma storage region and low Vp/Vs anomaly is based on work by Greenfield et al. (2016).The pale blue ellipse near the surface represents Lake Öskjuvatn.Note: non-linearity between depths of 10-20 km.