1 Introduction

As laid out in Scherneck and Rajner (2019, henceforth SR19) multi-campaign reduction at the drop level using a Superconducting gravimeter (SG) instead of a priori models for the temporal variation of gravity affords us a range of advantages in scrutinizing the performance of FG5 AG’s. Trends that appear unrelated to physical variations of g during the AG setups can be discerned and the SG scale factor obtained more robustly, as Fig. 1 shows. For further reading into FG5 instruments and applications, cf. Niebauer et al. (1995).

For the main points of this paper, we shortly summarize the data and the methods employed to expose the major obstacles on the way towards an early, reasonably uncertain inference of the rate of change of gravity as an enterprise in a postglacial setting at large, and at Onsala Space Observatory (OSO), Sweden, in particular. According to SR19, the hydrological impact at Onsala seems benign, in particular those related to groundwater variations. At other stations, we would expect that an SG would catch such signatures and reduce the AG measurements of them in the multi-campaign approach.

Fig. 1
figure 1

Parameter correlation in the multi-campaign adjustment (bottom) compared with a single campaign (top; the one with greatest duration, 201304a, was chosen). In both cases, setup slopes were estimated. Correlation between monument ties turns out large regardless which, however, calibration factor and setup slope parameters attain high correlation in the latter case, while the multi-campaign case gets much closer to orthogonality in general and on slope parameters in particular. Notably in single-campaign reduction, slopes and scale factor cook a dangerous brew together

Concerning the uncertainties of the result we present writing \(v\;\pm \sigma _v\), the value of propagated class-A errors (due to randomly distributed measurement errors or perturbations) is given by \(\sigma _v\) meaning \(1-\sigma \) deviations, also known as STD (short for standard deviation). Where we diverge from this use of \(\pm \,\sigma _v\), we indicate it in due context. Concerning the term offset, we use it in the sense of relative offset (between instruments and campaigns) throughout this article, with respect to a unique reference value as the observation equation will show.

1.1 Added value to IGRF realization and maintenance

For the maintenance of the International Gravity Reference Frame (IGRF), the International Association of Geodesy (IAG) in the framework of its flagship program GGOS (Global Geodetic Observing System) resolved in 2019 to urge for national reference stations to be established.Footnote 1 The Onsala site’s dedicated gravimetric laboratory co-located with reference stations for IVS and IGS in VLBI and GNSS (International VLBI, Very Long Baseline Interferometry, Service; International GNSS, Global Navigation Satellite Systems, Service, respectively, see Haas et al. 2019) appears suitable, not least with regard to the continuously operating SG and an array of recorded environmental sensors (meteorological station, borehole water level gauge, broadband seismometer and tide gauges, see Elgered et al. 2019 and SR19). More frequent visits with more AG instruments would certainly strengthen its robustness, as would a stationary, continuously recording next-generation quantum gravimeter. The latter could be achieved as a joint effort between Chalmers University of Technology, Lantmäteriet (the Swedish mapping, cadastral and land registration authority), and Research Institutes of Sweden (RISE, including the Swedish body for accreditation in metrology). What we can do in this paper is to buttress this prospect with a careful and, in parts, innovative analysis of the 16 AG campaigns that took place in the Onsala gravimetric laboratory from 2009 to 2018, including one paralleled with the GAIN quantum gravimeter operated by Institute of Physics at Humboldt University, Berlin, for 18 days in February 2015. We also dedicate our efforts to GGOS (Haas et al. 2015).

1.2 Estimating the secular rate of change of surface gravity

The primary contribution of this paper is to determine the secular trend of gravity in the context of maintenance of a gravity reference system (Olsson et al. 2019) and for geophysical research applications (Olsson et al. 2015). The main reason for a secular change, a decrease of little-g, in the region of the Baltic Shield is attributed to glacial isostatic adjustment (GIA). Models predict a value between \(-\,4.7\) and \(-\,3.5\) nm/s\(^2\)/year at this site (Olsson et al. 2012, 2015, 2019). Being located near the coast of Kattegat, sea-level change may play a though secondary role. Outrightly founded on crystalline bedrock, the influence of local water storage is much less there than at many other sites. We find a range of the admitted effect based on ECMWF’s ERAin and ERA5 (Dee et al. 2011; Hersbach et al. 2019) into the record of the SG (GWR OSG-054) at 4.0 nm/s\(^2\) RMS and 22 nm/s\(^2\) peak-to-peak range, which is on the order of the gravity residual in our most comprehensive variant of adjustment of tides and environmental parameters.

As a short introduction to the terminology used for AG measurements, campaigns are arranged as a visit of an FG5 at a station on, for example, a yearly basis lasting a couple of days. Each campaign is divided into so-called projects also known as setups, where the instrument is, for example, moved to different platforms or at least oriented in two opposite directions due north and due south for mitigating Coriolis and Eötvös effects. In each setup, the FG5 drops a prism in free fall hundreds to thousands of times. Usually, a number of drops, say 50, are combined into a set, and the drop interval within a set may be a value between five and 30 seconds, usually the same during the whole campaign. The repetition rate of sets may be chosen as a compromise such that the length of a setup is covered with an amount of drops that keeps wear in the instrument low. The drop interval is crucial in the presence of microseismic noise; if the level is high, the noise peaking in power between three and eight seconds period can be detrimental to the repeatability of sets due to sampling aliasing. The formation of sets is rather a means of convenience, reducing the amount of data and its noise as the set averages are computed. An overreaching criterion for the number of drops per setup is the attained repeatability of the setup’s collected measurements. A typical drop determines g with an uncertainty of at least 200 nm/s\(^2\). In order to arrive at 20 nm/s\(^2\) or less, which is the order of accuracy of an FG5 according to both producer and intercomparison of instruments, a setup should collect at least 1000 drops, a margin that anticipates unfavourable conditions like passing seismic waves. For a review, see Van Camp et al. (2017, chap. 2.1.1).

Table 1 Absolute gravity campaigns at Onsala Space Observatory from 2009 on, after the superconducting gravimeter had been installed

1.3 Data analysis

In this study, we process only the original drop measurements, almost 200,000 in number, in one stroke; we call this strategy henceforth multi-campaign analysis. However, not all drops are admitted. The drop measurements are accompanied by an uncertainty with which we weight design matrix and right-hand side of the observation equation. Drop values that exceed their uncertainty by a factor of three are removed at single-campaign level and those few permil that make their way into the multi-campaign inversion are further down-weighted in order to effectively bar them from influence.

Concerning the SG, the samples we have obtained up to this date (2019-10-27) evaluate to 327,014,220 with 182,580 lost (a 0.56 permil leakage). Forming hourly ordinates and analysing for signals, the series is 74,204 h long with 1423 samples missing due breaks in SG operation, due to outlier rejection (criterion \(4-\sigma \) and reiteration), and as an inevitable consequence of the whitening filter as it widens data gaps.

The very fact that we arrive at a residual RMS of 5–7 nm/s\(^2\) (SR19) is taken as vindication of a well-achieved breakdown of perturbing effects in terms of ancillary observations and models. Key enabling features are that

  • the instrumental drift can be parametrized by a few jumps accompanied by changes in linear rate (four events), and two events followed by an exponential decay signature. The rate uncertainties depend on the lengths of the segments and range between 0.15 and 1.1 nm/s\(^2\)/year;

  • the stochastic noise can be approximated with a low-order prediction error filter (PEF, with less than ten coefficients).

Thus, a strategy to reduce the measurements of the 16 AG campaigns at OSO could be laid out which might appear as unorthodox as bold: all drop observations, 199,542 in total, are reduced using the SG data and its inferred drift function; orientation biases individually for each meter (180\(^\circ \) reorientation from north to south is to eliminate the Eötvös effect, Křen et al. (2018)); a meter bias, i.e. a relative offset between meters (FG5X-233, FG5-220, FG5X-220) and the declared master meter (FG5-233) after applying the so-called degrees of equivalence (DoE) resulting from International Comparison Campaigns (henceforth ICC); monument ties (OSO has four observation points); and setup trend rate for each setup of the AG’s. We refer to the following ICC’s: Jiang et al. (2012), Francis et al. (2013), Francis et al. (2015), Pálinkáš et al. (2017) and Falk et al. (2020). For a review of the concept of DoE and their uncertainties, refer to Burns (2003) and any of the ICC publications. (In short, DoE represents the offset of a specific AG with respect to the ensemble of meters participating in an ICC.) The least-squares fit delivers also the scale factor of the SG and the secular rate of gravity, the latter being afflicted by shortcomings in the cross-hair of this investigation. Table 1 shows campaigns, instruments and drop yield.

Fig. 2
figure 2

AG setup slopes and SG scale factor. The multi-campaign’s solution for the scale factor is shown as the yellow horizontal line outlined by a red box to show its uncertainty. The value determined from the GAIN campaign in parallel with campaign no. 12 (201502b) is shown as a short horizontal grey bar (zero line), its thickness representing the uncertainty; all other scale factor ordinates represent differences to GAIN. Residual, campaign-related signatures are determined from covariation with the multi-campaign residual evaluated within the campaign subsets (red crosses). Also shown are scale factors determined from campaigns one at a time (blue crosses), then not reducing for setup slopes. Setup slopes are shown as small black dots and their error bars (standard deviation) in grey. Campaigns 10 and 11 in 2014 were partially conducted in parallel with the two FG5’s (220 and 233, marked as p.c.). There appears to be a relation between scale factor deviations and the presence of at least one significant setup slope in a campaign (most clearly at no. 16, 201806a)

Fig. 3
figure 3

An example for trends in FG5(X) setups, here in the parallel campaigns 10 and 11 (cf. Fig. 2) labelled 2014mmnn(a,b) (cf. Table 3). The measurements are shown in the two diagrams to the left, AG as black dots with grey error bars, SG as blue dots. In the column to the right, the best-fit linear slopes determined in the multi-campaign analysis are shown in red. The drop measurements reduced by the simultaneous SG data using \(-\,773.18\) nm/s\(^2\)/V (black dots) are shown with their uncertainties (grey bars). In the case of FG5-233, the first and third slope are highly significant (5 resp. 7 standard deviations), whereas none of the trends in FG5X-220 revealed statistically significant slopes. Offsets due to changes in orientation and platform and between and constant offsets between AG’s and SG have not been reduced for this display

Table 2 Drift model for the SG and estimated parameters
Fig. 4
figure 4

Seismic surface wave train after an earthquake in the Philippines. The FG5X-233 picked up the vertical accelerations at 1 \(\mu \)m/s\(^2\) RMS difference with respect to the broadband seismometer at Onsala, station ONA of the Swedish National Seismic Network. The seismometer’s calibration factor had to be honed a little, and the time stamps of the AG series were adjusted to yield maximum correlation at zero lag, owing credit to the seismometer’s GPS-controlled clock

The observation model (all time-dependent terms are weighted with the uncertainties of \(g_a\)) is stated as

$$\begin{aligned}&\beta _c (g_{s,i} - g_{d,i}) + \beta _m b_{m,i} + \beta _{r,m} r_{m,i} + \beta _p p_i \nonumber \\&\quad + \beta _s \frac{t_i-T_s}{\Delta _s} + \beta _0 + \beta _G (t_i-t_0) \nonumber \\&\quad = g_{a,i} + \epsilon _i,\;\;1\le i\le N_{drops}, \end{aligned}$$
(1)

where t denotes time; the \(\beta \)’s are the parameters to be estimated; i enumerates the drop measurements, \(\beta _c\) is the SG’s scale factor; \(g_s\) and \(g_d\) the SG readings and drift, respectively; \(g_a\) are the AG’s drop records reduced for the meter bias applying to the respective campaigns; \(b_m\) is unity if instrument m is observing (creating one matrix column for FG5(X)-220 as FG5-233 is chosen as the reference instrument); \(r_{m}\) is unity if instrument m is operated at south azimuth else zero (two matrix columns in our case); p is a platform tie (three columns as platform AA is the reference); the slopes of the AG trends are given by the time-dependent terms factored by \(\beta _s\) (one column for each of the 81 setups) with \(\Delta _s\) the duration of the setup and \(T_s\) it’s central time; \(beta_0\) is a mean offset, a nuisance term; and finally, the secular rate is factored by parameter \(\beta _G\), G short for GIA. In campaigns with significant microseismic noise levels, acceleration derived from the broadband seismometer is added to \(g_s\) with an amplitude factor and time shift to fit the variations and time-stamp offsets of \(g_a\) in short sections, typically four to five, of the setups; \(\beta _m\), \(\beta _{r,m}\) and \(\beta _p\), respectively, factor the signal that the first subscript indicates; and finally, \(\epsilon _i\) is the error of the model at each reading i. We introduced one more AG term, \(\beta _2\) for observations in June/July 2009. This campaign with FG5-233 was conducted soon after establishing the SG mid-month of June. In order to not affect \(\beta _G\), enhanced by the leverage that the first campaign would have, this campaign could thus at least contribute to the SG’s scale factor; and eventual trends in the setups could be examined.

Fig. 5
figure 5

Frame a shows the determination of AG time-stamp corrections using the seismometer. The cross-covariance between the two instruments’ time series has been calculated from 2-h segments (vertical axis). Frames b and c show a scatter plot before and after regression, respectively, determining a coefficient for each segment. The colour of the dots change from blue to red as campaign time passes along. Note that the correlation between acceleration measured with the AG and with the seismometer is negative. More on that in the main text, the reduced drop series loses 50% of its original RMS scatter as indicated by the error bar’s vertical line (the horizontal shows the repeatability of measuring seismometric acceleration)

Fig. 6
figure 6

Adding constant biases to the DoE’s from the ICC’s at a typical AG uncertainty value (20 nm/s\(^2\)) in every one of the 16 campaigns, once with a minus sign (left half of the abscissa), once with a plus (right), key parameters are affected, the secular rate (top frame), the SG’s scale factor (mid; the scale factor in GAIN campaign \(-\,773.82\pm 0.21\) nm/s\(^2\)/V serving as a reference). Weighted RMS of the overall fit (bottom). Of the 16 campaigns, those showing a lower RMS with offset increments added to the DoE’s (coloured bars) are candidates for continued inspection

The SG’s drift \(g_d\) is determined in a long-term tidal analysis involving reduction of environmental effects in the least-squares inversion (“extended tidal analysis”, see SR19 for details) and fixed subtraction of the nodal tide using a zonal gravity factor (Dehant et al. 1999).

The SG drift coefficient can alternately be solved (introducing an additional parameter \(\beta _d\)) or applied fixed as shown in (1). The SG scale factor can be fixed to the GAIN campaign’s result. These two fixed-parameter options will be applied when we inject meter biases in Sect. 3 to obtain their bearing upon key parameters.

However, not all of the time series in regression are free from an effective long-term linear slope. We decompose the SG readings

$$\begin{aligned} \begin{aligned} g_s(t)&= d(t) + \tilde{h}(t) + \zeta f(t) \\&= d_1\cdot t + \tilde{d}(t) + \tilde{h}(t) + \zeta \cdot (f_1\cdot t) + \zeta \tilde{f}(t), \end{aligned} \end{aligned}$$
(2)

where f contains a linear trend \(f_1\,t\), the tildes denote variables free from linear trends (and thus separable in inversion) and \(\zeta \) is the least-squares adjusted admittance of f, the linear rate of which, \(f_1\), biases the drift’s \(d_1\); they cannot be separated in inversion as they would create a singularity. Thus, the supposed instrumental d picks up the secular term \(\zeta \,f_1\,t\) of physical origin; our rate bias \(g_b\,t\) must be subtracted from \(d_1\). In practice, we sum up all such linear trends in the J columns of the design matrix of the SG analysis and multiply them with the corresponding admittances of the least-squares solution \(\zeta _j, 1\le j\le J\), in order to not miss any of these entrances through the back door.

In general, the SG’s linear drift rate \(d_1\) includes an instrumental and an inseparable physical part. Let us denote the estimated drift function by \(\bar{g}_d\). Removing the linear slope from \(g_s\), the AG’s \(g_a\) in (1) determines \(\beta _G\) alone. The virtue of \(g_s(t)-\bar{g}_d(t)+g_b\,t\) in (1) is in providing a multi-year zero-mean component of observed gravity variations, slope-free except for the slope bias, ready to replace the standard reduction of \(g_a\) on the basis of models. And \(g_b\,t\) will tend to cancel, for example, a secular trend in a local ground water effect picked up by the AG’s—how well it gets cancelled will depend on the accuracy of the estimated admittance coefficient or—more advanced: transfer spectrum—of the hydrological model. This and all other models appear instead at the analysis stage of \(g_s\), much richer and more coherent in its data base and with more favourable statistics with respect to deterministic and stochastic processes. The long-term residual RMS of the SG analysis at 5–7 nm/s\(^2\) around a constant mean is our strongest argument.

A weak part remains in (1), the stability of the AG measurements over the course of years. The aim of the ICC’s is to maintain a gravity reference system by determining the offsets of individual instruments so that the reference can be ported and deviations be traced. Questions arise like, if a significant change occurs from one ICC to the next when 2 years has lapsed, at what time between did it happen, before or after a field campaign? One can try either alternative and notice which assumption returns a more consistent result similar to Olsson et al. (2015b). Or as we will attempt in the last section, to obtain partial derivatives with respect to meter bias variations, which could be the basis for a further step of regression were it not for the dependence on a linear rate \(\beta _G\) to be known a priori. As much as we cannot resolve this ambiguity, we can at least illuminate the ICC effort as being crucial for the success of secular gravity rate projects and demonstrate a need for improvement.

Investigations from other points in Sweden observed annually or semi-annually with FG5 (Olsson et al. 2015b; Engfeldt et al. 2019) and investigations at Herstmonceaux in England observed once every week (Smith 2018) have proved that the gravity trends fit much better when using the meter biases from the ICC’s. If using the latter only for FG5-233 and not for FG5-220 (as in Olsson et al. 2019), the gravity value changes by 15 nm/s\(^2\) for Onsala AA, the most frequently observed point at Onsala (Engfeldt 2019). That is why ICC biases are used in the new Swedish gravity reference frame, RG 2000, for the used absolute gravimeters (Engfeldt 2019).

In regard to references and citations concerning the analysis methods and strategies, we refer to SR19. This article will devote itself to apply them to the specific case in the first hand, to discuss the matter in its own frame and to abstain from comparison with studies employing more orthodox methods at sites elsewhere. We think there is sufficient in content for the reader’s takeaway. New in this report is an account of the rate bias at the stage of determining the drift function of the SG. This spurious signal derives from the finite linear rate that is in principle present in all, yet mostly pronounced in environmental gravity time series and long-term tidal components. We elaborate on the drift bias in a subsection below.

The strength of the use of unbroken SG data to jointly evaluate all AG campaigns is, besides the ability to search for trends during a setup, that all environmental effects and the tidal variation of gravity can be considered equal as the AG platforms are at close range (the oft-visited AA and AC at only two metre).

1.4 SG scale factor

We do not expect the SG scale factor to change over time. There are good reasons, one being that the factor determined with quantum gravimeter GAIN (Freier et al. 2016) is almost compatible with the average of the set of FG5-campaigns; in those—second reason—setup slopes are not robust enough for estimation. If such slopes do occur in the AG, they regularly correlate with the tide signal, so that in our experience the spread of the scale factor over time is greater than the formal error (and the normalized error, unity \(\chi ^2/N\), too). In order to shed light on a potentially hidden problem, we investigated the residual of the multi-campaign adjustment, with setup slopes being adjusted, campaign by campaign by regressing the drift-reduced SG signal, DC levels removed (Fig. 2). The individual campaigns resolve the scale factor with 0.4–5.2 nm/s\(^2\)/V uncertainty, while coefficients solved in the campaign-by-campaign inspection often make a close match. AG trends and their slopes are shown in Fig. 3, showing FG5X-220 and FG5-233 drop data during the parallel campaign in May 2014.

Did the scale factor change, a test on the residual, which we shall detail below, would detect it. The setup slopes we investigate differ in character with respect to Imanishi et al. (2002) in so far as the latter study resolves it from basically daily averages, while the AG’s trending behaviour is inferred from a regression of tides and atmospheric loading, assuming a noise model that yields uncorrelated samples after double-differencing in time. They determine the scale factor from the ratios of the tidal coefficients using a parallel recording of their SG T011 and compare both instruments’ trends. Generally, our SG’s residual, reduced by the linear drift segments we fit (Table 2), appears much smaller in range (except in situations of fast changes of air pressure).

Fig. 7
figure 7

Estimating the uncertainty of the rate bias. Shown are periodograms (unless indicated otherwise) for the reason that here important long-period end is inaccessible to windowed stacking methods like Bartlett’s procedure. The spectrum of the sum of environmental effects determined in least-squares regression is shown in dark blue, the residual of the adjustment in the extended analysis in black, the residual where only the local barometer and the tide gauge are used (no hydrology, no non-tidal ocean loading) in green. The extended analysis’ residual is modelled with a prediction error filter of length 300 (using the MEM method of Burg), shown in yellow. With this filter and a random Gaussian deviate generator, 200 artificial series are produced in which linear trends are found. A histogram of the rates of these trends is shown in the inset. The spectrum of one of the artificial series is shown in grey. The MEM spectrum shown is yet unscaled; in the simulation the reproduction of the residual’s RMS is duly warranted

The results do not suggest a systematic variation of the SG’s scale factor over time. However, what appears to happen in those cases where scale factors deviate is that significant setup slopes exist in those campaigns, sometimes only within one of the setups. Whence, we propose that inferred variations of the scale factors reported in a range of studies may be biased due to setup slopes that were neglected, or better, are beyond what one can achieve when campaigns are treated in isolation. Yet, significant scale factor variations may arise, e.g. after an SG’s cryogenic system suffered discontinuity, like reported in Meurers (2018). Less surprising are the large excursions in those single-campaign solutions when the number of drops is small (campaign 8, 201106c), signalled by large uncertainties of the scale factors (cf Fig. 2).

The conflation method of Crossley et al. (2018), also that treating each campaign without entrusting an SG to tie them together over long stretches of time, would not be able to detect such slopes; however, it would attribute enlarged uncertainties that encompass both sources of deviation, the mentioned instrumental contribution—which is a systematic error—and actual variations of gravity unable to discriminate. It is the best you can do.

When we said above “almost compatible”, we admit that the two values differ more than their uncertainty. However, what will result from the discussion in Sect. 3 is the sensitivity of the multi-campaign result to apparently significant excess offsets in some campaigns, so that suggesting GAIN v FG5 multi-campaign compatibility is still tentative at this point. Now, had setup slopes of the FG5X-220 during the GAIN campaign in February 2015 been significant, the close fit achieved (FG5: \(-\,774.83\pm 3.0\) nm/s\(^2\)/V, GAIN: \(-\,773.82\pm 0.21\) nm/s\(^2\)/V) would still lack the crispiness of FG5 measurements in low microseismic noise conditions (the record breaking example being 201606a with the FG5-233 at only \(\pm \,0.32\) nm/s\(^2\)/V uncertainty). This notion holds also for the detection of setup slopes in February 2016, so that the alleged close fit would not only be tentative but also weak in resolution.

A revisit of the least-squares solution of the GAIN campaign showed that it was not free from a trend either. Neither sign nor range could be attributed to the SG. We found a slope range of \(-\,1.43 \pm 0.30\) nm/s\(^2\) while the long-term SG slope range amounted to 0.078 nm/s\(^2\) over the 140 h of data used. Parameter correlation between scale factor and slope was \(-\,0.026\). Not all of GAIN’s acquisition turned out useful: after the first week, the group from Humboldt University, Berlin, made adjustments to the ground-vibration isolation control circuit. Lastly, the scale factor was determined from 140 h of GAIN recording with the least amount of breaks in these two weeks. In the segment used were 119 breaks of 11 to 13 samples length, one of 19, and one of 148, 40,011 samples in total. The sampling interval (it is an integration time, not some sort of drop interval) was 12 s. As Freier et al. (2016) showed, microseismic perturbations could effectively be rejected owing to the active floor mounting system.

The large signal that the residual covariance analysis reaps from campaign 16 (red cross in Fig. 2) is worth further exploration for principle reason. So far, we have approximated sagging AG behaviour with straight lines only. The lower-than-usual value of little-g obtained from this particular campaign is not low enough to explain a 35 nm/s\(^2\)/V difference, not by orders of magnitude, as it would imply a 4.5% lower sensitivity on the AG’s part. It appears more plausible to search for higher-order polynomials to describe AG trends during the two setups; yet, it will be difficult to prove as it seems to be a rare condition.

1.5 Teleseismic perturbation, microseismic noise

At the high end of the frequency scale, microseismic noise might cause high drop scatter. At the frequencies of teleseismic waves, the AG proved to be a reliable seismometer, only that the sampling rate is too low to be useful in seismology (see Fig. 4). Surprisingly, while motion at periodicities around 30 s is picked up by the FG5 with high fidelity including the positive sign of the response, the situation reverses in the microseismic band. With careful adjustment of the time stamps issued by the FG5 control program, the maximum absolute value of cross-covariance between seismometer and FG5 implies a negative sign. In campaigns with an elevated microseismic noise level, a significant part of the drop noise can be reduced as shown in Fig. 5. Thus, in the multi-campaign regression system the seismometer’s acceleration series is subtracted from the AG drop sequence. Microseisms are blocked from invading the system via the SG series, owing to low-pass filtering (combined with the compensation of the GGP filter’s group delay) as described in SR19.

In the 9 years of operation, two campaigns gained from applying this procedure: 201502a when the North Atlantic was in typical late winter state of unrest and 201806a when just poor conditions had presented us of a 15 dB noise level upped above what’s typical for a summer season. We have included this experience as it vindicates our method working with drop-stage data; it would not work on sets.

2 Results and discussion

The results from the multi-campaign analysis are given in Table 3. In the subsequent subsections, aspects of our findings will be discussed.

2.1 Remarks on the uncertainty of rates and rate bias

The rate bias is a linear slope gone missing in the estimated drift function at the stage of SG data’ tidal analysis. It originates in the set of environmental effects in regression. If any of these harbour a trend, it is subtracted from the SG’s drift with each of the respective series’ admittance coefficient. This leakage is principally at work at all the signals in regression and trades slope signal with the drift terms. However, uncertainty in the sum of these terms arises only in the case of signals that contain stochastic errors; that is, tides take no part in this. In order to calculate the rate bias uncertainty we need to quantify these stochastic components, which would multiply owing to each series’ inherent difficulty. Instead, we try to argue as follows long a more heuristic path: let us assume the residual of the regression carries the spectral character of the modelling errors, probably dominated by the atmospheric model, itself conveying the largest gravity effect second only to the luni-solar tides. In Fig. 7, the result of a Monte Carlo exercise is shown. It starts with the estimation of a power spectrum using the MEM method to represent the residual. The prediction error filter at the core of the method is applied inversely to generate 200 innovations of this noise. Principally drift-free by construction, the fit of a straight line to the noise returns a set of nonzero rates, the histogram of which is used to determine the standard deviation of the slope rate. We find 0.2 nm/s\(^2\)/year, which would be the upper bracket of the \(1-\sigma \) range. More realistic values would scale with the fraction of signal (essentially: error) amplitude that the environmental series actually contribute. For the dominant term, Atmacs, the fraction is near 90%, so we choose to content ourselves with the uncertainty’s face value of \(\pm \,0.2\) nm/s\(^2\)/year.

Fig. 8
figure 8

Residual averaged over each setup. The official DoE’s, shown as yellow lines, are applied (Jiang et al. 2012; Francis et al. 2013; Francis et al. 2015; Pálinkáš et al. 2017; Falk et al. 2020). The secular rate of \(-\,3.53\) nm/s\(^2\)/year given in the diagram has an uncertainty of \(\pm \,0.32\) nm/s\(^2\)/year in the least-squares fit; the rate bias, however, adds \(\pm \,0.2\) nm/s\(^2\)/year to it

Fig. 9
figure 9

Official DoE’s and estimated campaign offsets assuming known key parameters: secular rate of gravity, SG scale factor, SG drift and rate bias. The SG’s residual during the campaigns is shown in purple. In the account of FG5-220, the variant nailing the meters’ differential offset to zero at the parallel campaign in 2015 is shown in pink, the unconstrained one in red

Table 3 Results of multi-campaign adjustment

A rate bias in the secular trend could emerge for the reason that the offset related to meter orientation changes with the upgrades of the dropping chambers (from FG5-220 to FG5X-220 in 2014, from FG5-233 to FG5X-233 in 2017). Notably, as we estimate instrument-specific orientation effects, they do turn out significantly different; see Table 3. However, Fig. 8 does not fly systematic, urging evidence: the residuals of the two orientations flip from setup to setup in apparently unrelated sequence.

2.2 SG drift solved or prescribed

In Table 3, a large impact can be seen between solved versus prescribed SG drift parameter. It turns out that correlations exists between SG drift, AG series and a secular rate. We can only argue by plausibility that the drift’s ratio of \(\beta _d/\beta _c\) ought to be near unity, \(1.150 \pm 0.004\) being a significant deviation, the SG’s scale factor would significantly differ from the result of the GAIN campaign and the secular rate almost double. More importantly, the excessive drift added with factor 0.15 to the residual of the SG analysis implies a range of 50 nm/s\(^2\) which is an order of magnitude greater than the RMS of the residual. All these notions, added together, encouraged us to prefer to take the extracted SG drift at face value.

2.3 A problem of limited information

Concerning the rates of apparent gravity change, i.e. the SG drift term, a drawback in the current state of affairs arises due to the discontinuation of the ECCO1 ocean series. The plan to replace it by CMEMS, a Copernicus product, had to be dashed, while the makers had to fix problems with the reference system. In order to accommodate the AG campaign of 2018, the trends of the extended SG analysis had to be extrapolated. (Note the drift rate changes between fits X and S in Table 2). The subsequent chapter’s efforts cannot warrant quantitative conclusions before these, and issues discussed next, become clearer.

3 Assessment of campaign offsets

One set of parameters the authors have little hand on is yet worth considering as a target of inquiry, the ICC’s. We devised a series of tests in which we add an incremental offset of 20 nm/s\(^2\) to each of the official values, well inside their uncertainty range, and observe its impact on key parameters, which are the multi-campaign SG scale factor, the secular rate, and the residual RMS. The results of this exercise are shown in Fig. 6. The weighted averages of the unmodified residual, evaluated over each setup, are shown in Fig. 8. For further information, cf. Engfeldt et al. (2019).

The quantities in the plot would provide a set of partial finite differences. From them, a best-fit solution could be constructed that yields a set of estimated campaign offsets. Notion of the danger to end up in a circle sporting nothing else, but internal consistency led us to a more radical approach.

Turning the problem on its head, we assume we know the secular rate from GIA modelling, (\(-\,3.6\) nm/s\(^2\)/year), apply the rate bias from SG analysis (\(-\,0.6\) nm/s\(^2\)/year), take the SG drift function at face value, use the scale factor of the SG from the GAIN campaign and instead solve for campaign offsets along with the other parameters (orientation per instrument, platform ties, setup slopes). What we find (see Fig. 9) is a set that in shape is not unlike the ICC series of meter biases, but tends to exceed its amplitude, and that not only slightly. To blame this on the SG performance is a vain prospect, remembering the analysis with the simplest regression set (tides, local barometer and tide gauge) renders us less than 8 nm/s\(^2\) RMS. The only place the SG residual has a conspicuous excursion of 20 nm/s\(^2\) is in 201009a. Note, however, that the SG series has been subtracted from the AG unaltered, and had the SG been offset for an instrumental reason, the extra offset of the FG5-233 with respect to the campaigns right before and after would have shown an inverted relation (the AG down if we had subtracted a false, positive excursion). A drop downward of FG5X-233’s offset in 2017 and continuing in 2018 has indeed been observed, but not to the extent we find.

At both ends of the timescale in Fig. 9, the inferred offsets might of course be exaggerated due to the coerced secular term. But within a bound of all but implausible GIA rates, \(\pm \,3\) nm/s\(^2\)/year, the range of variation at the ends could not be more than 30 nm/s\(^2\), i.e. dwarfed by the deviations that appear to persist throughout. As a test, we imposed a secular term of zero (but kept the rate bias). The first and the second offsets changed by \(-\,25\) nm/s\(^2\), the last by \(+7\) nm/s\(^2\). For obvious reason, the residual RMS of the fit did not change (as a single campaign would not be sensitive to such a small secular rate).

The decrease in residual RMS of the forced solution with respect to the one with the key parameters estimated is not spectacular: the normalized \(\chi ^2\) decreases from 0.475 in the latter to 0.472. The numbers also tell us that the weights, i.e. the error figures as stated in the drop files, are too large by a factor of 1.4 on average.

4 Conclusions

We used a simultaneous adjustment of all 16 absolute gravity (AG) campaigns at Onsala Space Observatory, Sweden, since June 2009, bypassing a priori models of gravity variation and, instead, using a superconducting gravimeter (SG). We extract its drift in a tide analysis with a range of time series in regression to account for non-tidal, environmental effects and a simple, quasi-deterministic function for the instrumental component of drift, and prediction error filters to whiten the noise. The residual of the AG–SG adjustment indicates remanent offsets during campaigns far exceeding (by a factor of two) the degrees of equivalence determined in international comparison campaigns. The multi-campaign adjustment working with single drop measurements instead of normal points (set averages) afforded us reduction of noise by employing a broadband seismometer in the cases of campaigns with elevated microseismic noise level. In both instruments, practically identical in the development state of FG5-technology, the efficiency of the “Super-spring” seemed limited; we propose that it over-compensated microseismic acceleration in the band with 2–10 s periods. At longer periods, the response to seismic surface waves (\(\simeq \)30 s) appeared almost undistorted.

The simultaneous adjustment of all setups facilitated the resolution of trending behaviour, most probably on the part of the AG’s, on timescales of less than one day. The low residual RMS of the multi-year SG analysis, order of 6 nm/s\(^2\), argues against asserting that SG could episodically drift at order of 1 nm/s\(^2\)/h in typical setup durations of 12 h or more. Whether these trends contain higher orders than a simple linear function of time could not be ruled out; it could explain conspicuous aberrations of the SG’s scale factor mistaken as tide-like signatures in the residual of the regression. We showed systematics of SG scale factor variations with AG trends during some of the setups in a campaign. In the multi-campaign analysis, the linear slopes of these trends could be determined, and the SG scale factor became almost indistinguishable from the GAIN Quantum gravimeter experiment as regards both value and uncertainty. Whence, we conclude that the SG’s scale factor remained constant during the 9 years of its operation. We caution other practitioners of AG–SG comparison to not jump to conclusions as to temporal SG scale factor variations lest its slope and trend correlations are tightly controlled (anomalous, unresolved SG trends on the order of 1 nm/s\(^2\)/h and more would lead to the same problem).

Two sources of systematic error have been found to influence the secular rate of gravity change, (1) portability of the meter biases from the most recent ICC; (2) biases in the determination of the SG’s drift, which conflates the instrumental with the geophysical, the AG’s role being to disentangle them. The secular trend thus determined consists of everything not exhausted in the analysis of 9 years of continuous SG measurements, i.e. not of atmospheric origin (we employ Atmacs, Atmospheric Attraction Computation Service, Klügel and Wziontek 2009, in an advanced scheme) nor due to continental hydrosphere loading (ERAin, Dee et al. 2011) and (ERA5, Hersbach et al. 2019) nor non-tidal ocean loading (ECCO1, Stammer et al. 2003) nor local ocean loading effects by proxy (using local tide gauges) nor leakage from (quasi-)periodic processes like Polar Motion. SR19 mentions the lunar nodal tide as one of the more important sources of a rate bias owing to the extrapolation of the body and loading tide response to its 18.6-year period. Provisionally, we subtract it using a zonal delta factor (Dehant et al. 1999), ignore its loading effect and the associated uncertainties’ entries to the rate bias and admit that more scrutiny could be devoted; yet, the problem ought to boggle every mind in the secular trade.

The test in which all parameters were originally fixed using the official DoE series and the change of the residual RMS was calculated adding constants the DoE’s one by one suggests that improvement could be expected from more precise porting of the ICC results. However, with only one site in the Nordic network of AG stations and two FG5’s a conclusive result as to the secular rate of change of gravity could not have been the goal. With a different choice of the stipulated secular rate a series of meter biases would result oppositely slanted. Improvement with respect to the DoE’s cannot be asserted this way, only their weakness could be assessed.

On the bottom line of this report stands the wish for future instrumental solutions in absolute gravimetry that are more robust as to varying mean levels, better isolated against ground vibrations, and capable for observing more samples, more tightly spaced over longer setups. Promising appears to be the concept of neutral atom interferometry.

The metrological approach of the determination of degree of equivalence of AG guarantees the long-term gravity reference better than 10 \(\mu \)Gal as shown in this investigation. But to tap the full potential of the existing gravity meters, the combined application of gravity reference stations equipped with superconducting gravimeters and the international AG comparisons as requested within the new IGRF will increase the global gravity reference into the few microGal level.