Environmental And Anthropogenic Gravity Contributions At The Þeistareykir Geothermal Field, North Iceland

doi:10.21203/rs.3.rs-692087/v1

Continuous high resolution gravimetry is increasingly used to monitor mass distribution changes in volcanic, hydrothermal or other complex geosystems. To quantify the often small target signals, gravity contributions from, e.g., atmospheric mass changes, global and local hydrology should be accounted for. We set up three iGrav superconducting gravity meters for continuous monitoring of the Þeistareykir geothermal field in North Island. Additionally, we installed a set of hydrometeorological sensors at each station for continuous observation of local pressure changes, soil moisture, snow and vertical surface displacement. We show that the contribution of these environmental parameters to the gravity signal does not exceed 10 µGal (1 µGal = 10^-8 m s^-2), mainly resulting from vertical displacement and snow accumulation. The seasonal gravity contributions (global atmosphere, local and global hydrology) are in the order of ±2 µGal at each station. Using the environmental observations together with standard gravity corrections for instrumental drift and tidal effects, we comprehensively reduced the iGrav time series. The gravity residuals were compared to groundwater level changes and geothermal mass flow rates (extraction and injection) of the Þeistareykir power plant. The direct response of the groundwater levels and a time-delayed response of the gravity signal to changes in extraction and injection suggest that the geothermal system is subject to a partially confined aquifer. Our observations indicate that a sustainable “equilibrium” state of the reservoir is reached at extraction flow rates below 240 kg s^-1 and injection flow rates below 160 kg s^-1.

Physical Geography

Planetary Geology

Geothermal monitoring

superconducting gravimetry

time series analysis

gravity reduction

Þeistareykir

Many geoscientific studies make use of continuous terrestrial gravity measurements to estimate subsurface mass changes associated with, for example, groundwater flow (Watlet et al. 2020), volcanic activity (Jousset et al. 2000; Carbone et al. 2019) or exploitation of geothermal reservoirs (Portier et al. 2018). Additional gravity effects from surrounding environmental processes, like Earth tides, atmospheric pressure changes, rain and snowfall superimpose onto the target signals of the gravity time series and have not been always considered, making interpretation inaccurate when gravity variations are small. The amplitudes of the environmental contribution in the recorded gravity signals range from a few nGal to several hundreds of µGal and vary from seconds to years (Boy and Hinderer 2006; Damiani 2014). Mikolaj et al. (2019) pointed out the need to reduce gravity contributions from atmospheric mass displacement, large-scale hydrology and nontidal ocean loading for local geophysical applications of terrestrial gravimetry. They quantified the uncertainties of these corrections to be in the order of several nm s^− 2, depending on the locations and the time scale of the mass variations of interest. Gitlein et al. (2013) modelled the combined gravity contributions of local and global atmospheric mass changes and applied them for reduction of superconducting gravity data, which improved the residuals by about 15% compared to the standard air pressure reduction with an admittance of 3 nm s^− 2 hPa^− 1. Voigt et al. (2021) used a superconducting gravity meter for hydrological monitoring of Mount Zugspitze and identified snowpack as primary contributor to seasonal water storage variations, with a snow-gravimetric footprint of up to 4 km distance around the gravity meter. In addition to environmental gravity contributions, artefacts of the gravity meter like instrumental drift or self-noise may decrease the accuracy of the target signal (Crossley et al. 2004; Rosat and Hinderer 2018). With respect to long-term drift rates, Schäfer et al. (2020) showed the impact of transporting a superconducting gravimeter. Therefore, an accurate estimation and reduction of instrumental drift and environmental contributions is essential before we can interpret the gravity residuals accurately with respect to a specific geophysical phenomenon.

With the aim of estimating subsurface mass changes associated with geothermal exploitation (geothermal fluid extraction and reinjection), we set up a network of three continuously recording and remotely operated iGrav superconducting gravity meters, together with GNSS stations and hydrometeorological sensors (humidity sensors, snow gauges, pressure sensors, temperature sensors, etc.) at Þeistareykir, a geothermal field in North Iceland. This region is part of the North Volcanic Zone and is located close to the Husavik fault (Gudmundsson et al. 1993) and the Krafla Caldera (Ármannsson et al. 1987). Surface deformation and seismic activity associated to the mid-oceanic ridge have been monitored in Iceland for over 50 years (Sturkell et al. 2006). Surface explorations at Þeistareykir in the 1970s and 1980s suggested beneficial reservoir temperatures (< 280°C), which was confirmed after exploration well drilling between 2002 and 2012 (Óskarsson 2015). We started our observations in December 2017 shortly after a new geothermal power plant (with 90 MWe total capacity) started operation in Þeistareykir. Figure 1 shows the location of the three gravity stations. The east station (iGrav032) is positioned in the vicinity of the geothermal extraction wells and the west station (iGrav006) at the reinjection wells. In June 2019, we started continuous measurements with iGrav015 at the central station, in the transient area between geothermal extraction and injection. Schäfer et al. (2020) provide further details about the gravity station setup and the continuous observation chronology. Gravity observations are subject to various environmental contributions (e.g. rain, snow, wind and soil humidity). Therefore, we deployed a series of instruments for measuring environmental parameters, which contribute to the total recorded gravity signal. Over more than three years, we acquired a unique data set of high-resolution gravity and environmental time series. We developed and applied an approach to reduce all environmental effects that may hide the geothermal mass changes.

In this work, we assess the contributions of different geophysical phenomena to the observed gravity signals, based on observations and models. Eq. 1 shows the gravity contributions examined in our study.

These include gravitational effects induced by solid Earth and ocean tides Δg_tide, local atmospheric mass changes Δg_atm, polar motion Δg_pm, global atmosphere, large-scale hydrology and nontidal ocean loading Δg_glob, local hydrology in terms of soil moisture changes Δg_hydr, snow Δg_snow, vertical surface displacement Δh with the vertical gravity gradient δg/δh, geothermal mass changes Δg_geoth and instrumental drift Δg_drift. The “Error” in equation 1 consists of further gravity contributions like magma-induced mass changes from volcanic activity that are neglected in this study. Before the gravity data can be used for interpreting geothermal related mass changes, we have to identify and reduce the interfering gravity effects that superimpose onto the target signal. To assess the contribution of the local environmental parameters in the gravity records, we analysed the continuous hydrometeorological measurements from remotely operated multiparameter stations (ROMPS; Schöne et al. 2013) at each site. Figure A1 in the Appendix shows top view sketches of the ROMPS sensors and gravity containers for the three gravity stations. In the following subsections, we describe details of the instruments, the methods and the models for these individual gravity contributions, and we show the observed environmental signals from Þeistareykir.

Earth tides (Δg_tide), local pressure effects (Δg_atm), polar motion (Δg_pm) and instrumental drift (Δg_drift)

The largest gravity signal results from the solid Earth and ocean tides and is estimated by tidal modelling (Agnew 2015). We calculated the barometric admittance factors for each iGrav to correct the gravity effects of local pressure changes (Merriam 1992). Both, local tidal models and barometric admittance factors were computed for each gravity station at Þeistareykir using the ETERNA 3.4 package (Wenzel 1996). Results of the tidal analyses are given in the Appendix (chap. A1, A2 and A3). The polar motion effect is provided as Earth orientation parameter file (EOP) by the International Earth Rotation and References Systems Service (IERS; https://hpiers.obspm.fr/eoppc/eop/eopc04/, Accessed 02 July 2021). We used the TSoft package (Van Camp and Vauterin 2005) for computation of the associated gravity contributions. Next, we examined and corrected the individual drift behaviour of the instruments by comparison to absolute gravity measurements (Hinderer et al. 2015). For this, we performed absolute gravity campaigns (with FG5#206) at each station once every year and adjusted the continuous time series to the absolute values. Detailed explanation of these standard corrections is given in Schäfer et al. (2020).

Global gravity contributions (Δg_glob)

In order to determine the contribution of global atmospheric mass variations in the gravity signal, we used a local air pressure measurements and the Atmacs model (Klügel and Wziontek 2009). We combined the results of local air pressure measurements at each gravity station with the Atmacs pressure model for Þeistareykir: We applied the remove-restore technique suggested by the Atmacs service (http://atmacs.bkg.bund.de/docs/computation.php, Accessed 02 July 2021) by removing the model surface air pressure from the total atmospheric effect (both calculated by Atmacs) and replacing it with the local pressure data from the gravity stations. For the computation of large-scale hydrological effects, we applied the land surface model NOAHv21 of the GLDAS model (Rodell et al. 2004). In addition, we computed nontidal ocean loading with the OMCTv06 model (Dobslaw et al. 2017). The simulated storage and mass variations from both models were converted to gravity effects using the mGlobe toolbox (Mikolaj et al. 2016).

Local hydrology (Δg_hydr)

At each station we recorded the variations of soil water content with in-situ soil moisture sensors. The sensors are arranged at different depths within soil profiles and at different distances to the gravity meter pillar (Table 1). From the soil moisture time series of all three gravity stations, we calculated the mean water content variations and their associated standard deviation at the four measurement depths (10 cm, 30 cm, 50 cm, 80 cm) and assigned those to the respective soil layers (0–20 cm, 20–40 cm, 40–65 cm, 65–95 cm). Assuming that the temporal variability of soil moisture decreases linearly with depth, we used the observed depth-dependence of the standard deviation to extrapolate at which depth it becomes zero, i.e., the threshold at which depth temporal variations of soil water content can be expected to vanish (Fig. A2, Appendix). The soil moisture variations of the deepest observation depth at 80 cm were accordingly extrapolated to this threshold depth (here 1.8 m, see Fig. A2, Appendix).

The calculated soil water content variations at the different depths, expressed in millimetre water equivalent were used as input variable to model the local hydrological gravity effects. Further, we included local digital elevation models (DEM) to account for topographic characteristics (i.e. relative height changes of the soil layers with regard to the gravity sensor). The hydro-gravitational modelling (HyGra) is based on the method of Leirião et al. (2009) and was adapted to SG observatory buildings by Reich et al. (2019). HyGra is a spatially distributed model that enables the setup of a nested component grid with adjustable radii around the gravity station. We chose a small lateral discretisation (of 0.1 m) for the model cells closest to the gravity sensor (radius of 50 m) and larger model cells with increasing distance (1 m cells for 50 to 300 m radius and 10 m cells for 300 to 3000 m radius), and a vertical discretization (cell height) of 0.1 m for every cell. The gravity sensor height above the ground surface of the DEM is 1.0 m. For the volume of the field container and the sub-surface column below the container, no (hydrological) mass changes were assumed because of the umbrella effect of the container.

Table 1

Distribution of soil moisture sensors installed at different depths around each gravity station; some soil profiles are deployed at equal distances to the iGrav pillar (e.g. two profiles with 12.1 m distance at the west station), in these cases the sensors are installed below different micro-topographic features (e.g. small hills or terrain depressions).
Gravity station	Nr. of sensors	Nr. of profiles	Profile distance to iGrav pillar [m]	Sensor depths [m]
West	14	4	6.5, 7.5, 12.1, 12.1	0.1 (4x), 0.3 (4x), 0.5 (4x), 0.8 (2x)
Central	13	5	10.3, 11.9, 11.9, 13.4, 13.4	0.1 (5x), 0.3 (4x), 0.5 (3x), 0.8 (1x)
East	13	4	8.1, 8.4, 11.4, 11.4	0.1 (4x), 0.3 (4x), 0.5 (3x), 0.6 (1x), 0.8 (1x)

Snow (Δg_snow)

To determine the mass changes of the snow cover around the monitoring stations in the course of snowfall and snowmelt, we continuously measured (every 15 minutes) the snow water equivalent (SWE), i.e., the amount of water that is stored in the snow cover in solid and liquid state. As independent observations of SWE, we installed two snow monitoring instruments at the eastern gravity station (Fig. 2). We used a snow scale to determine the SWE by weighing the column of snow that is on top of a pressure pillow of 6.72 m² in size. Additionally, we used a snow pack analyser system (Sommer Messtechnik) equipped with strap sensors that measure (with an electro-magnetic approach) the specific volume contents of ice, water and air within the snow cover, from which the SWE is calculated in combination with a snow height measurement of an ultra-sonic device. We used the mean SWE of the time series of both measuring systems (snow scale and snow pack analyser) as input to calculate the gravity effect of snow with the HyGra model. It should be noted that the calculation of the snow gravity effect considered the actual thickness of the snow cover relative to the gravity meter sensor height, so that both positive and negative gravity contributions may occur, depending on whether parts of the snow cover are below or above the sensor, respectively. Piling up of snow on top of the gravity container was minimal because of the windy conditions at the sites, as confirmed by daily photos of automatic cameras (see Fig. A1, Appendix) deployed at each station. Also, based on the camera observations, no snow mass accumulation was observed within the first two meters around the container due to snow drift by wind. Thus the snow mass in the near field of the gravity meter was neglected when computing the gravity effect of the snow cover.

Vertical surface displacement (Δh)

We used GNSS data observed at each gravity station at Þeistareykir to estimate the vertical surface displacement. The GNSS processing was performed using GFZ’s EPOS.P8 software based on a classical network approach while introducing satellite orbits, clock corrections, and Earth rotation parameters from GFZ repro3 solution (Männel et al. 2020, 2021). We processed undifferenced GPS observations in the ionosphere-free linear combination and fixed the ambiguities to their integer values. Finally we estimated daily coordinates, Earth rotation parameters, hourly troposphere delays using a piece-wise linear function, and daily troposphere gradients. The geodetic datum was realised by applying no-net-translation and rotation conditions to the core stations of the IGS14 frame (Johnston et al. 2017). According to the current IERS Conventions 2010 (Petit and Luzum 2010) nontidal surface loading was not corrected.

To account for the contribution of the observed height changes to gravity we used the vertical gravity gradient. As a first approximation, the free-air vertical gravity gradient (FAG) (Eq. 2) is given by:

with the Universal gravitational constant G (6.67 x 10^− 11 N m² kg^− 2), the mass M and radius R of the Earth (~ 5.98 x 10²⁴ kg and ~ 6.37 x 10⁶ m) (Brown and Rymer 1991). As we measured the FAG using a Scintrex CG5 gravity meter (Portier et al. 2020), we used 319 µGal m^− 1 for the west station (iGrav006), 330 µGal m^− 1 for the central station (iGrav015) and 307 µGal m^− 1 for the east station (iGrav032). However, Eq. 2 does not account for subsurface mass or density changes, which we expect to occur in an active geothermal reservoir (see Sect. 4).

Local environmental observations at Þeistarykir

Figure 3 displays the measured environmental parameters for the three years study period at the Þeistareykir site. Table 2 gives the minimum and maximum values, and standard deviations for each parameter. All time series show daily values, representing the smallest common time interval of the obtained data.

Figure 3a shows the soil water content variations at different depths averaged from all soil moisture measurements at the three Þeistareykir gravity stations. The short and long-term soil water content variations are similar in their overall dynamics in all depth layers, but their amplitudes decrease with depth. We observe largest soil water content variations at 10 cm depth with a standard deviation of 1.45 Vol% and decreasing variations with increasing depth (SD of 1.10 Vol% at 80 cm depth, see also Figure A2, Appendix).

The results of the snow measurements are shown in Fig. 3b. The SWE time series from the snow pack analyser (black line) and the snow scale (red line) show simultaneous responses to snow mass accumulation during the three winter periods. However, there is a large difference in amplitude between the two monitoring systems, with signals more than ten times larger for the snow pack analyser than for the snow scale. This can partly be explained, by the different positions of the two instruments at the east station (Fig. A3, Appendix). As a result of different wind exposure, the snow cover can be expected to be different at the two installation sites. This is also confirmed by time-lapse photos taken by the automatic cameras at the site. Snow conditions at the east station in the course of one year are shown in Figure A4 in the Appendix. A further possibility of the differences between the two time series are systematic measurement errors of the devices. The snow scale is known to underestimate SWE because of the internal cohesion of the snowpack, i.e., the formation of snow or ice bridges between the snow pack on top of the scale and the surrounding snow cover (e.g., Grossi et al. 2017). On the other hand, the lowest monitoring strap of the snow pack analyser used here for deriving SWE may tend to overestimate SWE because of the more compacted, denser snow pack at this depth, compared to the lower density snow in upper parts of the snow cover. In order to account for both measurements at equal weight, we used the mean SWE (blue line) from both time series for modelling the snow mass effect on the gravity monitoring.

In Fig. 3c we show vertical surface displacement observed from GNSS monitoring at the three gravity stations at Þeistareykir (blue, green and red dots). The periodical variations of approximately ± 5 mm per 3 months, clearly visible in the 90 days moving average, are caused by nontidal loading. For calculation of the vertical velocities, we applied a linear fit on the three years GNSS data. This revealed a significant subsidence of 11.1 mm yr^− 1 for the west station (blue line) and small trends of 1.5 and 0.3 mm yr^− 1 for the central and east stations (green and red lines). These values coincide with InSAR observations from the Iceland GeoSurvey (ÍSOR), which show increased subsidence rates of 7 to 8 mm yr^− 1 in the injection zone (west station) between summer 2018 and 2019 (Drouin et al. 2020). For the subsequent gravity reductions, we only used the contributions from the vertical velocities (linear trends) of the GNSS data.

Table 2

Quantities of environmental parameters (except tides) determined for the three gravity stations at Þeistareykir; including minima, maxima and standard deviations for soil water content variations at different soil depths, snow water equivalent and vertical surface displacement.
Environmental parameters		Unit	Min	Max	SD
Soil water content variation	10cm depth	Vol%	-4.98	3.41	1.45
	30cm depth		-4.00	3.89	1.26
	50cm depth		-3.12	3.63	1.19
	80cm depth		-3.33	3.52	1.10
Snow water equivalent	Snow pack analyser	mm	0.00	758.06	159.10
	Snow scale		0.00	91.55	16.33
	Mean snow water equivalent		0.00	408.13	86.02
Vertical surface displacement	West station	mm	-51.51	11.89	11.58
	West station 90 days average		-33.08	3.42	9.78
	Central station		-37.39	18.21	7.03
	Central station 90 days average		-14.83	2.75	4.17
	East station		-43.78	15.42	6.98
	East station 90 days average		-15.17	3.28	4.37

Environmental gravity reductions

Figure 4 shows the gravity contributions from global models (for the entire study area) and the observed environmental parameters (from Fig. 3) expressed in µGal for each of the three gravity stations at Þeistareykir. Table 3 gives the minimum, maximum and standard deviation for each parameter.

Among the global gravity contributions (Fig. 4a), the atmospheric effect (green line) shows the largest gravity contribution of more than 2 µGal. The gravity contributions modelled for global hydrology (GHE; red line) are dominated by a seasonal component, whereas the gravity signal of nontidal ocean loading (NTOL; blue line) is of higher frequency. Both components have a minor effect for the reduction of the local gravity residuals. We expect mean uncertainties of 0.38 µGal for GHE and 0.17 µGal for NTOL, based on the uncertainty assessment of Mikolaj et al. (2019). The total effect (black line) of these large-scale contributions to the local gravity observations has a seasonal amplitude of up to 4 µGal and a standard deviation of about 1 µGal (Fig. 4a, Table 3).

From all environmental contributions vertical surface displacement at the west station (11.1 mm yr^− 1; see also Fig. 3) causes the largest gravity effect of up to 10 µGal over the entire observation period (Fig. 4d). Global gravity contributions and variations in soil water content only exhibit a very small contribution to the gravity variations observed by all three iGravs. For the reduction of snow, the gravity effect for the east and central stations are larger than for the west station, visible between December 2019 and May 2020 (Fig. 4c). This can be explained by the different topography surrounding the three gravity stations, considered as input parameter for the snow effect modelling. At the west station with several little hills around it, a larger amount of snow cover is located at elevations higher than the gravity sensor of iGrav006. This creates a partly negative gravity effect and reduces the net gravity reduction of snow mass changes for the west station.

Table 3

Quantities of environmental parameters (in µGal) determined for the three gravity stations at Þeistareykir; including minima, maxima and standard deviations for global gravity contributions, soil water content variations, snow water equivalent and vertical surface displacement.
Environmental parameters	Station	Min	Max	SD
Global atmosphere	Entire area	-1.76	2.60	0.63
Global hydrology	Entire area	-0.71	1.49	0.55
Nontidal ocean loading	Entire area	-1.35	1.39	0.37
Total global effect	Entire area	-1.62	3.45	0.91
Soil water content variation	West	-1.60	2.05	0.58
	Central	-2.07	2.27	0.69
	East	-1.97	2.21	0.66
Snow water equivalent	West	0.00	1.87	0.47
	Central	0.00	11.07	1.85
	East	0.00	8.52	1.25
Vertical surface displacement	West	0.00	10.92	3.26
	Central	0.00	1.51	0.45
	East	0.00	0.28	0.08

Gravity residuals

In Fig. 5 we compare the time series of gravity residuals of iGrav006, iGrav015 and iGrav032, before and after the combined environmental reductions of global effects, soil water content, snow water equivalent and vertical surface displacement. The initial residuals (before the aforementioned environmental reductions were applied) have been derived by reducing Earth and ocean tides, polar motion, local pressure and instrumental drift, as well as removing spikes and offsets (caused by earthquakes and other disturbances) from the original gravity time series (Schäfer et al. 2020). Reductions of the individual environmental gravity contributions for each iGrav are shown in Figures A5 to A7 in the Appendix.

From linear regression of the reduced time series, depicted by dashed lines in Fig. 6, we determined a gravity increase of 3.2 µGal yr^− 1 for the west station (iGrav006) and a gravity decrease of -8.4 µGal yr^− 1 for the east station (iGrav032). Figure 6 also shows the gravity differences between the iGrav residuals after reduction of the environmental contributions. The gravity differences of two instruments exhibit an additional reduction of variability and amplitudes, i.e. they are smoother than the iGrav residuals themselves. This is because gravity changes that similarly appear at two stations cancel out to some extent in the differences. The remaining gravity variations can be attributed to effects that are locally different at the two stations and have not been reduced by the applied reduction models. Assuming that the local and large-scale gravity reductions as described above could reasonably be quantified based on the in-situ observations and the standard models, these remaining variations can be mainly attributed to geothermal mass changes (Δg_geoth in Eq. 1) by mass extraction in the vicinity of the east station (iGrav032) and reinjection of effluent water close to the west station (iGrav006). For iGrav015, the gravity signal is more similar to iGrav032 than to iGrav006, which is clearly visible in the lower residual amplitudes of the differences between iGrav032 and iGrav015 (red line) compared to the differences between iGrav006 and iGrav015 (orange line). This suggests that the central station receives a larger contribution to the gravity signal from mass extraction (depth ~ 2 km) rather than from mass reinjection (depth ~ 400m).

One particular phenomenon is the large gravity increase of about 10 µGal between May and June 2020, which is most pronounced for iGrav006. We assume this is due to local accumulation of snowmelt water (see also higher values of soil water content for that time; Fig. 4b) which drains off more easily at the other stations. At the west station there is a small nearby lake which may contribute to this pronounced water storage increase for a period of several weeks.

Geothermal signals

We investigate the subsurface mass changes that are primarily caused by geothermal fluid extraction and reinjection activities at the Þeistareykir power plant. For this purpose, we compare the gravity residual time series (from Fig. 6) with relative groundwater levels (GWL) and geothermal flow rates at Þeistareykir. The latter two we obtained as local monitoring data from the power plant operator (Landsvirkjun). The locations of the GWL monitoring wells and the injection well pad are depicted in Fig. 1. Figure 7 shows the combined time series of the iGrav residuals and differences between iGrav006 and iGrav032 (black line), GWL, geothermal extraction (red dots) and reinjection (grey dots). In the initial period, when injection flow is increased to just above 120 kg s^− 1 (May 2017 until January 2018), we observe only slight increase of GWL of about 2 m. After January 2018, when injection flow is increased above 160 kg s^− 1, we can see much stronger increasing GWL by about 6 m until July 2018. For the subsequent periods of reduced geothermal flow rates (extraction and injection drops marked by dashed arrows in Fig. 7) we notice an immediate descending response in the GWLs followed by a rebound when the flow rates are increased again. The direct response in the GWLs during periods of high extraction and injection flow rates suggests that the subsurface hydrology consists of pressure controlled systems. When the injection flow rates are reduced below 160 kg s^− 1, pressure reduces instantaneously indicating that the natural outflow of the aquifer in the west is in the order of the same amount. We observe similar effects for the gravity time series, like the long-term increasing trend in the differences (iGrav006 - iGrav032) and short-term gravity responses (decrease in the differences and increase in the iGrav015 and iGrav032 residuals; marked by A, B, C and D in Fig. 7) after reduced geothermal flow rates in March to April 2019 and July to September 2019. However, the gravity responses show time delays of several days up to weeks in comparison to the GWL changes and injection rates. For iGrav006 these responses are missing (or barely visible) in the gravity residuals.

Following temporary reduction in extraction and injection rates, the GWL data show a direct response, indicating that the reinjection aquifer is confined and that the observed response is relative to pressure change in the aquifer. The delayed responses in the gravity data indicate that subsequent to the pressure reduction, groundwater is moving. This time delays may also be attributed to a season depended recharge in the area of the extraction wells. In May 2020, injection and extraction were more permanently reduced, and the response of the gravity meters is also delayed. In addition, the gravity difference (iGrav006 - iGrav032) is not raising any more, suggesting that the mass difference may have become stationary between extraction and recharge, and injection and outflow. It would mean that the system may have attained “equilibrium / steady state” conditions. The missing (or barely visible) responses in the iGrav006 residuals suggest that mass changes of reinjection cause a much lower gravity contribution, compared to production, which influences the gravity signals of iGrav015 and iGrav032. In order to understand and quantify those data, one would require to model the hydraulic features and adjust hydraulic parameters such as hydraulic diffusivity and consider other influences, such as the temperature and density of the fluid injected, and the detailed geology of the reservoir, as known from exploration drilling at Þeistareykir (Kewiy 2013; Óskarsson 2015; Gudjónsdóttir 2018).

In a first order approximation, from the GWL data observation, it seems that the “equilibrium / steady state” conditions for the injection aquifer may be reached when injection rate is around 140 kg s^− 1 after about three years of power plant operation. Above that injection value, the pressure in the aquifer increases. This indicates that the aquifer is a partly open system with additional storage capacities for the reinjected fluids. Correlation between changes in reinjection and GWL changes west of the Tjarnarás fault system (yellow dashed line in Fig. 1) suggest that there is good pressure connection in the western compartment of the Þeistareykir geothermal reservoir. When injection rates are very low, the aquifer system tends to return to the initial state, prior the start of reinjection. Extraction rate changes (even below 180 kg s^− 1) are not clearly visible in the GWL data but have a “signature” in the gravity residuals at the stations east of the reinjection area with time delays of some days to some weeks. This is an indication, that there is a reservoir boundary between the reinjection area in the west and the extraction area in the east. The location of this boundary may be indicated by the surface appearance of the Tjarnarás fault (see Fig. 1). The data after June 2020, when power production was reduced by approximately 30% and extraction rates are in the order of 220 kg s^− 1 (and ~ 140 kg s¹ injection) suggests “steady state conditions” for such operation modus after about 2.5 years of operation. This hypothesis is also supported by the gravity observations in Fig. 7: The difference between the gravity variations at the west and the east stations that was continuously increasing changes slope and remains more or less constant for the time interval observed after extraction decrease (June 2020).

Our continuous gravity monitoring results were also used in a hybrid gravity study including measurements of the three iGravs and a Scintrex CG5 gravity meter at Þeistareykir (Portier et al. 2021). Both sensor types coincide in the observation of long-term gravity decrease at the extraction zone and minor gravity changes at the injection zone. As benefit from the combined hybrid gravity method we obtain gravity observations in high temporal resolution from the iGravs and high spatial resolution from the CG5.

We deployed and maintained a network of multiparameter stations including three iGravs superconducting gravity meters for more than three years at the Þeistareykir geothermal field in North Iceland. This allows us to monitor and detect mass changes associated with environmental and geothermal phenomena. The continuous gravity records were reduced from environmental contributions by means of accurate measurements and modelling of global effects, hydrology, snow and vertical surface displacement. In the case of Theisrareykir, the amplitudes never exceeded 10 µGal (except tides). Residuals were interpreted in terms of groundwater level and geothermal flow rates. The analysis of the gravity residuals with time revealed:

Gravity increase of 3.2 µGal yr^-1 at the site above injection (as expected due to the injected water mass).
Gravity decrease of 8.4 µGal yr¹ at the site above the extraction of mass.
Direct correlation between groundwater levels and injection flow rates, which indicates that the reinjection system is subject to a confined aquifer.
Delayed responses in the gravity signals only above the extraction zone, indicating that there is a reservoir boundary between the reinjection and the extraction areas.
The observed responses in groundwater levels and gravity indicate that the hydraulic response of the hydrothermal system is mainly governed by a partially filled aquifer up to injection flow rates of 140 kg s^-1 for which the system may have reached steady state conditions between production, natural recharge, injection and outflow. Above this value the system is over-pressurised.

This study demonstrates that a network of multiple gravity stations is useful for monitoring of a geothermal reservoir and for finding steady state values. In addition, it shows that observations of environmental parameters are fundamental to obtain accurate estimations of geothermal signals in volcanic environments. This study may also be useful for understanding gravity changes associated to active volcanoes, where environmental parameters and mass inflow are not known.

DEM: digital elevation model

EOP: Earth orientation parameter file

FAG: free-air vertical gravity gradient

G: Universal gravitational constant

Δg_atm: gravity contributions from local atmospheric mass changes

Δg_drift: gravity contributions from instrumental drift

Δg_geoth: gravity contributions from geothermal mass changes

Δg_glob: gravity contributions from global effects

Δg_hydr: gravity contributions from local hydrology

Δg_pm: gravity contributions from polar motion

Δg_snow: gravity contributions from snow

Δg_tide: gravity contributions from solid Earth and ocean tides

δg/δh: vertical gravity gradient

GHE: global hydrology effect

GNSS: global navigation satellite system

GWL: groundwater level

Δh: vertical surface displacement

HyGra: hydro-gravitational modelling

IERS: International Earth Rotation and References Systems Service

M: Earth mass

Max: maximum

Min: minimum

NTOL: nontidal ocean loading

R: Earth radius

ROMPS: remotely operated multiparameter stations

SD: standard deviation

SPA: snow pack analyser

SSG: snow scale

SWE: snow water equivalent

Availability of data and materials

The datasets analysed during this study are available from the corresponding author on reasonable request. A persistent web link in the GFZ repository will be implemented after embargo dates of the follow-up project in 2024.

Competing interests

The authors declare that they have no competing interests.

Funding

Funding for this project was provided by the German Federal Ministry for Education and Research (BMBF, grant: 03G0858A), the Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences and Landsvirkjun.

Authors’ contributions

FF: corresponding author, writing original draft, field work, data acquisition, formal analysis, methodology, visualisation, validation

AG: methodology, conceptualisation, supervision

PJ: field work, conceptualisation, supervision

MR: field work, methodology

BM: methodology, visualisation

JH: field work, methodology

KE: field work, data acquisition, project administration

All authors read and approved the final manuscript.

Acknowledgements

We gratefully acknowledge Anette Kærgaard Mortensen and her colleagues from Landsvirkjun the National Power Company of Iceland for the friendly collaboration, supportive provision of borehole monitoring data and for providing the necessary facilities and infrastructure at the Þeistareykir geothermal field. Many thanks to GWR Instruments, Inc. for their ongoing assistance during operation and maintenance of the iGrav superconducting gravity meters. We also thank Tilo Schöne, Nico Stolarczuk, Julia Illigner and Cornelia Zech for installation of the remotely operated multiparameter stations at Þeistareykir.

Authors’ information

Affiliations

FF, KE: Section 4.8 Geoenergy, GFZ German Research Centre for Geosciences, Telegrafenberg, 14473, Potsdam, Germany

AG, MR: Section 4.4 Hydrology, GFZ German Research Centre for Geosciences, Telegrafenberg, 14473, Potsdam, Germany

PJ: Section 2.2 Geophysical Imaging, GFZ German Research Centre for Geosciences, Telegrafenberg, 14473, Potsdam, Germany

BM: Section 1.1 Space Geodetic Techniques, GFZ German Research Centre for Geosciences, Telegrafenberg, 14473, Potsdam, Germany

JH: ITES (Institut Terre et Environnement de Strasbourg), Université de Strasbourg/CNRS, 5 rue René Descartes, 67084 Strasbourg Cedex, France

Corresponding author

Correspondence to Florian Forster: [email protected]

Agnew DC. Earth tides, in Treatise on Geophysics. 2nd ed. Vol. 3: Geodesy, pp. 151–178, ed. Herring T. Schubert G, editor-in-chief, Elsevier. 2015. https://doi:10.1016/b978-0-444-53802-4.00058-0.
Ármannsson H, Gudmundsson Á, Steingrímsson BS. Exploration and development of the Krafla geothermal area. Jökull. 1987;37:13–30. https://www.researchgate.net/profile/Asgrimur-Gudmundsson/publication/284698694_Exploration_and_development_of_the_Krafla_geothermal_area/links/56e9697208ae47bc651c716e/Exploration-and-development-of-the-Krafla-geothermal-area.pdf. Accessed 02 July 2021.
Boy JP, Hinderer J. Study of the seasonal gravity signal in superconducting gravimeter data. J Geodyn. 2006;41(1):227–33. https://doi.org/10.1016/j.jog.2005.08.035.
Brown GC, Rymer H. Microgravity monitoring at active volcanoes: A review of theory and practice. Cahier du Centre Européen de Géodynamique et de Séismologie. 1991;4:279–304.
Carbone D, Cannavò F, Greco F, Reineman R, Warburton RJ. The benefits of using a network of superconducting gravimeters to monitor and study active volcanoes. Journal of Geophysical Research: Solid Earth. 2019;124(4):4035–50. https://doi.org/10.1029/2018JB017204.
Crossley D, Hinderer J, Boy JP. Regional gravity variations in Europe from superconducting gravimeters. J Geodyn. 2004;38(3–5):325–42. https://doi.org/10.1016/j.jog.2004.07.014.
Damiani TM. Dynamic Effects in Gravimetry: An Assessment of the Current State of Knowledge. NOAA-National Geodetic Survey, 1315 East-West Hwy, SSMC3, Silver Spring, MD 20910, G51B-0362. 2014. https://ftp.ngs.noaa.gov/web/science_edu/presentations_archive/files/damiani_agu_2014_36x65.pdf. Accessed 02 July 2021.
Dobslaw H, Bergmann-Wolf I, Dill R, Poropat L, Thomas M, Dahle C, Flechtner F. A new high-resolution model of non-tidal atmosphere and ocean mass variability for de-aliasing of satellite gravity observations: AOD1B RL06. Geophys J Int. 2017;211(1):263–9. https://doi.org/10.1093/gji/ggx302.
Drouin V, Sigmundsson F, Li S. Ground deformation at the Theistareykir volcanic system, Iceland 535 following onset of geothermal utilization. Proceedings World Geothermal Congress, Reykjavik, Iceland. 2020.
Gendt G, Dick G, Söhne W. GFZ Analysis Center of IGS – annual report 1998, in: Gowey, K., Neilan, R., Moore, A, editors, IGS 1998 Technical Reports. Jet Propulsion Laboratory, Pasadena. 1999;79–87. https://boris.unibe.ch/80305/8/2013_techreport.pdf. Accessed 02 July 2021.
Gitlein O, Timmen L, Müller J. Modeling of atmospheric gravity effects for high-precision observations. International Journal of Geosciences. 2013;4(4):663–71. https://doi.org/10.15488/1531.
Grossi G, Lendvai A, Peretti G, Ranzi R. Snow Precipitation Measured by Gauges: Systematic Error Estimation and Data Series Correction in the Central Italian Alps. Water. 2017;9(7):461. https://doi.org/10.3390/w9070461.
Gudjónsdóttir SR, St. Ásgeirsdóttir R, Sigurgeirsson M, Gudmundsson Á. Drilling of nine high temperature wells in the Þeistareykir geothermal field, NE-Iceland, 2016–2017. Overview of the project, geology of the area and interpretation of the lithology and alteration of the sub-surface strata. In EGU General Assembly Conference Abstracts. 2018;1369.
Gudmundsson A, Brynjolfsson S, Jonsson MT. Structural analysis of a transform fault-rift zone junction in North Iceland. Tectonophysics. 1993;220(1–4):205–21. https://doi.org/10.1016/0040-1951(93)90232-9.
Hinderer J, Crossley D, Warburton RJ. Superconducting gravimetry, in Treatise on Geophysics, 2nd ed. Vol. 3: Geodesy, pp. 66–122, ed. Herring T, Schubert G, editor-in-chief, Elsevier. 2015. https://doi.org/10.1016/B978-0-444-53802-4.00062-2.
Johnston G, Riddell A, Hausler G. The International GNSS Service. In: Teunissen PJ, Montenbruck O, editors. Springer Handbook of Global Navigation Satellite Systems. Springer Handbooks. Cham: Springer; 2017. https://doi.org/10.1007/978-3-319-42928-1_33.
Jousset P, Dwipa S, Beauducel F, Duquesnoy T, Diament M. Temporal gravity at Merapi (Java, Indonesia) during the 1993–1995 crisis: an insight into the dynamical behaviour of volcanoes. J Volcanol Geotherm Res. 2000;100:289–320. https://doi.org/10.1016/S0377-0273(00)00141-4.
Kewiy WR. Injection and production well testing in the geothermal fields of southern Hengill and Reykjanes, SW-Iceland and Theistareykir, N-Iceland. Geothermal Training in Iceland. 2013;747–768. http://os.is/gogn/unu-gtp-report/UNU-GTP-2013-31.pdf. Accessed 02 July 2021.
Klügel T, Wziontek H. Correcting gravimeters and tiltmeters for atmosphericmass attraction using operationalweather models. J Geodyn. 2009;48(3–5):204–10. https://doi.org/10.1016/j.jog.2009.09.010.
Leirião S, He X, Christiansen L, Andersen OB, Bauer-Gottwein P. Calculation of the temporal gravity variation from spatially variable water storage change in soils and aquifers. J Hydrol. 2009;365(3–4):302–9. https://doi.org/10.1016/j.jhydrol.2008.11.040.
Männel B, Brandt A, Bradke M, Sakic P, Brack A, Nischan T. Status of IGS Reprocessing Activities at GFZ. International Association of Geodesy Symposia. Springer. 2020. https://doi.org/10.1007/1345_2020_98.
Männel B, Brandt A, Bradke M, Sakic P, Brack A, Nischan T. GFZ repro3 product series for the International GNSS Service (IGS). GFZ Data Services. 2021. https://doi.org/10.5880/GFZ.1.1.2021.001.
Merriam JB. Atmospheric pressure and gravity. Geophys J Int. 1992;109(3):488–500. https://doi.org/10.1111/j.1365-246X.1992.tb00112.x.
Mikolaj M, Meurers B, Güntner A. Modelling of global mass effects in hydrology, atmosphere and oceans on surface gravity. Comput Geosci. 2016;93:12–20. https://doi.org/10.1016/j.cageo.2016.04.014.
Mikolaj M, Reich M, Güntner A. Resolving Geophysical Signals by Terrestrial Gravimetry: A Time Domain Assessment of the Correction-Induced Uncertainty. Journal of Geophysical Research: Solid Earth. 2019;124(2):2153–65. https://doi.org/10.1029/2018JB016682.
Óskarsson F. Exploration and development of a conceptual model for the Theistareykir geothermal field, NE-Iceland. Short Course VII on Surface Exploration for Geothermal Resources”, organized by UNU-GTP and LaGeo, in Santa Tecla and Ahuachapán, El Salvador. 2015. . Accessed 02 July 2021.
Petit G, Luzum B. IERS conventions (2010). Bureau International des Poids et mesures sevres (france). 2010. https://apps.dtic.mil/sti/citations/ADA535671. Accessed 02 July 2021.
Portier N, Hinderer J, Riccardi U, Ferhat G, Calvo M, Abdelfettah Y, Heimlich C, Bernard J-D. Hybrid gravimetry monitoring of Soultz-sous-Forêts and Rittershoffen geothermal sites (Alsace, France). Geothermics. 2018;76:201–19. https://doi.org/10.1016/j.geothermics.2018.07.008.
Portier N, Hinderer J, Drouin V, Sigmundsson F, Schäfer F, Jousset P, Erbas K, Magnusson I, Hersir GP, Águstsson K, De Z Van Dalfsen, Guðmundsson E, Bernard Á J-D. Time-lapse Micro-gravity Monitoring of the Theistareykir and Krafla Geothermal Reservoirs (Iceland). Proceedings World Geothermal Congress, Reykjavik, Iceland. 2020.
Portier N, Forster F, Hinderer J, Erbas K, Jousset P, Drouin V, Li S, Sigmundsson F, Magnússon I, Hersir GP, Ágústsson K, Guðmundsson Á, Júlíusson E, Hjartasson H, Bernard J-D. Hybrid microgravity monitoring of the Theistareykir geothermal reservoir (North Iceland). Pure and Applied Geophysics. 2021. submitted.
Reich M, Mikolaj M, Blume T, Güntner A. Reducing gravity data for the influence of water storage variations beneath observatory buildings. Geophysics. 2019;84(1). https://doi.org/10.1190/geo2018-0301.1.
Rodell M, Houser PR, Jambor U, Gottschalck J, Mitchell K, Meng CJ, Toll D. The Global Land Data Assimilation System. Bull Am Meteor Soc. 2004;85(3):381–94. https://doi.org/10.1175/BAMS-85-3-381.
Rosat S, Hinderer J. Limits of detection of gravimetric signals on earth. Sci Rep. 2018;8(1):15324. https://doi.org/10.1038/s41598-018-33717-z.
Schäfer F, Jousset P, Güntner A, Erbas K, Hinderer J, Rosat S, Voigt C, Schöne T, Warburton RJ. Performance of three iGrav superconducting gravity meters before and after transport to remote monitoring sites. Geophys J Int. 2020;223(2):959–72. https://doi.org/10.1093/gji/ggaa359.
Schöne T, Zech C, Unger-Shayesteh K, Rudenko V, Thoss H, Wetzel H-U, Gafurov A, Illigner J, Zubovich A. A new permanent multi-parameter monitoring network in Central Asian high mountains – from measurements to data bases. Geosci Instrum Method Data Syst. 2013;2:97–111. https://doi.org/10.5194/gi-2-97-2013.
Sturkell E, Einarsson P, Sigmundsson F, Geirsson H, Olafsson H, Pedersen R, De Zeeuw Van Dalfsen E, Linde AT, Sacks SI, Stefánsson R. Volcano geodesy and magma dynamics in Iceland. J Volcanol Geoth Res. 2006;150(1–3):14–34. https://doi.org/10.1016/j.jvolgeores.2005.07.010.
Van Camp M, Vauterin P. Tsoft: graphical and interactive software for the analysis of time series and Earth tides. Comput Geosci. 2005;31(5):631–40. https://doi.org/10.1016/j.cageo.2004.11.015.
Voigt C, Schulz K, Koch F, Wetzel K-F, Timmen L, Rehm T, Pflug H, Stolarczuk N, Förste C, Flechtner F. Introduction of a Superconducting Gravimeter as Novel Hydrological Sensor for the Alpine Research Catchment Zugspitze. Hydrol Earth Syst Sci Discuss. 2021. https://doi.org/10.5194/hess-2021-78. in review.
Watlet A, Van Camp M, Francis O, Poulain A, Rochez G, Hallet V, Quinif Y, Kaufmann O. Gravity monitoring of underground flash flood events to study their impact on groundwater recharge and the distribution of karst voids. Water Resour Res. 2020;56(4). https://doi.org/10.1029/2019WR026673.
Wenzel HG. The Nanogal software: Earth tide data processing package ETERNA3.30, Bull. Inform.Marees Terrestres. 1996;124:9425–9439. http://www.eas.slu.edu/GGP/ETERNA34/MANUAL/ETERNA33.HTM. Accessed 02 July 2021.

Appendix.docx

Environmental And Anthropogenic Gravity Contributions At The Þeistareykir Geothermal Field, North Iceland

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Methods For Quantifying The Gravity Contributions

Gravity Contributions And Geothermal Signals

Geothermal Interpretation

Conclusions

Abbreviations

Declarations

References

Supplementary Files

Status:

Journal Publication

Version 1