1 Introduction

Since decades, laser ranging to passive spherical satellites contributes to the monitoring of the Earth’s rotation and the establishing and maintenance of global reference systems. The launch of such satellites, beginning with Starlette in 1975, was indeed motivated by geophysical questions, among them the measuring and understanding of long period variations of the gravity field (Pearlman et al. 2019). Throughout the 1980s and 1990s, a multitude of studies was dedicated to this topic, mostly applying classical perturbation theory and focusing on low-degree zonal harmonics which are best accessible from this type of analysis. The phenomena searched for in the orbit perturbations included solid Earth and ocean tides (Williamson and Marsh 1985), post-glacial rebound (Yoder et al. 1983; Cheng et al. 1989) and mass redistributions in the atmosphere (Chao and Au 1991; Nerem et al. 1993). Also mass variations in the outer core were discussed as causes (Cox and Chao 2002). At the end of the century, such studies could be based on eight satellites and continuous data spanning a period of more than 20 years (Cheng et al. 1997).

SLR lost its uniqueness in this area when in 2002 the twin satellites of the Gravity Recovery and Climate Experiment (GRACE) mission were launched. GRACE was the second dedicated gravity mission utilizing microwave ranging and accelerometers and designed for recovering the time-variable part of the field (Tapley et al. 2004). The standard science product of the mission are monthly snapshots of the gravity field resolved up to spherical harmonics degree 60 and above. GRACE thus established the mass effect of physical processes on and near the Earth’s surface as a new observation type in geosciences. Its results are widely used to close the sea level budget, to monitor groundwater storage or to quantify the ice mass losses in polar regions (Tapley et al. 2019). After 15 years, the GRACE satellites were decommissioned in October 2017 and replaced by the nearly identical twin GRACE Follow-On (GRACE-FO) launched in May 2018. Since the GRACE observations became insufficient months before decommissioning, the transition to GRACE-FO left a gap in the monthly solutions of more than one year.

Despite GRACE, the interest in gravity information from SLR did not disappear, mainly for two reasons: First, the GRACE results quickly gave rise to the idea to extend them backwards with the longer SLR time series. Second, it turned out that the GRACE solutions had deficiencies in the longest wavelengths. In particular, the \(C_{20}\) values exhibited a strong variation with a tide-like period whose causes are still under discussion (Cheng and Ries 2017). It became standard practice, therefore, to replace the \(C_{20}\) in the GRACE solutions by respective estimates from SLR. Such replacement values are provided by the Center of Space Research (CSR) as part of the official GRACE releases, currently as Technical Note 11 (TN11). These values are extracted from a solution complete up to degree and order 5 using five satellites and taking the higher degrees from a static field (Cheng and Ries 2017).

As for TN11, all recent SLR analyses solve for full sets of spherical harmonics which became feasible with the SLR constellation around 1990. The high altitude of the satellites and the low number of stations, however, impose strong limitations to the maximum degree of expansion. Matsuo et al. (2013) confine themselves to degree 4 for their monthly solutions from which they conclude on an accelerated ice melting in Greenland. Bloßfeld et al. (2015) compute with ten satellites monthly fields up to degree and order 6. In Sośnica et al. (2015), the maximum degree is raised to 10 which required applying constraints to the upper degrees. Bonin et al. (2018) use again a CSR solution up to degree 5 to reconstruct the mass history of Greenland and Antarctica with a global inversion approach. As Matsuo et al. (2013), they encounter serious problems to focus properly on their target area due to the low SLR resolution, placing their hopes finally in long-term means from which observation and leakage errors average out.

It is obvious that SLR will never compete with GRACE in resolution regardless of the progress that will be achieved in this technique in future. Any attempt to retrieve more spatial details from SLR will have to introduce prior information for which GRACE is the best source available. An approved instrument to condense this knowledge into an applicable form is the method of principal component analysis (Lorenz 1956; Preisendorfer and Mobley 1988) which allows to disassemble a time series of spatial data into a set of spatial patterns, the empirical orthogonal functions (EOFs), and time series of associated factors, the principal components (PCs). A number of studies have shown that this procedure is well suited to reveal the dominant signals in the GRACE time series (e.g., Schrama et al. (2007); Schmidt et al. (2008); Forootan and Kusche 2012). In fact, these signals are few as the four major EOFs alone explain more than 80% of the total variability of the GRACE time series. As one might expect, the EOFs contributing most can be assigned to geophysical processes such as the secular mass loss at the polar caps or seasonal changes within the hydrological cycle. As a consequence, we would expect that improved knowledge of these modes prior to the GRACE era would be beneficial for improving process models, e.g., via data assimilation (Eicker et al. 2014).

Up to now, only one study has applied EOF analysis to export the knowledge from GRACE to another space technique. In Talpe et al. (2017), EOFs are used to create a GRACE-like time series from low-resolution fields from SLR and tracking data acquired by the Doppler Orbitography and Radiopositioning Integrated by Satellites (DORIS) system. The strategy applied is based on a projection of the low-degree parts of the EOFs onto the SLR/DORIS fields which is done by a weighted least-squares adjustment. The estimated factors are then used to assemble the untruncated EOFs to monthly fields with the full GRACE resolution. The study is again focused on the polar regions and shows in Antarctica a good agreement with the GRACE results. In Greenland, the match with GRACE is less perfect, and the mass balance shows strong variations in the years preceding GRACE that are not confirmed by the validation data used.

It is easy to see that the potential of the EOF representation is not fully exploited by this method. The most apparent drawback consists in fitting the EOFs to previously determined fields whose resolution might or might not be appropriate for the tracking data. The choices made in the first step of this two-step procedure may thus influence the final result in a way that also a full error propagation cannot compensate for. This dependency can be circumvented only by integrating the two steps to one and adjusting the EOFs directly in the processing of the tracking observations. The present study will show that this can be done straightforward by using the EOFs as base functions in the gravity field modelling. In fact, this strategy opens up a variety of solution types since the EOF representation can be combined with the adjustment of individual spherical harmonics. Five of these solution types will be investigated in this study by real data analyses performed with an own software development. The investigations are based on SLR observations to five satellites going back to 1992.

The remainder of the paper is structured as follows: In Sect. 2, the method is developed starting with an analysis of the consistency of SLR and GRACE and an assessment of the spectral sensitivity of the SLR satellites. In Sect. 3, the data and models used in the numerical computations are described, and the resulting time series are compared with those from GRACE and GRACE Follow-On. Finally, in Sect. 4, the results are validated against independent data in selected case studies.

2 Development of the method

The approach presented in this study is entirely based on the assumption that the monthly gravity fields from GRACE are consistent with the SLR observations. There is a broad consensus that both techniques see the same temporal changes of the field, but also differences have been suspected (Bonin et al. 2018). Hence, it appears appropriate to start with a residual analysis that should provide a clear picture of this point.

As throughout this study, we confine the analysis to the five satellites usually considered as the most valuable for gravity field determination: the two Lageos satellites with altitudes around 5800 km, Ajisai at 1500 km and the pair Stella and Starlette from which Stella moves in a circular polar orbit at 810 km and Starlette in an inclined elliptical orbit with an altitude varying from 800 to 1100 km. We analyse the Lageos data in arcs of 10 days, whereas, for the three low-orbiting satellites, the arcs are shortened to 3 days because of the limitations in modelling the atmospheric drag. The orbit fit for the residual analysis is performed by iteratively adjusting the initial state vectors together with daily scale factors for the thermosphere model and some additional parameters. The model environment and the parameter setup are basically the same as used in the following sections for gravity field recovery with the details given in Tables 2 and 3.

Fig. 1
figure 1

Annual means of post-fit range residuals of the five satellites with different gravity field models applied, averaged from arc WRMS values

Full size image

Figure 1 shows the residuals of the five satellites averaged to annual mean values for six different scenarios: In the first run, only the static field GOCO05c (Pail et al. 2016) up to degree and order 120 is applied (black line). In the GRACE time frame, the static field is then replaced up to a degree \({\bar{n}}\) by the corresponding part of the monthly GRACE solutions from ITSG-Grace2018 (Mayer-Gürr et al. 2018) with the \(C_{20}\) values from TN11. The cases considered here are \({\bar{n}}=2\) (red line), \({\bar{n}}=5\) (blue), \({\bar{n}}=10\) (green) and \({\bar{n}}=60\) (orange). For completeness, the figure also shows the residuals when the gravity field is recovered with one of the approaches from this study (black dotted line).

As a first striking result, it can be seen that the static field matches the SLR observations best around 2010 which corresponds to the midterm of the GRACE mission. This is expected since GOCO05c is partly based on GRACE data and depicts, as a consequence, the mean state of the field during the GRACE period. We speculate that another positive impact is related to the low solar activity and thus quiet thermosphere at that time, in particular for Stella which moves in a sun-synchronous orbit. From 2010 onwards, the residuals increase for all satellites, but fastest at low altitudes, which indicates that the true field strongly deviates from the static field. With the aid of the GRACE solutions, this increase is effectively reduced, leading for the Lageos satellites and Ajisai to a rather balanced curve in the period considered. For Starlette and Stella, the SLR residuals are reduced with similar efficiency, but remain increased around the strong solar maximum in 2001 and, to a lesser extent, the maximum in 2014. These two peaks are most likely introduced by the thermosphere model despite the measures to absorb its errors by appropriate parameters. In particular, the first of these maxima will propagate into all gravity field solutions presented in the following.

In summary, Fig. 1 shows quite unambiguously that the monthly gravity fields from GRACE fit well with the SLR observations. Since the residuals attain their minimum at different stages of our test series, it also provides estimates for the spectral sensitivity of the satellites: Obviously, the Lageos satellites experience not more than degree 2 of the time-variable field while Ajisai and Starlette are sensitive up to degree 5. In case of Stella, the GRACE solutions enhance the fit even up to degree 10. This increased impact is certainly related to Stella’s low polar orbit which is strongly influenced by the ice mass losses at high latitudes. Hence, Stella’s high sensitivity might apply mainly to harmonics that are affected by these changes, not to the full set of coefficients.

Despite this caveat, it can be stated that SLR as a whole is sensitive up to degree 10. Together with the low spatial coverage of the observations, this makes it anything but easy to properly invert the signal from SLR into a time series of global gravity field solutions. The common approach would be to represent the time-variable potential in a month i by a truncated spherical harmonic expansion

$$\begin{aligned} V_{i} = \sum _{n=2}^{n_{max}^{SH}} V_{n}^{SH_{i}} \end{aligned}$$
(1)

with \(V_{n}^{SH_{i}}\) being the potential of degree n given at spherical coordinates \(r, \lambda , \vartheta \) by

$$\begin{aligned} V_{n}^{SH_{i}} = \frac{GM}{R}\sum _{m=0}^{n} \left( \frac{R}{r}\right) ^{n+1} \left( c_{nm}^{SH_{i}} \; C_{nm} + s_{nm}^{SH_{i}} \; S_{nm} \right) \end{aligned}$$
(2)

where

$$\begin{aligned} C_{nm} = P_{nm}(\cos \vartheta )&\cos (m \lambda ) \\ S_{nm} = P_{nm}(\cos \vartheta )&\sin (m \lambda ) \, . \end{aligned}$$

Here, GM and R are the gravitational constant and the radius of the Earth, respectively, \(P_{nm}\) the Legendre polynomials of degree n and order m and \(c_{nm}^{SH_{i}}\), \(s_{nm}^{SH_{i}}\) the coefficients of the expansion. These coefficients would then be solved for. The obvious difficulty in this approach is the proper choice of the truncation degree \(n_{max}^{SH}\): Setting it to a moderate value such as 5 will neglect all signal from the degrees above or, even worse, enable it to alias into the lower degrees. A maximum degree of 10, in contrast, is far to be supported by the observations entailing the need to regularize the solution as done in Sośnica et al. (2015).

A way out of this dilemma can be found when SLR is allowed to “learn” from GRACE. As mentioned in Sect. 1, the method of choice for that is a principal component analysis of the GRACE time series. Such analysis is performed by a singular value decomposition (SVD) typically applied to spatial data. However, it is more appropriate here to apply it to a matrix populated columnwise with the potential coefficients of the GRACE solutions; the orthogonal vectors provided by the SVD are then sets of potential coefficients as well, with the same spatial resolution as the GRACE solutions. The number of GRACE solutions being \({\bar{k}}\), the SVD then yields \({\bar{k}}\) such new coefficient sets, the empirical orthogonal functions, and to each of them a time series of factors, the principal components. The potential in the month i in the GRACE time series can thus be represented either, analogously to Eq. (1), by

$$\begin{aligned} V_{i} = \sum _{n=2}^{n_{max}^{GRACE}} V_{n}^{GRACE_{i}} \end{aligned}$$
(3)

where \(V_{n}^{GRACE_{i}}\) is computed with the GRACE solution for month i, or, equivalently, by

$$\begin{aligned} V_{i} = \sum _{k=1}^{{\bar{k}}} \sum _{n=2}^{n_{max}^{EOF}} V_{n}^{EOF_{k}} \; p_{ik} \end{aligned}$$
(4)

with \(V_{n}^{EOF_{k}}\) synthesized with the coefficients of the k-th EOF and \(p_{ik}\) being the principal component associated with month i and EOF k. Since the EOFs inherit the resolution from GRACE, it holds: \(n_{max}^{EOF} = n_{max}^{GRACE}\).

The major advantage of representation (4) stems from the fact that, without losing any spatial resolution, the EOF series can be truncated at a certain point so that only the most significant temporal signals are retained. Provided that the series can be shortened that way to a reasonable length, it can be transferred to other techniques to be used as base functions for the gravity field modelling. Recovering the field from SLR thus can be based instead of (1) on the truncated expansion

$$\begin{aligned} V_{i} = \sum _{k=1}^{k_{max}} \sum _{n=2}^{n_{max}^{EOF}} V_{n}^{EOF_{k}} \; s_{ik} \end{aligned}$$
(5)

where the parameters to be estimated are the \(k_{max}\) scaling factors \(s_{ik}\). These factors have to be adjusted, as any force parameters, in a linearized model using partial derivatives computed by means of the variational equations. The gravity functional applied in this computation is the complete potential defined by the respective EOF, \(V^{EOF_{k}}\), since the \(s_{ik}\) scale en bloc all degrees and orders.

Table 1 Solution types considered in this study (\(k = k_{max}\) in the text)
Full size table

The truncation of the EOF series in (5) may lead to a neglect of signal as the truncation of the spherical harmonics series does in (1). Another shortcoming of the model may appear when operating outside the GRACE time frame since the EOFs describe the time-variable signal as it is observed by GRACE. The point addressed here is the stationarity of this signal which cannot be assumed in general. For mitigating both these effects, it seems appropriate to give the approach more flexibility by including additional parameters. An obvious choice for that would be spherical harmonics as in (1) which can be truncated now at any degree without the risk of signal losses. By merging (1) and (5), we thus end up with the hybrid approach

$$\begin{aligned} V_{i} = \sum _{n=2}^{n_{max}^{SH}} V_{n}^{SH_{i}} + \sum _{k=1}^{k_{max}} \sum _{n=n_{max}^{SH}+1}^{n_{max}^{EOF}} V_{n}^{EOF_{k}} \; s_{ik} \end{aligned}$$
(6)

with a set of unknowns consisting of the coefficients \(c_{nm}^{SH_{i}}, s_{nm}^{SH_{i}}\) and the scale factors \(s_{ik}\). Note that the evaluation of the EOFs now begins at the degree succeeding the spherical harmonics series. Since the adjustment of individual spherical harmonics keeps getting critical at higher degrees, only the cases with \(n_{max}^{SH} \le 5\) are considered here. Together with the pure EOF approach (5), this results in five types of solutions that are specified in Table 1. The choice of the number of EOFs, \(k_{max}\), will be discussed in the following section.

The use of EOFs as base functions in gravity field determination may remind of the concept of lumped coefficients developed in the early days of satellite geodesy to reduce the multitude of spherical harmonics to observable quantities (e.g., King-Hele et al. (1980)). The effect achieved here is, similarly, a dramatic reduction of the parameter space replacing thousands of potential coefficients by a handful of scaling factors. The benefit from that is both a more stable solution and an enormous gain in spatial resolution. It is needless to say that the latter is completely due to constraints so that no single spatial detail can be credited to SLR.

3 Implementation and validation with GRACE and GRACE Follow-On

Table 2 Data and models used in SLR processing
Full size table
Table 3 Parameter set for SLR processing
Full size table

For investigating the performance of the EOF approach, a series of monthly solutions is generated covering completely the period from November 1992 to June 2019. The computation is performed separately for each month by accumulating the normal equations from a total of six 10-days arcs from the two Lageos satellites and thirty 3-days arcs from Ajisai, Stella and Starlette. The input data are the normal points from ILRS screened here by a median-based procedure and weighted at first individually according to the accuracies provided. In a subsequent step, the data weighting is refined by variance component estimation in order to balance the different noise levels of the laser systems. For this, the normal points are grouped by station, aggregating all observations to all satellites within a month, and for each station a variance component is estimated every month, following the iterative procedure outlined in Förstner (1979). Noise levels estimated in this way (not shown here), e.g., clearly reflect improvements in technology in the late 1990s.

As detailed in Tables 2 and 3, the observations are processed in agreement both with current SLR and GRACE standards. The station coordinates are taken from SLRF2014 and corrected for all known temporal displacements including atmospheric loading. The force modelling for the orbit integration is based on GOCO05c refined by short-term variations from the atmosphere and ocean dealiasing product from the GRACE release 06. For the three low-orbiting satellites, the atmospheric drag is computed with NRLMSISE-00 and, as mentioned, adjusted by daily scale factors. The parameter set further includes empirical accelerations for the Lageos satellites, scale factors for the solar radiation pressure and a range bias for each station-satellite pair intended to compensate for station-dependent center-of-mass corrections. The Earth rotation parameters are held fixed as well as the station coordinates except for a few sites exhibiting exceptionally large residuals.

Fig. 2
figure 2

Empirical orthogonal functions from ITSG-Grace2018. Individually scaled view

Full size image

A decisive point in the model setup is, of course, the computation of the EOFs. As the base for them, we choose again the time series from ITSG-Grace2018 after replacing the \(C_{20}\) by the TN11 values. The series is then reduced by GOCO05c and subjected to EOF analysis exactly as it is, without any filtering or smoothing. This strategy is different from that applied in Talpe et al. (2017) and aims to prevent that regional mass balances are biased by the mass redistribution effectuated by such operations. In doing so, it is accepted that the notorious stripe patterns in the GRACE fields propagate into the EOFs (Fig. 2) and from these onwards to the SLR solutions. In an overall view, the SLR solutions are contaminated that way neither more nor less than the GRACE series itself. However, since the stripes are coupled in the EOFs with signal patterns, their temporal distribution is completely different. Though spatial filtering is in general avoided in this study, this fact makes it appropriate when comparing SLR and GRACE directly on grid points.

Table 4 summarizes a first such comparison with GRACE on a global grid intended to assess the relative performance of the different realizations of the EOF approach. The comparison is carried out with the five solution types from Table 1 with the number of EOFs varied from 4 to 12. Because of the global averaging, the deviation from GRACE is expressed by RMS values of pointwise differences which requires, as mentioned, a spatial filtering of the grids. For the sake of comparison, the table includes a solution following the approach in Talpe et al. (2017). This solution, denoted as “EOF post-fit”, applies the two-step procedure described in Sect. 1 estimating in the first step a gravity field of degree 5. The adjustment in the second step is performed again with the unfiltered EOFs; hence, the results must differ from those in the referenced source.

As a first outcome, it is obvious from Table 4 that SLR accepts only a rather limited number of EOFs. In the majority of cases, the difference to GRACE attains its minimum with six EOFs. If the spherical harmonics are fully expanded up to degree 5, this number drops to four. While it remains true that the first few EOFs retain 80% of the temporal signal, one has to realize at this point that 20% of the signal thus are lost for the recovery. This quota may be even higher when focusing on smaller regions. The mass balance recovered then cannot be other than incomplete, regardless of whether or not the full signal is experienced by SLR.

Table 4 Mean monthly RMS of differences between SLR solutions and ITSG-Grace2018 (155 months, April 2002 – August 2016). Differences evaluated on a global \(1^\circ \times 1^\circ \) grid Gauss-filtered with radius 300 km. Values in cm of equivalent water height
Full size table
Table 5 RMS of regionally averaged differences between SLR solutions and ITSG-Grace2018 (155 months, April 2002 – August 2016). Differences evaluated on a \(1^\circ \times 1^\circ \) grid without spatial filtering. Values in cm of equivalent water height
Full size table
Table 6 RMS of regionally averaged differences between SLR solutions and GRACE Follow-On solutions from JPL (9 months, June 2018 – April 2019). All further details as in Table 5
Full size table
Table 7 RMS of regionally averaged differences between SLR solutions using EOFs from different GRACE time series and ITSG-Grace2018. Solution type is S2+6E (S2+8E in the last case). All further details as in Table 5
Full size table

The second lesson from Table 4 seems to be that the SLR data are not reliant on the spherical harmonics: The differences to GRACE are smallest when the solution is based on the EOFs alone and increase with each spherical harmonics degree estimated in addition. This picture changes, however, when leaving the global perspective and focusing on areas with a spatially homogeneous signal. This is done in Table 5 for some important cases using a statistics based on regionally averaged differences. To simplify matters, the number of EOFs is now fixed to the optimum values from the global analysis; the statistics further includes the ensemble mean of the five solutions as well as a typical SLR solution without EOFs, denoted as “S5+0E” following the logic used for the other shortcuts. It is obvious from the table that spherical harmonics up to degree 3 or 4 bring SLR closer to GRACE in Greenland, West Antarctica and the Amazon basin. Leaving out the spherical harmonics is the best choice only for East Antarctica, but also appropriate for the ocean area where it yields the same low difference as the spherical harmonics of degree 2. The best solution in the ocean is, in fact, obtained with the ensemble mean, with an average difference of 1.9 mm. The difference in trend found with the ensemble mean amounts to 0.08 mm/year which is at the error level of current estimates for sea level change (e.g., Uebbing et al. (2019)). Incidentally, it is seen from Table 5 that in the three first mentioned regions the post-fit approach is outperformed by any of the other solutions.

The regional statistics is continued in Table 6 with the first results from GRACE Follow-On. This still small dataset is of great use for the present study since it provides a fixed point beyond the GRACE period that should be met by the SLR solutions. This challenge addresses explicitly the ability of the EOF approach to cope with possible changes in the signal patterns. The table shows that the GRACE-FO solutions are indeed reproduced with similar accuracies as GRACE itself. In particular for Greenland, however, the spherical harmonics are indispensible for that, with the lowest difference obtained with an expansion up to degree 4. Figure 3 illustrates the situation by showing some of the underlying time series. It reveals that the culprit in Greenland is a deceleration of the ice mass loss after 2017. This development which is beyond the expectations from the GRACE perspective is significantly missed by the pure EOF approach but rather well captured with the aid of the spherical harmonics.

All results reported so far may change in detail with other settings for the data processing, e.g., different arc lengths or different choices for the auxiliary parameters that shall absorb model biases. A further impact may come from the GRACE time series used as source for the EOFs for which ITSG-Grace2018 is only one option among others. To assess this impact, the solution S2+6E was recomputed twice using EOFs from the official GRACE solution RL06 from CSRFootnote 1 and the RL06 mascon solution from the NASA Goddard Space Flight Center (GSFC) (Luthcke et al. 2013). The validation on the global grid indicates a slightly increased error in both cases, yielding an RMS of 4.2 cm with CSR RL06 and 4.1 cm with the GSFC mascons while, for ITSG-Grace2018, 3.6 cm are read from Table 4. As shown by the regional statistics given in Table 7, the error increase is restricted to the ocean area in case of CSR RL06, while the solution degrades in Greenland and West Antarctica when based on the GSFC mascons. By contrast, the mascons allow for a remarkably low error in the ocean which might be related to the absence of striping errors which is a feature of this solution type. A further question might then be if this peculiarity of the mascons also allows for submitting more EOFs to SLR. This is not confirmed when using eight instead of six EOFs which leads to no improvement at the global scale and a significant degradation in one of the regions.

It is not surprising that all statistics in the tables are widely underestimated when propagating the errors from the parameter adjustment. As frequently encountered in satellite gravimetry, these formal errors turn out to be too optimistic. By considering GRACE as the reference, they can be raised to a more realistic level by replacing their global mean in the GRACE period by the respective value of the comparison with GRACE. This is done by a simple scaling of the raw errors \({\hat{\sigma }}\) to calibrated ones

$$\begin{aligned} {\tilde{\sigma }} = \frac{RMS(\varDelta _{SLR}^{GRACE})}{RMS(\sigma _{SLR})}{\hat{\sigma }} \end{aligned}$$
(7)

where \(RMS(\varDelta _{SLR}^{GRACE})\) is one of the values in Table 4 and \(RMS(\sigma _{SLR})\) its counterpart computed from the propagated errors at the grid points. These calibrated errors should cover the influence of the observational noise as well as the truncation error discussed above. Due to the scaling with one global factor, this works well in the ocean area where the calibrated errors are \(\sim \)2.5 mm in the mid-1990s and then fall and level off at \(\sim \)1.5 mm until 2006; for the GRACE period, the average found is 1.6 mm which is in good agreement with the RMS of 1.9 mm found in the comparison with GRACE in Table 5 (ensemble mean). At other regions, the calibrated errors develop in parallel over the decades, but remain too optimistic in absolute terms. For Greenland, as an example, the average during GRACE is 2.7 cm while, from Table 5, a value of 6.6 cm would be expected. This mismatch limits the applicability of the calibrated errors. Another limitation comes from the fact that these errors do not cover the impact of changes in the signal patterns. No effort is made here to assess this effect otherwise since occurrence and magnitude of such changes cannot be deduced from the GRACE period as is evident from the Greenland example.

Fig. 3
figure 3

Monthly mass anomalies in Greenland from GRACE, GRACE-FO and SLR solutions

Full size image
Fig. 4
figure 4

Time series of \(C_{20}\) from TN11 and this study

Full size image

It is worth noting that the SLR analysis suggested here not only benefits from GRACE, but also could support GRACE processing by providing alternative replacement values for \(C_{20}\). Such values accrue from any application of the method, either as direct estimates or, when omitting the spherical harmonics, as composites of the \(C_{20}\) values in the EOFs. The resulting series are close to each other and have in common to be assessed as part of a time-variable field with full resolution. This background lets expect more reliable values than the low-degree field which is the source for the TN11 series. As shown in Fig. 4, the \(C_{20}\) from this study deviate indeed systematically from the TN11 values after 2010, falling below them within ten years by nearly the annual amplitude. Such a difference will have a substantial impact on sea level studies and the mass balance for Antarctica. This and further issues related with TN11 were recently studied in detail by Loomis et al. (2019) who finally recommend an own \(C_{20}\) series from SLR applying GRACE-derived monthly fields up to degree 10. Though this series does not fully coincide with ours, we agree with this study that the TN11 series is conceptually outdated.

4 Regional case studies

To obtain a more detailed picture, the validation is now continued by considering regional mass balances in the time domain. The comparison with GRACE and GRACE-FO will be essential here as well. The main focus, however, will be set on the time prior to GRACE for which a validation is, naturally, more challenging. The datasets suited for comparison mainly originate from altimetry or models which have its own uncertainties. Moreover, none of these datasets can be referred directly to the mass anomalies from satellite gravimetry. In many cases, for instance, a model for the glacial isostatic adjustment (GIA) will be required. In the ocean area, the comparison with altimetry also requires the knowledge of the steric effect due to changes in ocean temperature and salinity. Such additional data introduce further errors so that the validation will not always achieve authoritative results.

A further issue arising here is the handling of the degree-1 coefficients which describe the motion of the geocenter. Such coefficients based on the approach developed in Swenson et al. (2008) are provided for GRACE as Technical Note 13 (TN13) in the official GRACE releases, but are not available for the years before. In principle, SLR is able to assess these coefficients on its own. In fact, their estimates from SLR exhibit a considerable scatter and seem to exaggerate the annual amplitude by a factor of 2 as compared to TN13. In the early years, also large interannual variations are present which may be caused by network effects (Pavlis and Kuźmicz-Cieślak 2008). To prevent the mass balances to be corrupted by such artefacts, we refrain completely from using degree 1. Where appropriate, its influence will be discussed based on the degree 1 series for GRACE.

For simplicity, we also refrain from further differentiating between the five solution types. All balances shown below are computed with their ensemble mean which turned out in Tables 5 and 6 as a good compromise in all regions with some limitations in Greenland.

Fig. 5
figure 5

Monthly mass anomalies in the ocean area. a Area means from GRACE, GRACE-FO and SLR. b Post-processed anomalies from GRACE and SLR and reference time series from steric-corrected altimetry. GRACE and SLR reduced by GIA, mean ocean mass from AOD1B restored, annual and semiannual signal removed. All time series smoothed with 400-day boxcar filter. c Time series from (b) fitted to detrended GRACE series

Full size image

4.1 Ocean

Considering the ocean mass from satellite gravimetry can contribute to sea level studies which are typically based on ellipsoidal sea surface heights observed via satellite altimetry since 1993. Expressed as water heights, the gravimetric mass anomalies should equal the altimetric heights when complemented by the steric height changes. Since we aim to validate the gravimetric time series, we rearrange this balance and reduce the altimetric height by the steric effect. The comparison is then carried out between the gravimetric mass on the one hand and what could be called an altimetric mass on the other.

To put the comparison on a broader basis, we use altimetric heights from two different sources: the global mean sea level (GMSL) product from the ESA Climate Change Initiative (CCI; Horwath 2019, personal communication) using ESA and NASA/CNES satellites, and a GMSL series from the University of Colorado (CU; Nerem et al. (2018)) which is based on NASA/CNES satellites alone. The steric effect is taken from the time series from Dieng et al. (2017) which is used in the CCI project and, as an alternative, from an update from Levitus et al. (2012) provided by the National Oceanic and Atmospheric Administration (NOAA). For emphasizing the scatter of these data, both GMSL series are combined with both series for the steric effect. We thus end up with four solutions for the altimetric mass that should represent approximately the bandwidth of all such solutions.

For being comparable with altimetry, the gravimetric time series are created with an ocean mask excluding a 300-km coastal zone. The ocean mass variations from the AOD product are restored as monthly means and GIA is reduced using the model from A et al. (2012). The annual and semiannual signals that are not present in the GMSL products are removed by adjusting corresponding harmonic functions. For Fig. 5b, c, all time series are further smoothed by a boxcar filter with a length of 400 days.

As obvious from Fig. 5b, the processing steps applied bring SLR and GRACE near to the altimetric time series with only a difference in trend of about 0.6 mm/year. This gap can be closed well for GRACE by using the degree-1 coefficients from TN13Footnote 2 which contribute a trend of the same value. Taking a degree-1 series from SLR,Footnote 3 however, reduces the gap only by 0.1 mm/year. This contradiction cannot be solved here; hence, none of these options is applied.

To enable nevertheless a comparison, the focus can be shifted to shorter timescales. For Fig. 5c, the time series are further manipulated in that sense by detrending the GRACE curve and fitting SLR and the steric-corrected altimetric series to GRACE by offsets and trends. With the restrictions coming along with this manipulation, it can be concluded from the figure that the SLR solution is mostly in or at the edge of the envelope spanned by the altimetric series. This envelope is significantly missed only in the period from 2000 to 2003 which we attribute to mismodelled drag as explained in Sect. 2.

Fig. 6
figure 6

Monthly mass anomalies in Greenland from GRACE, GRACE-FO and SLR and reference data. Gravimetric time series smoothed with 400-day boxcar filter

Full size image

4.2 Greenland

As widespread retreat of its peripheral glaciers is reported, Greenland has been the subject of numerous quantitative surveys. For validating our results, we rely in the first place on the recent study Mouginot et al. (2019) which applies the mass budget method comparing glacier discharge data with ice accumulation from regional climate modelling. As a second source, the IMBIE 2012 datasetFootnote 4 is used which provides a shorter time series but is based on a variety of methods including the mass budget approach and GRACE.

While geocenter variations and GIA play a minor role in Greenland, the area is known for signal leaking out in the surrounding ocean. To compensate for this, we fetch back the signal from a 100-km offshore zone following the approach in Baur et al. (2009). The GIA effect is again reduced using the model from A et al. (2012). For Fig. 6, the mass balances are smoothed, as in the ocean case, by a 400-day boxcar filter. It should be recalled that the unsmoothed case was shown in Fig. 3.

As one might expect, Fig. 6 shows an excellent agreement between IMBIE and GRACE, but, in fact, the agreement is good or excellent between all time series until 2012. This result contains the message that SLR confirms that the Greenland ice mass was rather stable until the late 1990s. In the years after, SLR duplicates truly the GRACE curve, in the shown ensemble mean solution with a slightly misled continuation in the gap before GRACE-FO. As pointed out in the previous section, this gap is bridged with greater success by the individual solution using spherical harmonics of degree 4.

Fig. 7
figure 7

Monthly mass anomalies in Antarctica from GRACE, GRACE-FO and SLR and reference data. Top: West Antarctica and Antarctic Peninsula, bottom: East Antarctica. No smoothing applied

Full size image

4.3 Antarctica

For validating the results for Antarctica, we chose as primary source the IMBIE 2018 study (Shepherd et al. 2018) which applies like its predecessor a synthesis of methods. Further reference data are taken from Rignot et al. (2019) authored by the same group as Mouginot et al. (2019) and, like this study, a dedicated representative of the mass budget method.

For a fair comparison with these data, the computation of the gravimetric mass balances is slightly changed as compared to Greenland. To avoid leakage-in from bordering ice shelves, the mass anomalies are computed without leakage correction. For the GIA reduction, the model ICE6G (Peltier et al. 2015) is used. Furthermore, we consider separately West Antarctica including the Antarctic Peninsula and East Antarctica. This separation is motivated by the diversity of the mass signal since the ice melting in Antarctica affects up to now mainly the western coastal areas of the continent. Besides, the two regions raise different problems for isolating the ice mass loss. While the west of Antarctica is barely sensitive to geocenter variations and the choice of the GIA model, both have a considerable impact in East Antarctica.

According to this, the results in Fig. 7 are easy to interpret only for West Antarctica. As in Greenland, we notice an excellent agreement between all time series from gravimetry and reference data. Again, it follows from this that SLR confirms the estimates for the ice mass loss in the 1990s from other techniques. The situation is quite different in East Antarcica, where Rignot et al. (2019) observe a strong mass loss starting around the year 2000. This hypothesis is definitely dismissed by the gravity solutions which follow closely the IMBIE scenario of a slow mass increase. This agreement, however, is quickly disturbed when changing the model for GIA. The trend of 42 Gt/year that can be measured with the GRACE curve would drop, for instance, to -1 Gt/year with the model A et al. (2012) applied. A change of similar order would be effectuated when using the degree-1 coefficients from TN13 which contribute an additional trend of 33 Gt/year. Taking degree 1 from SLR, in contrast, would increase the trend only by 4 Gt/year. The uncertainties in GIA and geocenter motion thus prevent a proper comparison with the reference data.

Fig. 8
figure 8

Monthly mass anomalies in river basins from GRACE, GRACE-FO and SLR and reference data. Annual signal removed

Full size image

4.4 River basins

The low sensitivity of SLR and the low number of EOFs suggest that mass signals are less well retrieved at smaller spatial scales at the level of hydrological basins. This assumption is verified with Fig. 8 showing the interannual mass variations in ten basins. As can be seen, the SLR solutions are still close to GRACE in a couple of cases, i.e., in the Amazon basin and those of the large African and Siberean rivers. The good fit for the Amazon is indeed not surprising since the Amazon leaves a massive signature in three out of six EOFs. Accordingly, also secular events are well represented, such as the low water storage in the strong El Niño seasons 1997/98 and 2015/16.

It is obvious that the same few base functions are less suited to represent secular anomalies in other basins. The EOF approach thus reaches limitations in regions where water storage is heavily influenced by the local climate or anthropogenic interventions. This failure is evident in Fig. 8 in the basins of Ganges, Volga and Danube where significant details in the temporal variations are missed. A particularly poor result is obtained for the Mississippi basin in which SLR detects little more than a global upward trend. A forward computation shows that around 20 EOFs would be needed to bring this result in the vicinity of GRACE. For Danube and Volga, even 40 EOFs would be required.

To give an impression of the signal missing in the time prior to GRACE, Fig. 8 also displays basin averages from the water storage changes reconstructed by Humphrey and Gudmundsson (2019). This reconstruction is based on a statistical model with meteorological forcing whose parameters were determined with GRACE mascon solutions.Footnote 5 For methodological reasons, this approach yields no annual variations and no trends except those introduced by the meteorological input. As obvious in several basins, this absence of trends is a serious drawback of the method. Nevertheless, the variations at shorter timescales are mostly well retrieved in the GRACE period so that these time series can claim some reliability also in the period before.

Unfortunately, there are few opportunities to complete the validation with less biased data. Comparing with time series from hydrological models would provide little clarity since such models are known for having large artefact trends due to biases in forcing data (Springer et al. 2014). A more conclusive, though heavily selective validation could be based on the observed vertical motion of GNSS stations which can be compared with the predicted gravity-induced deformation. In Fig. 9, such comparison is done for three IGS stations for which relatively accurate positions are provided even for the very first years. In all three cases, we notice a good agreement in phase and amplitude of the annual variations, at the stations Brasilia and Bangalore also the interannual changes roughly agree over the whole period. In contrast, the GNSS series for Cape Ferguson suggests a downward movement in the first few years which could point to a gravity change but may be caused by a lowering of the monument as well. Such ambiguities are almost unavoidable when analysing early GNSS data; hence, this type of comparison will not lead to definitive results either.

Fig. 9
figure 9

Vertical motion of GNSS stations and predicted vertical deformation from GRACE, GRACE-FO and SLR. GNSS positions from weekly IGS solutions (second reprocessing). For consistency with gravimetric time series which do not contain a degree-1 part, GNSS series were reduced for load-induced degree-1 variations determined from the IGS solutions by spherical harmonic expansion as in Wu et al. (2002) and Blewitt (2003)

Full size image

5 Conclusions

It was shown in this paper that empirical orthogonal functions from the GRACE solutions can be used as base functions when recovering time-variable gravity field from SLR. This approach can be combined with the estimation of low-degree spherical harmonics which can help to exploit signal not representable by the EOFs.

The advantage of the approach consists in providing the full spatial resolution when considering extended regions with a strong mass signal such as the ocean area, Greenland or Antarctica. It was found in these cases that the results are close to GRACE and mostly fit well with the first results from GRACE Follow-On. For the time prior to GRACE, a good agreement with reference data from altimetry or other techniques was achieved. By this means, present assumptions about the ice mass loss in Greenland and West Antarctica during the 1990s could be confirmed.

As expected, the approach does not succeed in breaking down the limits of the SLR technique. These limits manifest here in the fact that only a small number of EOFs can be reliably determined. Accordingly, the recovered mass signal must be incomplete in areas where significant parts of the signal are stored in higher EOFs. This shortcoming is evident in smaller river basins while large basins such as the Amazon are captured with high accuracy.

Further improvements of the method may be achieved by replacing the EOFs by spatial patterns from more advanced decomposition techniques like rotated EOF analysis or independent component analysis (Forootan and Kusche 2012). Most likely, however, the low sensitivity of SLR will be the limiting factor here as well. A more rewarding topic for future work should be, therefore, to apply the method to satellites at lower altitudes tracked by GNSS or DORIS.