Inflation-Deflation Episodes in the Hengill and Hrómundartindur Volcanic Complexes, SW Iceland

Non-eruptive uplift and subsidence episodes remain a challenge for monitoring and hazard assessments in active volcanic systems worldwide. Sources of such deformation may relate to processes such as magma inflow and outflow, motion and phase changes of hydrothermal fluids or magma volatiles, heat transfer from magmatic bodies and heat-mining from geothermal extraction. The Hengill area, in southwest Iceland, hosts two active volcanic systems, Hengill and Hrómundartindur, and two high-temperature geothermal power plants, Hellisheiði and Nesjavellir. Using a combination of geodetic data sets (GNSS and InSAR; Global Navigation Satellite Systems and Interferometry Synthetic Aperture Radar, respectively) and a non-linear inversion scheme to estimate the optimal analytical model parameters, we investigate the ground deformation between 2017–2018. Due to other ongoing deformation processes in the area, such as plate motion, subsidence in the two geothermal production fields, and deep-seated source of contraction since 2006, we estimate 2017–2018 difference velocities by subtracting background deformation, determined from data spanning 2015–2017 (InSAR) or 2009–2017 (GNSS). This method highlights changes in ground deformation observed in 2017–2018 compared to prior years: uplift signal of ∼10 km diameter located in the eastern part of the Hengill area, and geothermal production-related temporal changes in deformation near Húsmúli, in the western part of the Hengill area. We find an inflation source located between the Hengill and Hrómundartindur volcanic complexes, lasting for ∼5 months, with a maximum uplift of ∼12 mm. Our model inversions give a source at depth of ∼6–7 km, located approximately in the same crustal volume as an inferred contracting source in 2006–2017, within the local brittle-ductile transition zone. No significant changes were observed in local seismicity, borehole temperatures and pressures during the uplift episode. These transient inflation and deflation sources are located ∼3 km NW from a source of non-eruptive uplift in the area (1993–1999). We consider possible magmatic and hydrothermal processes as the causes for these inflation-deflation episodes and conclude that further geophysical and geological studies are needed to better understand such episodes.


INTRODUCTION
Uplift and subsidence episodes without eruptive activity are common in several active volcanic systems (e.g., Etna, Campi Flegrei, Yellowstone, Okmok, Alutu; Italy, United-States, Ethiopia Lima et al., 2009;Chang et al., 2010;Biggs et al., 2011;Walwer et al., 2019). The exact physical origin of the responsible processes (e.g. hydrothermal or magmatic) is often not clear. In volcanic areas, magma movements are often the primary explanation for causing deformation. However, studies have shown that more complex mechanisms, often linked to highenthalpy fluids, volatiles, and their interaction with magma and magmatic gases, can also generate transient uplift motions followed by gradual subsidence (e.g. Lima et al., 2009;Chang et al., 2010;Chiodini et al., 2015). Eruptions are often preceded by accelerated uplift and seismicity (Voight, 1988). Uplift episodes can thus resemble precursory signals for an eruption. It is important for hazard assessments of volcanic systems to better understand the sources of deformation, motivating detailed documentation of inflation-deflation episodes at volcanic systems.
The Hengill area is at the junction of the North American plate, Eurasian plate and Hreppar micro-plate (Einarsson, 2008). The surface structures of the plate boundaries within our area of interest are an oblique rift at the Reykjanes Peninsula (RP), a transform zone named the South Iceland Seismic Zone (SISZ), and the rift zone of the Western Volcanic Zone (WVZ; Figure 1). The topographic high of the area, Mt. Hengillpeaking at 803 m a.s.l.is located approximately 25 km ESE from the capital city of Iceland, Reykjavík (Figure 1). The area is a centre of volcanic and tectonic activity hosting three successive volcanic systems, from oldest to youngest: the extinct Grensdalur volcanic system and the active Hrómundartindur and Hengill systems (Saemundsson, 1967;Foulger, 1988). The last eruption in the area was approximately 2000 years ago (Saemundsson, 1967).
The dominant surface lithology in the Hengill area is basaltic and, with the exception of few intra-glacial lavas, mainly from subglacial eruptions as evidenced by the hyaloclastite pillow breccia, tuff and pillow-lava (Saemundsson, 1967;Arnason et al., 1969). Minor andesitic and rhyolitic units can be found, which may suggest partial recycling of magma in the crust (Trønnes, 1990).
The numerous fumaroles and mud pots found in the Hengill area are manifestations of subsurface high-temperature geothermal systems resulting from young cooling intrusions at depth and permeable lithologies (Stimac et al., 2015). The Nesjavellir geothermal power plant, located north of the Hengill central volcano (Figure 1), has produced hot water for district heating in Reykjavík since 1990. In 2006, the Hellisheiði power plant, located southwest of Hengill (Figure 1), started thermal and electric production (Franzson et al., 2010). The Hellisheiði power plant was expanded in 2016, with new production wells in Hverahlíð (SE of Hellisheiði, Figure 1). As of end of 2019, the Hellisheiði power plant has an installed capacity of 303 MW in electricity (MWe) and 200 MW in thermal energy (MWth), whilst Nesjavellir has an installed capacity of 120 MWe and 290 MWth. The extraction of geothermal fluids by the power plants causes localised subsidence of up to approximately 2.5 cm/yr in the production areas (Juncu et al., 2017). Injection of geothermal fluids in Húsmúli (Figure 1), that started in September 2011, triggered a 4 month episode of localized uplift and ∼20,000 earthquakes (Juncu et al., 2020). The sources of the anthropogenic subsidence and uplift signals, as well as most of the associated seismicity, reach 2.5-3 km depth, approximately the maximum depth of the geothermal wells (Juncu et al., 2017(Juncu et al., , 2020. In addition to tectonic and anthropogenic deformation, the Hengill area displays episodes of deformation caused by deep sources (>5 km), which are the focus of this study. During 1993During -1999 an unrest episode took place in the SE part of the Hengill area. The uplift rate was around 2 cm/yr, observed using GPS, levelling (Sigmundsson et al., 1997;Hreinsdóttir, 1999) and InSAR data (Feigl et al., 2000). A simple point source of pressure within an elastic halfspace model (Mogi, 1958) was used by Feigl et al. (2000) to constrain a source at 7 ± 1 km depth with a volume increase of 3.9 × 10 6 m 3 per year, which was located approximately 3 km north of the town of Hveragerði ( Figure 1). This depth coincides broadly with the brittleductile boundary of the Hengill area, estimated using maximum earthquake depth, approximately 6-7 km below the surface, thinning towards the western end of the Reykjanes Peninsula and thickening through the eastern part of the SISZ (Foulger, 1995;Li et al., 2019). This 6 year uplift episode was associated with new surface fractures (Clifton et al., 2002) and intense seismic activity (with earthquakes reaching M L > 5) in the whole area of Hengill (Sigmundsson et al., 1997;Feigl et al., 2000;Li et al., 2019;Parameswaran et al., 2020;Blanck et al., 2020). A seismic tomography study of the 1993-1999 uplifting volume from Tryggvason et al. (2002), shows a lower ratio of P-to S-wave velocities in its estimated source location, which is interpreted as the presence of supercritical fluids at these depths rather than a large partially melted magmatic body. Recent seismic and resistivity studies seem to agree with the findings of the aforementioned seismic tomography (Jousset et al., 2011;Gasperikova et al., 2015;Li et al., 2019). The integration of resistivity models based on electromagnetic data and seismic tomography highlights a layer with low electric resistivity anomaly and high Vp/Vs at 2.5-4 km depth, which possibly extends deeper (Jousset et al., 2011;Árnason et al., 2010). Under the center of a 2006-2017 subsidence episode (∼3 km NNW of the 1993-1999 inflation center; Figure 1), the resistivity structure is different. That area shows a high resistivity core layer around 2-5 km depth (Gasperikova et al., 2015).
Between 2006 and 2017, Global Navigation Satellite Systems (GNSS) and Interferometry Synthetic Aperture Radar (InSAR) time series indicate a broad-scaled subsidence of the Hengill area (Juncu et al., 2017). Their study indicates a source of contraction located at ∼5-6 km depth, in the vicinity of Ölkelduháls, between the Hengill and Hrómundartindur volcanic systems ( Figure 1) and an estimated volume decrease of 2.4 × 10 6 m 3 per year. Juncu et al. (2017) discuss the possible magmatic and hydrothermal contracting processes behind the observed subsidence. The overall natural and man-made ground deformation of the Hengill area appears to induce heightened strain to the neighbouring SISZ (Árnadóttir et al., 2018).
The Hengill area is a location of interacting anthropogenic, magmatic, and tectonic processes, which are manifested at the surface by spatiotemporally complex ground deformation. Therefore, small-scale deformation transients in the Hengill area may be obscured by the multiple deformation processes. In this study, we analyse geodetic GNSS and InSAR data sets over the Hengill area to investigate transient deformation in 2017-2018. This research was sparked from observations in our continuous GNSS times series displaying, amongst other signals, an unexpected uplift in the area ( Figure 2). We explore the possible origins and source geometry of the uplift using analytical models. We further discuss their possible relation to our current geological and geophysical knowledge of the area and compare our findings with seismicity and borehole geothermal data in the uplift region.

Geodetic Data
We map the crustal deformation over the Hengill area by using a combination of GNSS and InSAR methods (GNSS and InSAR Data). The GNSS and InSAR data sets are complimentary: while the InSAR data allow for excellent ground coverage, the campaign GNSS data capture better 3D motion and the continuous GNSS network (cGNSS), with six stations in the Hengill area (Figure 1), helps constrain the start and end times of the 2017-2018 uplift episode. To simplify our study of the transient ground deformation between 2017 and 2018, we calculate the change in deformation relative to previous years. In practical terms, this means we subtract the pre-inflation velocity field from the data covering the inflation period. This approach allows us to cancel out what we assume to be steady background deformation (i.e., plate motion and subsidence at Nesjavellir, Hellisheiði, and the deep source of contraction). Subtracting the preinflation deformation field, however, causes some complications, which we address specifically in Isolating and Modelling the 2017-2018 Uplift Signal.

GNSS
We analyse of a total of 62 GNSS time series, composed of data from seven continuous sites (HVER, GFEL, OLKE, HLID, NVEL, HUSM and SELF; see detrended time series in Figure 3) and 55 campaign sites located within approximately 20 km of the main study area (Figure 1). Campaign GNSS data are typically acquired annually in early summer. Each benchmark is usually measured for at least 48 h, in order to obtain accurate position; however, a subset of data from 2009, 2012, 2016, and 2018 were measured for a shorter time. The campaign and continuous GNSS data between 2009 and 2019 are processed using the GAMIT/ GLOBK suite of programs (version 10.6, Herring et al., 2015) to estimate daily positions in the ITRF14 reference frame (Altamimi et al., 2016). We solved for satellite orbits and Earth rotation parameters, estimating atmospheric zenith delay every 2 hours. We corrected for ocean loading using the FES2004 model (Lyard et al., 2006) and applied IGS14 azimuth and elevation dependent phase center variation model for all antennas (Montenbruck et al., 2015;Schmid et al., 2016). For details regarding the analysis of the GNSS campaign data and processing steps, the reader is referred to (Árnadóttir et al., 2018, and references therein).
The deformation associated with the transient inflation in 2017-2018 can be seen as northward motion of station NVEL, southeastward motion at HVER, and eastward motion at OLKE (Figure 3). The inflation signal is, however, not as obvious in many of the campaign GNSS time series (Figure 4). To estimate deformation from the inflation period, we start by correcting the time series for instrumental offsets resulting from antenna changes (see Supplementary Material for more details). We then delete outlier data, defined as points with position offset by more than 15 mm in the horizontal components (for cGNSS 10 mm) and more than 20 mm in the vertical components relative to the other positions recorded within the same campaign. We also remove data points with standard deviation above 10 and 30 mm in the horizontal and vertical components, respectively. We further remove data that are clearly affected by post-seismic motion following the 2008 SISZ earthquakes (Figure 1) or local transient motions in the Húsmúli area (e.g. 2011 uplift; Juncu et al., 2020), whose characterizations are outside the scope of this study. For cGNSS stations we estimate and remove seasonal motion (annual and biannual) in each component. This is done using a least-square fit of cyclic components (Blewitt and Lavallée, 2002), using time series between 2012 and 2017, or shorter if the data of a station do not cover this time span (e.g. NVEL, Figure 3).
We use the Tsview module from the GGMatlab toolbox (version 2.02; Herring, 2003;Herring et al., 2015) to estimate the displacements between 2017 and 2018 for the horizontal components of each of our cleaned campaign and continuous GNSS time series. We found that the vertical component was too noisy to allow for reliable displacement estimates. For stations with more than 30 data points, such as cGNSS and most frequently observed campaign stations, the displacement uncertainty was estimated with the Real Sigma option (Herring, 2003); otherwise, the uncertainty estimate assumes uncorrelated (white) noise. The velocity of each station is calculated, allowing for an offset (i.e., the displacement we are after) in the time series at the time of the uplift. Our estimated horizontal displacements range from ∼0-15 mm (Supplementary Material) and show an overall outward radial pattern from the Ölkelduháls region ( Figure 4).

InSAR Data
We use all available Sentinel-1 (S1; C-band sensor) and TerraSAR-X (TSX; X-band sensor) data over the Hengill area, which consist of three different tracks from the S1 constellation and data from a single track of the TSX satellite ( Table 1). Snow-coverage limits useful images to the summer months, from early June until late September of each year. The relatively short and small-amplitude uplift signal of interest ( Figure 5) is partly overprinted by other ongoing processes, such as plate motion and deformation in the geothermal production areas. We thus perform a two-fold time series analysis for each of the tracks. First, we estimate the mean velocity between 2015 and 2017 (time period prior to the uplift episode). Second, we estimate the mean velocity between 2017 and 2018 using the same methodology as the analysis for 2015-2017. As the focus of our study is the 2017-2018 uplift episode seen in the Hengill area, we correct for the long term and local signal from the aforementioned local and broad scale deformation by subtracting the mean velocities of 2015-2017 from the ones in 2017-2018. As the number of SAR images at our disposal depends on the constellation used (more than 20 images per year for Sentinel-1 and two to four images per year for TerraSAR-X), we adopt different methodologies for each of the time series analysis as detailed in the section 2.1.2.1. and 2.1.2.2. For both S1 and TSX processing, we use the intermediate resolution TanDEM-X digital elevation model (Rizzoli et al., 2017) to remove the topographic phase contributions.  Our S1 InSAR data set is generated from images acquired in interferometric wide-swath mode provided by the European Space Agency (ESA), from June 2015 to September 2018. Three tracks, two ascending (T16, T118) and one descending (T155), cover our area of interest. Interferograms were generated using the InSAR Scientific Computing Environment software (ISCE, version 2.2.0; Rosen et al., 2012), all referred to a single image in 2017 for each track ( Table 1). The topsApp function of the software initially uses the single look complex bursts from the Terrain Observation with Progressive Scanning SAR (TOPSAR) acquisition of the S1 constellation. The bursts are merged before coarse and fine co-registration and formation of filtered wrapped and unwrapped interferograms. The interferograms generated for each track are then divided in two time spans for the generation of our time series; 2015-2017 and 2017-2018, and cropped to our common area of interest. We use the cropped stack of common referenced interferograms to visually select interferograms suitable for time series generation, removing interferograms with clear atmospheric disturbances or missing bursts. The time series analysis method is the same as described by Drouin and Sigmundsson (2019). Our reduced selection of interferograms is then referred to a common area near the GFEL GNSS station (location in Figure 1 and red time series in Figure 3) which shows little to no deformation that can be associated with local signals due to transient anthropogenic and broad-scale uplift and subsidence ground motions of the Hengill area (Juncu et al., 2017). Our stacks of interferograms (cropped and referenced to a common area) are then inverted to estimate the mean constant velocity in each pixel, accounting for their FIGURE 5 | InSAR mean LOS velocities for TerraSAR-X (TSX) track T41 and Sentinel-1 (S1) tracks T16, T118 and T155 (Table 1)  Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 725109 weight estimated from the coherent maps. Pixels with a coherence below the threshold of 0.8 (1 being a perfectly coherent area and 0 being a non-coherent area) are judged unreliable and not taken into account in the velocity estimate. The time series analysis minimizes the effect of transient atmospheric disturbances. Figure 6 shows the resulting mean LOS velocities, in which we masked out areas with a normalized standard deviation above 0.7 mm/yr. This masking particularly affects the junction of the Hengill and SISZ, where more vegetated and wetland surfaces introduce noise in our mean LOS velocity estimates.

TerraSAR-X
To complement our ground deformation analysis of this uplift episode, we use a third independent data set which is our longterm time series (2009-2019) of SAR images from the high resolution X-band TSX satellite from the German Aerospace Center (DLR). TSX images from the ascending track T41 between 2015 and 2018 are limited to 2-4 StripMap SAR acquisitions per year. We adopt a different methodology than described in InSAR Data to exploit the higher signal-to-noise ratio scatterers from the reduced number of images at our disposal for our time series analysis. We thus use the Small-Baseline time series analysis available from the Stanford Method for Persistent Scatterers software (StaMPS, version 3.3b1; Hooper, 2008;Hooper et al., 2012) to generate short-temporal baseline interferograms. In this methodology, interferograms, referenced to a single image (September 2017, Table 1), are generated from single-look complex images generated using the DORIS software (Kampes et al., 2003). The small-baseline processing then generates interferograms which minimize the perpendicular baseline and the time between acquisitions. This method maximizes the correlation in the interferogram stack and therefore allows to find more permament scatterers. They are used for time series estimations on this track between 2015 and 2018. The StaMPS processing further allows estimation and removal of errors, for example orbital ramps and spatially correlated look angle errors to limit the propagation of these errors to our final InSAR time series. To avoid phase jumps in our unwrapped mean LOS velocities, the default value of the spatial unwrapping grid is set to 100 m. We minimize the introduction of noise and decorrelated scatterers from the complex topographic escarpments and vegetated area E-SE of our study area by using the weeding option of the StaMPS software. Table 1 details the main parameters and information for this track. The mean LOS velocities from 2015 to 2017 and 2017-2018 we obtain from our analysis are presented in Figure 6.

Isolating and Modelling the 2017-2018 Uplift Signal
The differential InSAR and GNSS velocities (Figure 4 and Figure 6) highlight two apparent sources, which henceforth we refer to as Húsmúli for the localized subsidence, NW of the Hellisheiði geothermal field, and Ölkelduháls for the broadscale uplift in eastern Hengill (named after the closest topographic high). The footprint (i.e. maximum amplitude location and radial wavelength) of the latter signal is similar to Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 725109 8 the eastern Hengill contracting source between 2006 and 2017 described by Juncu et al. (2017) and observable in the vertical time series of the cGNSS station OLKE (Figure 2). Initial inversions of our data sets implied that a similar source location could explain the subsidence during 2006-2017, and the uplift episode in 2017-2018. Physically, a volume cannot simultaneously inflate and deflate. It is thus necessary to isolate the transient displacements associated with the uplift in our geodetic data sets. In our prior approach (Geodetic Data), the data used to correct for background deformation contain longer term subsidence signal. If we assume that the subsidence was not on-going during the period of uplift, then we need to evaluate the time-span of the uplift in our 2017-2018 data (using cGNSS time series, see Uplift Time Span). The estimated time-span of uplift is then used to correct our geodetic displacements for the time-span in which the area did not subside in 2017-2018 (5 months; from approximately November 2017 until March 2018). We calculate the ground motion using the source of contraction (model of Juncu et al., 2017) for each of the InSAR LOS and GNSS horizontal components used in our study. The predicted displacements are then scaled relative to the 5 months time span when no subsidence occurred in 2017-2018, and added to our previously corrected data sets ( Figure 6 illustrates the data sets before and after correcting for the deep source of contraction).

Uplift Time Span
To determine the time span of the uplift, we use the cleaned linear and seasonal detrended time series from the two cGNSS stations OLKE and HVER (location in Figure 1; details on the GNSS data in GNSS) from January 1, 2012 until June 30, 2019. Both stations are proximal (within ∼5 km) to the main uplifting area and show clear motions in horizontal and vertical components associated with the uplift episode ( Figure 3). The North component from the OLKE station, located ∼1 km due east of the center of the uplift, is not included in this analysis as it shows little to no motion in the North direction within this time span ( Figure 5). In order to study the time-span of the uplift we fit a linear piecewise function (Jekel and Venter, 2019), with the fixed number of three segments (i.e. prior linear trend, uplift motions and subsequent linear trend) to the north, east and vertical components of the selected and detrended cGNSS data. We use a custom optimization algorithm that includes the individual uncertainties of our data points in our likelihood function which is then used by the algorithm to retrieve the optimal breakpoints of the piecewise function. Additional details on the method are in the Supplementary Material. The vertical component has larger uncertainties and is more affected by variable seasonal signals than the horizontal motions observed by GNSS. Therefore, we estimate the time interval of uplift using only the horizontal displacements. The optimisation results show that the uplift lasted for 5 months.

Inverse Modelling
We attempt to explain the observed deformation with geometric sources embedded in uniform elastic half space to reproduce the observations. The inflation trends appear quite linear with sudden stop, with no relaxation of trends observed as could be indication of poro-elastic or viscoelastic deformation ( Figure 5). Furthermore, the time scale of inflation is short (5 months). For these reasons we use elastic host rock rheology (shear modulus G 10 GPa and Poisson's ratio ] 0.25; in accordance with the previous studies carried out in the SW of Iceland, e.g. Keiding et al., 2010;Juncu et al., 2017). We considered several geometries for the deformation sources. The observed surface deformation pattern does not fit well to a vertical dike; therefore we considered other geometries: point sources of pressure (Mogi, 1958, see Lisowski (2007 for formulation in terms of volume change), prolate spheroids (alike Yang et al. (1988) following the formulation of the MATLAB package dMODELS; Battaglia et al., 2013), or horizontally aligned penny shaped cracked (similar to Fialko et al. (2001) following the formulation of the MATLAB package dMODELS; Battaglia et al., 2013).
We jointly invert all geodetic data (i.e. four maps of LOS SAR displacements and horizontal displacements from 62 GNSS stations covering the 5 months uplift deformation of 2017-2018; Geodetic Data details the data sets used) to estimate the optimal parameters for two sources: a deflating source in Húsmúli and an inflating source in Ölkelduháls. To limit the computational time of the inversion, we use a variancebased quadtree algorithm from Jónsson et al. (2002) to subsample our InSAR data set. We opt for an increased number of InSAR data points in areas of steeper gradient in the deformation field, between the Ölkelduháls and Húsmúli sources, with the aim to optimize the fitting of our models to the observed complex deformation caused by two adjoining sources. The algorithm subsamples adequately the steep gradient between the main broad inflation signal and local subsidence signal in Húsmúli; however, we note that data points are also increased in areas of higher noise (e.g. topography-related noise in eastern Hengill). As this noise is different between the four tracks used for inversion (i.e. T155, T16, T118 and T41), we minimize any potential biases introduced in our inversion. To quantify the limits of spatial correlation and subsequent weight of our individual SAR tracks, we use the bootstrapping approach of Bekaert et al. (2015) to estimate the nugget, range and sill parameters, which are then used to form the variance-covariance matrix.
The posterior probability density function (PDF) is then formulated using the residuals of InSAR and GNSS (difference between the predicted and observed displacements) and their respective weights. We use the cascading adaptive transitional metropolis in parallel (CATMIP) algorithm from Minson et al. (2013) to sample the parameter space. The CATMIP algorithm is similar to simulated annealing (Kirkpatrick et al., 1983), iterating through intermediate PDFs towards the final posterior PDF, spawning additional Markov chains in the process for better coverage of the parameter space. To explore the location of the Húsmúli and Ölkelduháls sources we initially use two point sources (Lisowski, 2007). We further test different geometries for each source: a prolate spheroid (Battaglia et al., 2013) to potentially explain the elongated deformation signal observed in Húsmúli (Figure 4), and a penny-shaped crack to explain the broad scale uplift in Ölkelduháls (Battaglia et al., 2013). These Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 725109 9 simple mathematical formulations of the sources are valid regardless of their magmatic or hydrothermal origin. The three combinations of sources tested are: two point sources, a prolate spheroid and a point source, a prolate spheroid and a penny-shaped crack.
For simplification, we fix the semi-minor axis length of the prolate spheroid to 50 m, as initial freer inversions showed that this parameter gravitates towards small values. This particularly helps constrain the aspect ratio (a; ratio between the values of semiminor axis and semi-major axis) and volume change related to the prolate spheroid source, calculated in the same manner as Battaglia et al. (2013). We estimate the preferred model parameters as mean of the sampled parameters distribution. Due to the large number of degrees of freedom in our inversions (∼4 × 10 6 ), a consequence of the number of data sets jointly inverted and increased subsampling in eastern Hengill, we compute the reduced Chi square (χ 2 ] ) and degree of freedom on a subsampled subset of our InSAR residuals containing the deformation signal. These statistics are then used to quantify the goodness of fit of the model combinations tested.

RESULTS
The InSAR and GNSS times series analysis (Geodetic Data) highlights two main ground deformation signals in the Hengill area in 2017-2018 that differ from the general trend of the prior years ( Figure 6). These two adjacent ground deformation signals are particularly evident in the near-up S1 InSAR displacement map and horizontal GNSS displacements plotted in Figure 4. A broad scale uplift with a diameter of ∼11 km, with up to 12 mm of vertical displacement between the Hrómundartindur and Hengill central volcanoes, close to the Ölkelduháls cGNSS, is the largest signal ( Figure 4). A more localized area of subsidence northwest of the Hellisheiði geothermal field, Húsmúli, with a diameter of 6 km and maximum vertical displacement of -8 mm is evident west of the main uplift signal. From our analysis of the time-span of uplift (Uplift time span) we estimate it started in early November 2017 and lasted until the end of March 2018, spanning approximately 5 months. Table 2 summarises the optimum parameters we obtain from inversion of the geodetic data using three combinations of models 2 | Best-fit solutions of the joint inversion of GNSS and InSAR displacements for the Ölkelduháls and Húsmúli sources. Results of three models are presented here: two point sources, a prolate spheroid and a point source, and finally a penny-shaped crack and prolate spheroid source. a signifies the aspect ratio of the prolate spheroid source. In brackets are the 95% confidence intervals for each parameter. The reduced Chi square (χ 2 ] ), weighted residuals sum squared (WRSS) and degrees of freedom presented here were calculated using GNSS and subsets of the subsampled InSAR displacements. The preferred model ensuing our statistical test is highlighted in grey.   Table 2) indicates that a contracting prolate spheroid (Húsmúli source) and an expanding point source in Ölkelduháls is the simplest model combination explaining satisfactorily the changes in deformation observed in the Hengill area between 2017 and 2018. Figure 7 shows the data, subsampling, model-predicted and residual displacements according to the best fit parameters ( Table 2). Results from the other model combinations tested, as well as preliminary models, using the time-differenced time series (assuming the subsidence was ongoing through-out 2017-2018, as described in Isolating and modelling the 2017-2018 uplift signal) are provided in Supplementary Material. The subsampled InSAR data have an increased number of points by the Ölfus river ( Figure 1) and at the topographic issues in the eastern part of Hengill for most of the tracks. We additionally considered and tested different masked options of data and even-sampling of the data sets, to test the robustness of our results. This led to similar results. We opted for finer subsampling as it better samples the high gradient of the deformation between the Húsmúli subsidence and Ölkelduháls uplift ground motions.

Source/Parameters
The source of inflation in Ölkelduháls appears to be deep-seated (∼6 km depth) with a volume increase of ∼3.1 × 10 6 m 3 if we consider a point source model. Some GNSS and InSAR residuals (Figure 7) are observed in the Hellisheiði, Hverahlíð and Nesjavellir geothermal fields, which may hint to localized changes in deformation signals within 2017-2018. These changes could be linked to harnessing of high-enthalpy fluids in these localities.
The deformation signal in Húsmúli is caused by the superposition of transient and localized uplift and subsidence episodes FIGURE 7 | Data, subsampled InSAR data, model-predicted and full residuals displacements of GNSS and InSAR in the LOS of the tracks T41 (TSX), T16 (S1), T118 (S1) and T155 (S1) according to our favoured model ( Table 2). The cross presents the location of best fit point source (following the expressions of Lisowski, 2007) for Ölkelduháls (black) and the elongated ellipse for the prolate spheroid source (Battaglia et al., 2013) for Húsmúli (white) sources.

Insights From Uplift and Subsidence Episodes in Other Volcanic Systems
Alternating subsidence and uplift episodes have been observed in numerous silicic and basaltic volcanoes around the world (e.g. Campi Flegrei, Alutu, Yellowstone, White Island, Krýsuvík, Krafla, Theistareykir Waite and Smith, 2002;Biggs et al., 2011;Cannatelli et al., 2020;Fournier and Chardot, 2012;Michalczewska et al., 2012;Metzger and Jónsson, 2014;Gudjónsdóttir et al., 2020). The nature of non-eruptive inflation and deflation sources in volcanic systems is often still debated. For example, decades of studies consider possible hydrothermal and/or magmatic processes as the explanation of meter-scale uplifts and subsidence episodes as well as mini-uplifts observed at Campi Flegrei (Italy). The studies particularly highlight the importance of integrative modelling and investigations to further understand hydrothermal and magmatic interactions at volcanoes (e.g. Gaeta et al., 2003;Lima et al., 2009;Troiano et al., 2011;Chiodini et al., 2015;Cannatelli et al., 2020).
Possible documented processes at origins of succeeding uplift and subsidence episodes in magmatic areas include: solidification of an existing magmatic body and new intrusive material, motion of hydrothermal fluids, as well as trapping of magmatic or hydrothermal fluids and subsequent phase changes (e.g. Lima et al., 2009;Chang et al., 2010;Biggs et al., 2011). Notable clues associated with migration of fluids or magmatic material from deeper to shallower depth in the crust can be observed via changes in: subsidence rates leading to the subsequent uplift, seismic rates and/or shallow geothermal reservoir conditions (e.g. Waite and Smith, 2002;Lu et al., 2003;Lima et al., 2009). In the following section we discuss the various geodetic, seismic and geothermal borehole information available during the 2017-2018 Ölkelduháls uplift.

Ölkelduháls
The 2017-2018 source of inflation and the 2006-2017 contracting source in eastern Hengill (Results and Table 2; Juncu et al., 2017), appear to be in proximity (∼3 km NNW) of the point source estimated for the uplift episode in 1993-1999, at a depth of approximately 5-7 km (Figure 4; Sigmundsson et al., 1997;Feigl et al., 2000). This depth range is near or within the inferred brittle-ductile transition zone of the area (Foulger, 1995;Li et al., 2019). The inferred 2017-2018 inflating and 2006-2017 contracting sources are located under Ölkelduháls, between the Hengill, Hrómundartindur and Grensdalur volcanic systems (Figure 4 and Figure 7). Similar cases of inflating sources located between central volcanoes have been observed, for example, in the Three Sisters volcanic center (United States) (Dzurisin et al., 2009;Ebmeier et al., 2018). The point source model used can be a proxy for sources of both hydrothermal and magmatic nature. We also tested a pennyshaped crack model, often used to model magmatic sources, to explore other geometries (Section 3). Although our statistical test does not favour this model over a simple point source model, it is not possible to disregard the results from the penny-shaped crack model as it explains the broad uplift signal well (see Supplementary Material for the model-predicted and residual displacements). The optimum depth associated with the pennyshaped crack model is ∼11 km, deeper than the point source model (∼6-7 km). This estimated depth is shallower than the inferred Mohorovičić discontinuity (∼14-15 km; Weir et al., 2001) but deeper than the inferred brittle-ductile boundary in the area (Foulger, 1995;Li et al., 2019), indicating a possible lower-crustal deformation source.
The interpretations from our modelling presented here is limited by the assumption of static and immediate response of the homogeneous medium to a volume change of a source (i.e. elastic half-space model). Post-seismic transients and timedependent deformation signals have been studied in the SW of Iceland (e.g. poroelastic or viscoelastic behaviour; Jónsson et al., 2003;Decriem and Árnadóttir, 2012). The temporal resolution of the deformation signal is limited to a few continuous GPS stations in the Hengill area ( Figure 3; mainly OLKE, HVER and NVEL). Moreover, the inflation episode appears short-lived. Hence, we do not attempt to isolate nor model transient deformation, such as due to a poroelastic and/or viscoelastic processes. More data to improve the spatial and temporal resolution, and longer time series are necessary to study such transient signals. While more complex rheological models and transient signals should be considered in further modelling of the brittle-ductile boundary sources of uplift and subsidence in the area, our main objective here, is to document a short-lived uplift episode during 2017-2018 in the Hengill area, and provide a simple model to explain the observations. The physical properties of the source(s) in Ölkelduháls, possibly causing the extended period of subsidence and a brief uplift episode, are unknown. Do the signals stem from a single source volume, or are two separate volumes acting at the same time, or possibly a layered source? A layered source, here, can be envisioned as vertically stacked sources, with a possible conduit connecting the sources. However, due to the low signal-to-noise ratio of our data, we have studied simple analytical models and do not explore more complicated sources or interactions of sources which can be investigated using higher order numerical parametric models (e.g. Chang et al., 2010;Walwer et al., 2019;Delgado and Grandin, 2020).
Magmatic intrusions can sometimes be more irregular in shape than the simple geometries we can apply analytically (e.g. sill, dike, or spheroid), thus the deformation might not be fully explained by the simple mathematical expressions tested in this study. Migrations of fluids towards the surface (e.g. of magmatic or hydrothermal nature), can be accompanied by decaying subsiding motions and seismic events (Waite and Smith, 2002;Lu et al., 2003). Integrative modelling relating Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 725109 ground deformation and hydrothermal processes, particularly the possible poro-elastic or thermal expansion of such fluids, could possibly explain the rapid cycling between uplift and subsidence near Ölkelduháls (e.g. Campi Flegrei (Italy), Whaakari/White Island (New Zealand); Gaeta et al., 2003;Troiano et al., 2011;Chiodini et al., 2015;Fournier and Chardot, 2012). In the following two subsections (Hydrothermal Information and Seismicity Rates), we discuss additional geothermal and geophysical information available in the Ölkelduháls area prior to and during the 2017-2018 and consider the implications of these data.

Hydrothermal Information
Pressure and temperature measured in geothermal systems vary with time under host-rock temperature fluctuations, fluid motions, and may change via natural (e.g. cooling, heating, percolation of meteoric water) or anthropogenic processes (e.g. injection, extraction of fluids). A temporal analysis of the pressure and temperature at depth measured in three wells (HE-02;800 m, HE-20;1200 m and HE-22;500 m) in Ölkelduháls is presented in Figure 8 (borehole locations in Figure 1) using monitoring data sets of Reykjavík Energy group (see 2015-2018 pressure and temperature profile in Supplementary Material). No anomalous changes in pressure or temperature are observed during 2017-2018 compared to the preceding years. All three boreholes seem to show a gradual decrease in pressure between 2011 and 2019. Data from wells HE-20 and HE-22 show a slight decrease in temperature both within that same time span. Measurement from well HE-22 should, however, be taken with caution, as the shallow depth of our available measurements cannot ensure that the temperature and pressure measured represent the geothermal reservoir conditions. Overall, these measurements do not indicate fluctuations in temperature and pressure at shallow depth which could be resulting from deeper upwelling volatiles or rise of temperature from the shallow geothermal systems.
A seismic tomography study from Tryggvason et al. (2002) carried out after the 1993-1999 inflation beneath Hrómundartindur suggests that the footprint of seismic velocities between 5 and 7 km depth is more consistent with seismic velocities of rocks which pore space is filled with supercritical fluids, contradicting earlier interpretations of a large melt body at those depths as the origin for the 6-year uplift episode. Small bodies of melt could, however, potentially be present. Microseismicity (depth range around 3 km), resistivity and tomographic studies were carried out recently in the Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 725109 eastern part of Hengill that focused on shallower depths (Árnason et al., 2010;Jousset et al., 2011). These studies supported the previous tomographic results. Fluid motions are possible within the brittle-ductile transition zone of the Icelandic oceanic crusts (Violay et al., 2012). If one assumes that the physical conditions of the fluids may be akin to a supercritical state of pure water, then one can speculate on the potential effect induced by the motions of such fluids. An isothermal (around 650 ± 100°C; Foulger, 1995) decompression of fluids from 7 km (an equivalent of ∼220 MPa of lithostatic pressure) to 6 km (∼175 MPa) would decrease the specific density of the fluid from ∼0.6 g cm −3 to ∼0.3 g cm −3 (Anisimov et al., 2004;Imre et al., 2012). The volume being the inverse of the specific density, one can infer that the decompression of the fluids between 7 and 6 km would approximately double their initial volume. Future pressure and temperature measurements in the boreholes in the Ölkelduháls area may potentially show delayed changes in the shallow hydrothermal reservoir, in response to 2017-2018 changes in the deeper part of the geothermal system in the eastern Hengill region.

Seismicity Rates
Several studies report intense seismic activity in the Hengill area, from prior deep-seated inflation episode in 1993-1999 (tens of thousands earthquakes; Sigmundsson et al., 1997;Feigl et al., 2000) to anthropogenic activity in the Hellisheiði geothermal power plant in September 2011 (Juncu et al., 2020). The estimated amplitude of the 2017-2018 uplift is a ninth of the amplitude of the 1993-1999 uplift; however, the rate of uplift during 2017-2018 (lasting 5 months) is approximately 2.9 cm/yr, similar and slightly above the uplift rate estimated during 1993-1999Feigl et al., 2000). We use seismicity as an indicator of possible stress changes during the transient 2017-2018 uplift compared to the background stress of the area (e.g. 2006-2017 deformation and plate motion). We consider seismic records of the Icelandic Meteorological Office (IMO) of earthquakes between 1 January 2015-30 June 2019 located near Ölkelduháls. Clusters of seismic activity associated with geothermal extraction in Hellisheiði and Nesjavellir are excluded from our data set. To minimize the effect of possible increase or decrease in the number of operational seismic stations by the IMO network during January 1, 2015 until the June 30, 2019, we exclude from our analysis earthquakes with a local magnitude (M L ) below 1. The seismic moment is calculated assuming that local magnitude of earthquakes below three is equal to the moment magnitude (M W ) (after the findings of Greenfield et al., 2020) according to the formulation of Kanamori (Kanamori, 1977, Figure 9). The largest event in the data set is a M L 2.9 occurring on September 25, 2016 at ∼5 km depth. A small number of seismic events (M L < 2) are located within 4-6 km depth during the uplift time span, proximal to the depth of our preferred inflating source, Table 2). However, earthquakes occurring at this depth range appear to be common and no clear seismic swarms resembling those recorded during the 1993-1999 uplift episode (Sigmundsson et al., 1997) can be observed. The lack of earthquakes below ∼6 km confirms the depth of brittle-ductile transition zone in the area, as found in previous seismic studies (e.g. Foulger, 1988;Tryggvason et al., 2002;Li et al., 2019).
A physical explanation for the 2017-2018 uplift in the Hengill area could be a new intrusion at depth, related or unrelated to the prior historical intrusions in the area, in a zone of relative weakness beneath the Hengill-Hrómundartindur-Grensdalur volcanic complex. Magmatic intrusions are often associated with increased seismicity, induced by stress changes (McNutt and Roman, 2015). However, we do not observe any marked increase nor decrease in seismicity during the time span of the 2017-2018 uplift episode (Figure 9). A possible explanation for the lack of recorded seismicity may be the proximal nature of the source to the brittle-ductile boundary; in which case, triggered events may be below the magnitude threshold of our analysis. Secondly, due to the close location of the inflating-deflating source(s), the 2017-2018 uplift is only a 5 months reversal of 11 years accumulated stresses from the 2006-2017 subsidence, and thus not enough to trigger significant seismic events. Magmatic or hydrothermal fluids may cause changes in pore pressure, which in turn could change seismicity rates. Moreover, the lack of change in seismicity we observe could be because the fluids have not migrated to active faults.

Origin of the Ölkelduháls Source
The cause(s) of the 2017-2018 inflation and 2006-2017 deflation episodes in the Ölkelduháls area, eastern Hengill, remain unclear. Intrusion, solidification of a magmatic body, motions of fluids, or more complex processes such as upwelling, accumulation and release of trapped magmatic or hydrothermal fluids to shallower depth as suggested by Lima et al. (2009) may be considered, although seismic and borehole information (Figure 8 and Figure 9) do not highlight particular changes during the time span of the uplift compared to prior years. Altogether, we cannot exclude that the 2006-2017 subsidence and 2017-2018 uplift episodes in Hengill are independently related to both magmatic and hydrothermal processes. Further integrative studies, including gravimetric, fluid chemistry or time-lapse resistivity measurements may help improve our understanding of such episodes as well as understanding the state of unrest of the Hengill-Hrómundartindur-Grensdalur volcanic systems (e.g. Tizzani et al., 2015;Pritchard et al., 2019). Nevertheless, the ground deformation documented here confirms temporal changes in the activity of the Hengill-Hrómundartindur-Grensdalur volcanic systems and the need to closely monitor the area for the foreseeable future.

Origin of the Húsmúli Source
The parameters (e.g. volume change, depth) resulting from the inversion of the Húsmuli source highlight a potential relationship of these transient episodes and the known structural tectonics of the area: the elongation in ∼N30°; of the prolate spheroid source corresponds to the trend of the fissure swarms, large fracture zones and faults of the Hengill area (Saemundsson, 1967;Clifton et al., 2002;Steigerwald et al., 2020). Adjacent transient uplift and subsidence ground deformation signals have been noted in other volcanic systems around the world and related to an event of magmatic or hydrothermal motions (e.g. Fernandina (Galapagos), Akutan, Yellowstone (United States); Pepe et al., 2017;Wang et al., 2018;Delgado and Grandin, 2020). The location and shallow depth (<3 km) of the Húsmúli contracting source coincides with the depth of the main extraction and injection boreholes of the Hellisheiði power plant and within similar crustal volume as man-made hydrothermal deformation sources documented in Juncu et al. (2017Juncu et al. ( , 2020. We interpret that the changes observed in 2017-2018 are likely related to localized hydrothermal fluid extraction and injection processes.

CONCLUSION
The Hengill area is the locus of documented deep-seated sources of inflation (6-7 km depth; 1993-1999, 5 months in 2017-2018) and deflation (2006-2017; 2018-2020). Using GNSS, InSAR and analytical models, we conclude that the 2006-2017 subsidence and following (2017-2018) transient uplift in the area seem to originate from a similar location in the crust, near the brittleductile transition zone of the area. The nature of the processes at the origin of these motions is inferred to be hydrothermal or magmatic in nature. Localized ground motions in Húsmúli between 2015 and 2016 (uplift) and 2017-2018 (subsidence) are the responses of shallow (<3 km) processes related to geothermal extraction and injection of the Hellisheiði geothermal plant. The 2017-2018 surface ground motions highlight the intricate connection between geothermal exploitation, crustal deformation and tectonic deformation in the area. Further geophysical and geochemical studies are needed to understand the origin of these deep-seated inflation-deflation episodes in the Hengill-Hrómundartindur-Grensdalur volcanic systems. This study on deep seated sources in the Hengill area demonstrates the need to better understand the driving processes of long-lived high temperature geothermal fields and behaviour of volcanoes.

DATA AVAILABILITY STATEMENT
Initial InSAR mean LOS velocities grids and GNSS displacements are detailed in Supplementary Material and are available on request to the authors. Original Sentinel-1 images are available on the Copernicus and ESA platform. Seismic data sets used in the study can be requested from the Icelandic Meteorological Office or downloaded from https://skjalftalisa.vedur.is. Borehole data sets used in the study can be requested from Reykjavík Energy.

AUTHOR CONTRIBUTIONS
CD, HG, TA designed the study. CD made the figures and wrote the paper. DJ and TA provided the original inversion code, then adapted by CD. VD provided codes and interferograms for Sentinel-1 time series. HG, SH and TA provided GNSS data sets. CD, DJ and VD created the InSAR time series. SH assisted CD with the GNSS processing and processed the data used in Figure 2 and Supplemental Material. BK, ST, and GG provided the borehole pressure, temperature and traces. HB provided the seismic data. All authors contributed via discussions, comments and edits on the drafts.

FUNDING
This work was funded via the Icelandic Research Fund, grant 174377-015 and NordVulk fellowship programme 2020.