1 Introduction

Multi-annual gravity time series offer a unique, noninvasive, way to monitor mass redistributions in the Earth. However, the gravity value at any time is always the sum of several sources such as solid Earth tides and other mass redistributions both at the surface and inside the Earth, e.g., oceanic and continental water, atmosphere or magma, to name a few. Therefore, focusing on one source require to remove all other potential sources of gravity variations. Despite remaining an active field of research (e.g., Sulzbach et al. 2022; Sun et al. 2023), the present methods to remove the effects of tides and atmospheric loading in gravity time series are efficient (Boy and Chao 2005; Schüller 2015). Thus, the residual gravity variations (i.e., without tidal and atmospheric loading) can be used to monitor other geophysical processes that may influence gravity. Difficulties arise when one want to decipher such processes and assess their respective amplitude in gravity residuals. For instance, gravimeters set in active volcanic regions are influenced at the same time by water and magma redistribution in the ground, but their respective contribution is indistinguishable in the gravity data. In such a case, a proper quantification of magma redistribution is necessarily linked to the assessment of the groundwater effect on the gravity measures (Kazama and Okubo 2009; Mouyen et al. 2016; Carbone et al. 2019).

Hydrological effects on gravity are notoriously challenging to assess (Van Camp et al. 2017), yet they contribute significantly to temporal gravity variations, whether they are measured by satellite- (Ramillien et al. 2008; Feng et al. 2018), or by ground-based gravimeters (Champollion et al. 2018; Hinderer et al. 2020). On the other hand, areas where groundwater redistribution is expected to be the unique source of gravity variations effects are useful to evaluate methods that aim to remove hydrological contribution from gravity time series, as we do in this study. We use a 11-years long gravity time series acquired by a superconducting gravimeter (SG) located at the Onsala Space Observatory, in Sweden (serial number SG054, Fig. 1a, Scherneck et al. (2022)).

Fig. 1
figure 1

a Location of the superconducting gravimeter SG054 at the Onsala Space Observatory. The sub-rectangular area around the instrument outlines the limit of the umbrella effect (see Sect. 2). One graduation on the topographic map is 500 m. b Aerial photography over the umbrella area outlined in red. The dashed red circle shows the edifice where the SG054 is located

Full size image

The objective of this study is to evaluate the relative performance of gravity variations computed from various global hydrological models, specifically GLDAS2, MERRA2 and ERA5 and satellite gravimetry data (GRACE), provided by the EOST Loading Service (Boy 2015), also combined with local, in situ, hydrological measurements. Hereafter, we use the term “hydrological loadings” to denote the effect of continental water on the gravity; it combines both the effects of the Newtonian attraction and of the Earth deformation induced by such mass redistributions.

First, we describe how we obtain the gravity residuals from the original gravity time series acquired by the gravimeter. We also introduce the hydrological loadings models and the groundwater time series that will be used in the analysis of the gravity residuals. We then describe and compare different methods to reduce the effect of groundwater redistribution in the gravity residuals. Finally, we discuss our results and the perspectives they open.

2 Site description

The Onsala Space Observatory is located on the south-west coast of Sweden. The gravimeter is installed in a rectangular edifice with sides measuring 8 by 6 m, made on purpose to host the instrument. During the construction, the soil layer was removed and the pillar on which the gravimeter stands was built onto the basement rocks. The basement is essentially composed of gabbro and diorite (SGU 2023a) that can withstand fracturing and have low value of hydraulic conductivity, ranging between \(10^{-9}\) and \(10^{-7}\) m \(\hbox {s}^{-1}\) (Knutsson 2008; Olofsson et al. 2001; SGU 2023b). These are considered as “non-aquiferous rocks or with local aquifers of minor spatial extension” (Kitterød et al. 2022). Such aquifers are referred to as fractured aquifers, where water circulates and is stored in fracture systems that are not well developed, with fracture porosity less than 1% (Knutsson 2008). In addition, the gravimeter house and the fact that about 5500 \(\hbox {m}^2\) of ground surface around the gravimeter is covered with asphalt are creating an umbrella effect (Creutzfeldt et al. 2008), preventing groundwater recharge and soil water storage directly below and in the vicinity of the gravimeter. This area is outlined and detailed in Fig. 1a, b, respectively.

Below we do a simple gravity modeling that will provide useful insight for a later discussion on the local hydrological effects on gravity at Onsala. Using a digital elevation model (DEM) with a horizontal spatial resolution of 2 m (Lantmäteriet 2023) and masking out the areas that constitute the umbrella, we compute water-thickness to gravity admittance factors when gravity variations are measured at the SG. The computation is done for a thin water layer mapped onto the topography (surface layer) and about 3.4 m below the SG, where stands the average head of the water layer (groundwater layer), according to piezometer measures done in a well that was drilled below the gravimeter house (more details about the piezometer data will be given in Sect. 3.3). The layer is discretized into prisms of 2-m spatial resolution (same as the DEM), and the gravity effect of these prisms is computed. For the surface layer, the computation is done with and without the umbrella effect. There is no umbrella case for the layer at 3.4 m depth since the piezometer shows water table variations, proving that underground lateral water fluxes are possible under the gravimeter houses. In practice, the computation is done for increasing radii of integration around the superconducting gravimeter, to have an upper bound of the radius of influence of the local hydrological signal. The results are shown in Table 1. They give us more accurate quantification than the Bouguer plate approximation (4.19 nm \(\hbox {s}^{-2}\) \(\hbox {cm}^{-1}\)), by accounting for topography and umbrella effects.

Table 1 Water layer to gravity admittance factors summary
Full size table

The admittance factor for a surface layer of water, taking into account the umbrella effect, is 0.10 nm \(\hbox {s}^{-2}\) per cm of water. Because of this very low value, the effect of water stored in the soil (as soil moisture) to the gravity variations will be significantly diminished.

Fig. 2
figure 2

a Gravity residuals using the drift computed within the IGETS L3 workflow. b Instrumental drifts adjusted to the IGETS level-3 (L3) gravity time series (the drift correction computed within the IGETS level-3 workflow is not applied). c Residuals after the refined drift correction

Full size image
Fig. 3
figure 3

Gravity residuals \(g_{\textrm{hydres}}\) with an atmospheric loading correction using ERA5 and assuming an inverted barometer ocean response (blue) or a dynamic ocean response (TUGO-m model) to the air surface pressure (orange), a in the time domain, b in the frequency domain

Full size image

The deeper layer of water, one that would monitored by the piezometer, has a larger admittance factor, slightly larger than the Bouguer plate approximation. But this factor is computed assuming a plain water layer, 1 cm thick, at 3.4 m depth. According to the previous hydrogeological description of the area, such a wide layer is unlikely to exist. Groundwater may rather circulate and be stored only in a few fractures scattered over the area.

3 Materials

3.1 Gravity data preparation

We use the level-3 gravity residual time series provided by the International Geodynamics and Earth Tide Service (IGETS) (Boy et al. 2020), with a few modifications. The level-3 residuals are obtained from the SG054 gravity time series after correction for spikes, solid and oceanic tides, polar motion and length-of-day, atmospheric loading and instrumental drift. The time series originally has a time resolution of 1 min, which we down-sample to 1 h, and runs from August 2009 to February 2021.

The original level-3 drift correction is replaced by a combination of linear and exponential drifts that suits the data better (Fig. 2). In particular, the new drift accounts for a drift break that occurred on 22 February 2011, when the control card of the instrument was replaced (Scherneck and Rajner 2019). We thus split the drift estimation on two time periods before and after that break: August 2009–February 2011 (linear) and February 2011–February 2021 (linear + exponential). Note that this drift estimation also includes a linear effect due to the global isostatic adjustment (GIA), which is estimated to be \(3.0 \pm 0.7\) nm \(\hbox {s}^{-2}\) \(\hbox {year}^{-1}\) (Olsson et al. 2019) at the SG054 location. We do not need to separate the linear instrumental drift from the GIA for our purpose.

The original atmospheric loading correction uses ERA5 hourly surface atmospheric pressure (Hersbach et al. 2020) and assumes an inverted barometer ocean response to that pressure. This correction is replaced with an atmospheric loading correction still based on ERA5 surface pressure but that assumes a dynamic response of the ocean to the air pressure and wind, using the TUGO-m model (Carrère and Lyard 2003). In addition, the local component of this correction, that is within a radius of 0.1\(^{\circ }\) around the gravimeter, is replaced by local surface pressure measurements converted to gravity units using an admittance of \(-2.2105\) nm \(\hbox {s}^{-2}\) \(\hbox {hPa}^{-1}\) (Boy 2015). This approach was shown to be better than the inverted barometer model, allowing a precise estimate of the ocean non-tidal loading contributions of major storm surges (Boy and Lyard 2008). That is also noticeable in our results (Fig. 3), for the frequency between \(10^{-2}\) and 2 cycles per day (cpd), where the noise power of the residual is up to 20 times lower when considering a dynamic response of the ocean to the air pressure instead of the inverted barometer hypothesis. The benefit of the TUGO-m model is exacerbated here by the proximity of the sea, which is only a few hundreds of meters from the superconducting gravimeter (Fig. 1a).

At this stage, the residual gravity time series \(g_{\textrm{hydres}}\) predominantly shows the effect of water storage variations on gravity. In the rest of this study, we will focus on removing such hydrological residuals.

3.2 Hydrological loading models

We use three hydrological loading products computed by the EOST loading service (Boy 2015):

  1. 1.

    The Global Land Data Assimilation System (GLDAS) version 2.1 with the Noah land surface model (Rodell et al. 2004). Space and time resolution: 3 h, 0.25\(^{\circ }\).

  2. 2.

    The Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA2) model (Gelaro et al. 2017; Reichle et al. 2017). Time and time spatial: 1 h, 0.625\(^{\circ }\).

  3. 3.

    The ECMWF Re-Analysis (ERA5) (Balsamo et al. 2009). Time and spatial resolution: 6 h, 0.7\(^{\circ }\).

In addition, we use also a modified version of the GRACE-derived hydrological loading, computed from the global iterated global mascons RL06v1.0 (Luthcke et al. 2013; Loomis et al. 2019), with monthly time resolution and 1\(^{\circ }\) spatial resolution. That modified version uses only the land mascons to compute the hydrological loading (such land regions are readily available from the mascons RL06v1.0 dataset). Among these four hydrological loadings, only the GRACE-derived one account for total water storage variations, including groundwater storage variations (Swenson and Wahr 2006).

Fig. 4
figure 4

Hydrological products of the EOST loading service used in this study, with their a total, b non-local and c local contribution to the gravity measured at the location of the SG054 gravimeter. The scale of the y-axis is kept identical to facilitate the comparison of each signal. In a, we also add the SG054 gravity residuals, which should only contain hydrological effects

Full size image

All models are in gravity units (nm \(\hbox {s}^{-2}\)). Their time resolution is changed to 1-h, expect MERRA2 since it is already sampled at 1 h. Thus, all hydrological models can be directly compared to the gravity residuals (Fig. 4). They account for both the gravitational attraction and the elastic deformation of the Earth due to the redistribution of water masses over the continents. Computational details are described on the EOST loading service web-page and references therein (e.g., Boy et al. 2009). Each hydrological product is divided into three parts: local, non-local and total. The local part is the Newtonian attraction of a layer of water mapped on the surface around the gravimeter, which is computed assuming an effect of \(a_{N} = 4.2677\hbox { nm s}^{-2}\hbox { cm}_{\textrm{water}}^{-1}\) (Boy 2015). Figure 4c shows that the local part maximal amplitude is about \(100\hbox { nm s}^{-2}\), which corresponds to about \(100/4.2677=23\hbox { cm}\) of equivalent water height. Such a layer reaches its maximum effect in a radius of 100 m a gravimeter, further than that, the Newtonian attraction of this water layer is negligible. Ideally, the local part should be replaced by measurements of all relevant storage compartments in the vicinity of the gravimeter, such as soil moisture, water content in the unsaturated zone, snow, and eventually surface water bodies such as lakes. In our case, we will only add the groundwater contribution, which hydrological models do not take into account. The total (global) part combines the elastic deformation of the crust and the Newtonian attraction effect due to the land surface water at the global scale. It is computed by the convolution of gravity Green functions with a hydrological model. The non-local part is obtained by subtracting the local part from the total part.

We see in Fig. 4a that all hydrological products overestimate the gravity residuals measured at Onsala. In fact, the amplitude of the gravity residuals is much closer to that of the non-local part of the hydrological products. This already suggests that local hydrology only has a minor influence on the gravity records.

3.3 In situ groundwater level data

In situ groundwater level data are recorded with a piezometer set at 6 m under the gravimeter house, in a water well drilled at about 2.2 m next to the SG054. Groundwater level data are available since 31 August 2015, but the SG054 records started in June 2009. For completeness, we use an artificial neural network (ANN) algorithm to model the missing groundwater level data, between June 2009 and August 2015. Modeling missing data is called data imputation and is done here using a Long Short-Term Memory (LSTM) recurrent neural network (Hochreiter and Schmidhuber 1997), a type of ANN (Géron 2019) that can be thought as an elaborated regression method making the link between groundwater level variations and at set of other meteorological data. Among the meteorological data available at the Onsala observatory since before the SG was set, the air temperature, pressure, relative humidity and accumulated rainfalls over the last 24 h show the largest coefficient of correlation with the groundwater level time series. The rainfalls are responsible for water inputs in the ground, while temperature, pressure and relative humidity are influencing evapotranspiration processes (Zhang et al. 2016), hence water outputs. The amount of water in the ground as measured by the piezometer is the balance between such water inputs and outputs. We use these parameters (air temperature, pressure, relative humidity and accumulated rainfalls over the last 24 h) to model the missing groundwater table variations in the period 2009–2015 (Fig. 5).

Fig. 5
figure 5

Hydro-meteorological data used in this study. a 24 h-cumulative rainfalls. b Relative humidity of the air. c Air temperature. d Air pressure. e Groundwater level variations (height of the water-table above the piezometer, which is located 6 m below the pillar where stands the gravimeter, in a well that was drilled 2.2 m next to it)

Full size image
Fig. 6
figure 6

a Train and test fits to the groundwater level data obtained with the artificial neural network. b Random close-up on the test fit, where the algorithm predicts groundwater level variations

Full size image

The first step is to train the ANN to return the measured groundwater level using air temperature, pressure, relative humidity and accumulated rainfalls over the last 24 h as inputs. Only the first 70% of the groundwater data are made available to the ANN, which makes its own model with the aim to fit the groundwater data using the meteorological inputs. In a second step, the remaining 30% of the groundwater data are then used as a test set: the ANN applies the model it just created to the inputs and the resulting groundwater level is compared to the observed one. It returns a root mean squared error (RMSE) that quantify how well the model fits the actual data. A few user-defined parameters are necessary to do such an ANN study but, to properly predict the groundwater level at a given time t, we find that the most influential parameters is how long before t must start the time series of meteorological inputs. By running systematic tests on that parameter, we eventually reach the lowest RMSE when using 190 h of meteorological data before the time of groundwater level prediction. The RMSE is 5 cm, which is about 10% of the maximum amplitude of the groundwater table variations. Comparing the model to the test set (Fig. 6), we see that our ANN is able to retrieve the main features of the observed groundwater variations:

  1. 1.

    Immediate increase in the water level during the rain.

  2. 2.

    Rapid/exponential decrease after the rain.

  3. 3.

    Seasonal variations.

4 Removal of the hydrological effects in gravity time series

Given the corrections applied so far, we proceed with the hypothesis that hydrology is the only geophysical signal remaining in the gravity residuals. In this section, we first remove the non-local part of the hydrological signal, for all the models described in Sect. 3.2. We then introduce several ways to further remove the local hydrological effects, using the local part of the models from Sect. 3.2, but also the local measurements of the groundwater level measurements introduced in Sect. 3.3.

4.1 Non-local effect of the hydrology on gravity variations

The non-local part of hydrological loading is subtracted from the gravity residuals. This is done for the four hydrological loadings, namely GLDAS2, ERA5, MERRA2 and GRACE (land contribution). The results are shown in Fig. 7, and the impact on the gravity residuals is summarized in Table 2.

Fig. 7
figure 7

a Time series of the gravity residuals before any hydrological correction (\(g_{\textrm{hydres}}\), top blue line) and of the gravity residuals after removing the non-local part of hydrological loading computed from GLDAS2 (labelled GLD2), ERA5, MERRA2 (labelled MER2) and GRACE. All series are offset for legibility. b Same as (a) but represented as the power spectrum density (PSD) and emphasizing the seasonal band, i.e., the periods between 346 and 415 days, where the gravity residuals without hydrological corrections (blue line) show the largest power

Full size image

Just accounting for the non-local contribution of the hydrological loading significantly reduce the gravity residuals. As one could anticipate from Fig. 4, this reduction is mainly visible on the seasonal component of the gravity residuals. These gravity residuals, corrected from the non-local hydrological loading, are labelled \(g_{\textrm{hydres}\_\textrm{Loc}}\), since they should only contain the local part of the hydrological contribution to the gravity.

4.2 Total hydrological reduction methods

Here, we describe how we evaluate the total (non-local and local) contribution of water redistributions on the gravity data. We propose different methods depending on whether we use (1) GLDAS2, MERRA2 or ERA5, or (2) GRACE. We recall that the local hydrological effect on gravity is the sum of the surface water (mainly the soil moisture) and of the groundwater storage variations. In addition, GLDAS2, MERRA2 or ERA5 do not account for groundwater storage variations while GRACE does. These methods all have the same constrain, that is to minimize the standard deviation of the final gravity residuals after all hydrological corrections, called \(g_{\textrm{res}}\).

Table 2 Performance of the non-local effect of the hydrology on gravity variations computed for different hydrological models at Onsala’s SG
Full size table

4.2.1 Method 1

This method combines the groundwater level data with either GLDAS2, MERRA2 or ERA5. The final gravity residuals \(g_{\textrm{res}}\) corrected from the local and non-local part of the hydrological loading are expressed as

$$\begin{aligned} g_{\textrm{res}} = g_{\textrm{hydres}\_\textrm{Loc}} - k_{\textrm{Loc}} g_{\textrm{Loc}} - k_{\textrm{gw}} a_{\textrm{gw}} \Delta h_{\textrm{gw}} \end{aligned}$$

where \(g_{\textrm{hydres}\_\textrm{Loc}}\) are the gravity residuals corrected only from the non-local contribution of the hydrological loading, as shown in Fig. 7a, \(g_{\textrm{Loc}}\) is the local part of the hydrological loading and \(k_{\textrm{Loc}} = a_{su}/a_N = 0.1/4.2677 = 0.023\) is the ratio between the local hydrological admittance for gravity, computed considering the local topography and the umbrella effect at the SG54 (Table 1), and the local admittance used in the EOST-loading products. It allows to adjust the local hydrological effect computed for each model to the local properties of the site (topography and umbrella effect). \(a_{\textrm{gw}} \Delta h_{\textrm{gw}}\) is the gravitational attraction of groundwater level variation \(\Delta h_{\textrm{gw}}\) in the water well under the gravimeter house, using the admittance factor \(a_{\textrm{gw}} = 4.34\) nm \(\hbox {s}^{-2}\) \(\hbox {cm}_{\textrm{water}}^{-1}\) (Table 1). However, it assumes that the groundwater level change is a plain layer of water, which is impossible since here the groundwater can only fill the fractures of the rock. Hence, the introduction of the factor \(k_{\textrm{gw}}\), which can primarily be interpreted as the aquifer average porosity. We estimate \(k_{\textrm{gw}}\) by minimizing the standard deviation of \(g_{\textrm{res}}\). As we only use GLDAS2, MERRA2 or ERA5 (not GRACE), the local hydrological contribution is computed by adding the surface water contribution (from the models) and the groundwater contributions (from the piezometer measurements). Nevertheless, the surface water contribution being modeled, it may differ from the actual surface water storage changes. Thus, \(k_{\textrm{gw}}\), although interpreted as the porosity, can also tend to adjust to such surface water modeling errors.

4.2.2 Method 2

This method only uses the product of the EOST loading service and thus, does not account for any in situ hydrological observations. We have:

$$\begin{aligned} g_{\textrm{res}} = g_{\textrm{hydres}\_\textrm{Loc}} - k_{\textrm{Loc}} g_{\textrm{Loc}} \end{aligned}$$

4.2.3 Method 3

This third method is almost identical to the first one, except that \(k_{\textrm{Loc}}\) is also estimated, rather than computed from the local topography and the estimated umbrella effects, as in Method 1. Therefore, the minimization of the residuals is constrained by the estimation of the two parameters \(k_{\textrm{Loc}}\) and \(k_{\textrm{gw}}\). Note that as methods 1 and 3 use local observations of the groundwater, they cannot be used with the hydrological products computed from GRACE, which local part already account for groundwater. Thus, only the method 2 uses GRACE products.

Fig. 8
figure 8

Percentage of reduction in the standard deviation (STD) of the gravity residuals achieved after applying the methods described in Sect. 4.2. Each method is tested on all loading product except GRACE, which can only be used in the method 2. The percentage of reduction is computed relative to the standard deviation of the gravity residuals without any hydrological correction is 10.7 nm \(\hbox {s}^{-2}\). Their values are also given in Table 3

Full size image
Table 3 Values of the scaling factor \(k_{\textrm{Loc}}\) and \(k_{\textrm{gw}}\) for the methods listed in Sect. 4.2
Full size table

5 Results

5.1 Performance of each method

The results of the different methods to remove hydrological effects in the SG time series are summarized in Fig. 8 and in Table 3. We report the standard deviation of the gravity residuals but also the PSD reduction in the gravity residuals in the seasonal band (periods between 346 and 415 days). Before analyzing these results in details, it is worth noticing that, comparing the STD reductions reported in Tables 2 and 3, the local hydrological contribution to the gravity measured at Onsala is minor, only improving the STD reduction by 0.1 nm \(\hbox {s}^{-2}\). We elaborate on this in Sect. 5.3.

The best hydrological correction is obtained with any method applied on the ERA5 model, leading to the lowest standard deviation of the residuals (7.8 nm \(\hbox {s}^{-2}\)) and to the highest noise reduction in the seasonal band (89%). This method combines the non-local and local parts derived from ERA5 (surface water) and the groundwater level variations measured below the gravimeter. This is quite representative of our results, since Fig. 8 and Table 3 show that, considering both the standard deviation of the residuals and to the noise reduction in the seasonal band:

  • Overall, hydrological corrections based partly or entirely on ERA5 model perform consistently better than those based on MERRA2, GLDAS2 or fully on GRACE, regardless of the method employed.

  • Focusing on the methods, method 3 always works best. Nevertheless, with ERA5, the improved performance of method 3 over methods 1 and 2 is barely noticeable.

The GRACE hydrological loading products are only used in Method 2 and are also those leading to the poorest reduction in the gravity residuals for that method. A reason is that, although GRACE gives a more complete view of water storage variations than models, because it senses all water sources, its spatial resolution will significantly impede the accuracy of local hydrological effects.

The scale factors on the groundwater effect found with the method 1 are about 0.003 on average, that is interpreted as a 0.3% porosity. Focusing in particular on the results using the ERA5 model, that porosity is 0.2%. This is compatible with the porosity of such fractured bedrock, expected to be less than 1% (Sect. 2 and Knutsson 2008).

But these scaling factors are only the result of numerical optimization and may thus also compensate for mismodelled hydrogeological processes. For instance, in method 3, some of the scaling factors of the surface water \(k_{\textrm{Loc}}\) and groundwater \(k_{\textrm{gw}}\) local are negative. Since both the non-local, local and groundwater gravity effect have strong, in phase, seasonal signatures, such negative factor have the role to dampen each others contributions. With ERA5, \(k_{\textrm{gw}}\) is slightly negative, likely to lower the local contribution and hence converging toward to method 1’s scale factors.

For MERRA2 and GLDAS2, \(k_{\textrm{gw}}\) are larger and positive, leading to porosity values of 2 and 1%, larger than what is expected regarding the local geology. On the other hand, their \(k_{\textrm{Loc}}\) are negative. Interestingly, the \(k_{\textrm{Loc}}\) values found in method 2 for MERRA2 and GLDAS2 are also negative. Yet method 2 only uses the non-local and local contributions of the models, not the groundwater data. Therefore, the non-local contribution computed from MERRA2 and GLDAS2 seems overestimated, and the negative scaling factors tend to reduce them.

Fig. 9
figure 9

a Initial gravity residuals, without any hydrological corrections (\(g_{\textrm{hydres}}\), blue) and the hydrological model computed method 1 and ERA5 (orange). b Gravity residuals after the complete reduction in hydrological effect using method 1 and ERA5 (\(g_{\textrm{res}}\), blue). It is overlayed with a singular spectrum analysis (SSA)-based smoothing done by summing the two first components of the SSA of the residuals, out of 200 (black)

Full size image

GRACE returns the largest \(k_{\textrm{Loc}}\), but still leading to the poorest results in method 2. This suggests that GRACE hydrological loading is not well suited to model the water mass variations at Onsala. Specifically, its local part is underestimated, which may first be explained by the too large spatial resolution of this method compared to the very local sensitivity of the SG (Van Camp et al. 2014).

5.2 Residual gravity using the best hydrological corrections

As seen before, using the ERA5 model, all methods perform equally well. Here, we show the results of method 1 because they were constrained using the local site properties (topography and umbrella effect) and also return a value for the bedrock porosity that is compatible with the local geology. Thus, we consider that method 1 is physically more relevant than the other methods.

The two gravity lows in 2018 and 2019 (Fig. 9a) correlate in time with the respective heatwaves that occurred in northern Europe over the same periods (Boergens et al. 2020; Wilcke et al. 2020) and can logically be attributed to a decrease in the water storage during these periods. The hydrological correction does not entirely match these gravity decrease yet, underestimating it by about 20 nm \(\hbox {s}^{-2}\).

Fig. 10
figure 10

Relative parts of the non-local hydrology (blue), local surface water (orange) and local groundwater (green) to the hydrological models obtained for each method and hydrological loading products

Full size image

5.3 Relative contributions of the non-local and local hydrology to the gravity residuals

Figure 10 shows that for all models, the non-local hydrological contributes most to the gravity residuals at Onsala. For the best models (ERA5; any method), this contribution ranges from about 80 to 90%. The remaining contributions are in majority that of the local surface water and, eventually the local groundwater. This is in agreement with the local settings, where a low-porosity bedrock, soil removal beneath the gravimeter house and umbrella effect all tend to lower local water storage variations, both in the shallow and deeper ground. In method 1, the scaling factor for the local surface water is fixed (\(k_{\textrm{Loc}} = 0.023\)), so the groundwater can only take the remaining contribution. In method 3, however, both \(k_{\textrm{Loc}}\) and \(k_{\textrm{gw}}\) are estimated, giving more freedom to the relative influence of the local surface and groundwater contribution. We have seen in Table 3 that these factors can be negative. Such negative \(k_{\textrm{Loc}}\) or \(k_{\textrm{gw}}\) factors mean that the local surface water or the local groundwater, respectively, do not physically influence the gravity. Instead, they have only an optimization role, mitigating another contribution in order to best reduce the STD of the residuals.

We also observe that, the scaling factor on the groundwater effect \(k_{\textrm{gw}}\) being quite small, the hydrological model can only explains the seasonal component of the gravity signal. Gravity variations at higher frequencies are not lowered by the model. Indeed, Fig. 4 shows that ERA5 is predominantly a seasonal model, both in local and non-local. Only the groundwater table variations have significant signals at periods between 1 h and a few days, which are significantly lowered once \(k_{\textrm{gw}}\) is applied (see the model in Fig. 9a).

6 Discussion

6.1 Perspectives for others sites and other models

Testing these methods at other SG sites of the IGETS network depends on the local hydrological measurements available. If no such data at all, one may just use method 2, which returned good results in our case, when used with ERA5. Yet method 2 does not account for groundwater, because the hydrological models do not contain that information. One could thus take the problem backward, and try to assess groundwater storage variation from the gravity residuals, rather than minimizing the residuals. Considering method 1, it is possible to compute \(k_{\textrm{Loc}}\) for each site, using a local DEM and mapping the umbrella area. In that case, assuming gravity variations (after the correction of tides and atmospheric loading) are only influenced by water redistribution, the residual gravity would represent the effect of groundwater storage variations, plus noise.

For all the methods described in this study, the non-local and local contributions are both taken as pairs from the same sources (ERA5, MERRA2, GLDAS2 or GRACE). However, it would be interesting to systematically correct the non-local contribution using the GRACE mascons (land contribution only). Indeed, GRACE records total water storage variations, including aquifers, which are missing in the hydrological models. Therefore, GRACE is better than the models at sensing the non-local contribution of hydrology to gravity measurements. In addition, the coarser spatial resolution of GRACE compared to the other models is a lesser issue when assessing the non-local hydrological contribution to the gravity. This is because the relative error in the distances of the masses from the gravimeter diminishes as these non-local masses are located further away. Practically, we could use the same methods as described above, but always using the GRACE-derived non-local contribution. The way the EOST loading products are presently computed does not permit such a mixing of models, but that could be a perspective for future studies.

6.2 Inter-annual gravity changes

The black line in Fig. 9b is the smoothed residuals of the gravity time series after removing hydrological loading effects, using method 1 with ERA5. This smoothing is obtained by summing the two first components the singular spectrum analysis (SSA) of the residuals, done to split them into 200 components. The other 198 components represent oscillations and noise. This SSA is only used as a visual smoothing to highlight residual gravity signal that do not have a seasonal component.

We observe a trough of about \(10\hbox { nm s}^{-2}\) amplitude starting in late 2017. The gravity residual time series decreases to a minimum in the middle of 2018 then rises until the end of 2020. This behavior spans the two European heatwaves of 2018 and 2019 but is not bounded by them. Global hydrological models do not properly capture inter-annual (non-seasonal) water storage variations (Scanlon et al. 2018), and our observation might be another example of this issue. Yet, Fig. 4b shows that the non-local part from GRACE does not display such a feature, and GRACE non-local is better than models at capturing inter-annual hydrological variations. Therefore, if of natural origin, this trough would come from local hydrological processes that are not properly described by our model. Another hypothesis could be an instrumental bias, perhaps a succession of offsets, individually unnoticeable but which, put together, reach this 10 nm \(\hbox {s}^{-2}\) amplitude. We dismiss a drift issue because SGs have a low and stable instrumental drift, often modeled by a linear or an exponential function (Van Camp and Francis 2007; Hinderer et al. 2015), and a combination of both in this case, which cannot create such a pattern.

7 Conclusions

Using hydrological products available from the EOST Loading Service (Boy 2015), we proposed and tested various methods to remove hydrological effects in long-term gravity time series measured by Onsala’s superconducting gravimeter. All the methods using ERA5 products, either alone or combined with in situ groundwater level data, returned equally good results that explained 89% of the seasonal hydrological signal and decreased the standard deviation of the gravity residuals by about 3 nm \(\hbox {s}^{-2}\), i.e., 27% of the gravity residuals without hydrological correction. Thus, in this case, even without in situ groundwater level data, it is possible to properly remove hydrological effects for this station. This hydrological correction is yet uniquely working for the seasonal component. Gravity variations of higher frequencies are left unchanged. Hydrological effects are usually very local, and our conclusions for Onsala will surely not be valid for gravity time series located in different areas. However, since the hydrological products on which are based our corrections are available for all SGs that belong to the International Geodynamics and Earth Tide Service (IGETS), testing the approaches introduced in this study to other sites is a reasonable perspective.