Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7
https://doi.org/10.1186/s40645-024-00608-z

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Impact of the combination and replacementof SLR‑based low‑degree gravity field
coefficients in GRACE solutions
Filip Gałdyn1*   and Krzysztof Sośnica1

Abstract
GRACE and GRACE Follow-On (FO) missions provide time-variable gravity field models of unprecedented qualitythat allow for the hydrological, oceanic, and ice mass change studies on a global scale. However, the very low-degree
coefficients derived from GRACE and GRACE-FO are of inferior quality due to thermal effects acting on satellites
and malfunctioning of the onboard accelerometers. Therefore, ­C20 and ­C30 coefficients describing the Earth’s oblateness and the pear shape of the Earth, respectively, are being replaced by values derived from satellite laser ranging(SLR) in the standard GRACE solutions. This study assesses the impact of the replacement of low-degree gravity
field coefficients in GRACE/GRACE-FO solutions by SLR data on the trend and seasonal signals of ice mass changes
in Greenland and Antarctica. We found that the replacement of the low-degree gravity field coefficients changes
the estimates of trends by 4, 8, and 22 Gt/year in Greenland, West, and East Antarctica, respectively, dependingon the source of SLR coefficients and period for which the coefficients are replaced. In SLR and GRACE solutions, allcoefficients of the same order and the same parity of degrees are strongly correlated. Therefore, replacing only two
selected coefficients may lead to a biased solution. Thus, we propose to combine GRACE with SLR solutions up to adegree and order 10 × 10 to properly consider the sensitivity of each of the techniques to gravity field coefficients,
instead of replacing two coefficients from SLR in GRACE solutions. The combined solution reduces the residual trendof post-glacial rebound from 1.2 to 0.9 Gt/year and from − 57.8 to − 57.0 Gt/year in Scandinavia and South Canada,
respectively, when compared to GRACE/GRACE-FO solutions with the replacement of coefficients. The SLR-GRACE
combination reduces the noise in the GRACE/GRACE-FO solutions by 8%, from 38 to 35 Gt, in the Fennoscandia
region. In the periods when GRACE is at the end of its mission and observations are disrupted, the weights adjust
the contribution from SLR and GRACE based on relative ratio of variances from each techniques. Thus, the combined
solutions are more consistent with independent geophysical models of glacial isostatic adjustment, and the combinations are affected by smaller noise than the standard GRACE solutions and properly account for different sensitivities
of SLR and GRACE techniques to low-degree time-variable gravity field coefficients.
Keywords GRACE, Satellite laser ranging, Ice mass loss, Greenland, Antarctica, GIA*Correspondence:
Filip Gałdyn
filip.galdyn@upwr.edu.pl
1
Wroclaw University of Environmental and Life Sciences, Institute
Geodesy and Geoinfromatics, 25 Norwida St., 50‑375 Wrocław, Poland

1 Introduction
The Gravity Recovery and Climate Experiment (GRACE;
Tapley et al. (2004)) and its successor, GRACE Follow-On(GRACE-FO; Landerer et al. 2020), have brought a significant advancement in our ability to monitor crucial
processes involving mass changes and environmental
transportation in Earth spheres. These missions provide© The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://creativecommons.org/licenses/by/4.0/.

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

monthly data at a spatial resolution of 200–500 km enabling us to gain a better understanding of long-termchanges in Earth’s large ice sheets, land hydrology, andclimate through time-variable gravity field and masschange products. Furthermore, the wealth of data gathered by the GRACE and GRACE-FO missions has revolutionized our understanding of geophysical processes,
particularly in the polar regions (Tapley et al. 2019). By
monitoring mass changes in the hydrosphere, cryosphere, ocean, and solid Earth, we can grasp the intricate
dynamics of these interconnected systems (Luthcke et al.
2013). This knowledge is instrumental in comprehending the long-term implications of climate change and its
impact on sea level rise (Intergovernmental Panel on Climate Change (IPCC) 2022). Ongoing monitoring of theGreenland and Antarctic ice sheets is crucial not only for
its contribution to global sea levels but also for the invaluable insights it offers into the overall health and stabilityof Earth’s polar regions (Otosaka et al. 2023).
Early in the mission, it became evident that the estimated values of the Earth’s oblateness term, C­ 20, derived
from GRACE were not reliable (Loomis et al. 2019).
These estimates exhibited an unexpected periodic signal of approximately 161 days, and their trends differedfrom the values obtained through satellite laser ranging(SLR). Several theories have been put forward to explain
the erroneous 161-day signal, e.g., ocean tide aliasing or
thermal effects affecting the satellite components (Chenet al. 2009). The common current practice is to replace
the GRACE C­ 20 estimates with values obtained through
SLR. This replacement ensures the accurate utilization of
GRACE data products in scientific applications (Chenget al. 2013; Cheng and Ries 2023; Loomis et al. 2020).
The GRACE and GRACE-FO missions are equippedwith an onboard accelerometers to accurately measure
non-gravitational forces such as atmospheric drag andsolar radiation pressure and to distinguish them fromgravitational perturbations. These accelerometers playeda critical role in obtaining precise inter-satellite ranging measurements. Unfortunately, during the last few
months of the GRACE mission, the accelerometer on oneof two satellites was deactivated most of the time to conserve power. To compensate for this loss, an algorithm
was developed for the replacement of the accelerometermeasurements in the form of the data transplant. Thisalgorithm took into account factors such as flight time,orientation, and satellite-specific events. The GRACEFO mission, which builds upon the experiences gained
from the GRACE mission, involves the use of two satellites that rely on precise inter-satellite ranging. After thelaunch of GRACE-FO, it was observed that the accelerometer on GRACE-D had higher levels of noise compared to GRACE-C. Nevertheless, the GRACE-FO data

Page 2 of 20

products have been determined using a similar accelerometer replacement approach as in the final months of
the GRACE mission, with adjustments made for differences in the attitude control system. As shown in Loomis
et al. (2020), estimates of C
­ 30 are less reliable when either
mission is operating in a single accelerometer mode.
Therefore, it is recommended that ­C30 values in GRACE
series should be replaced by SLR-based values sinceAugust 2016, but SLR values of C
­ 30 have also excellent
agreement with all GRACE solutions for March 2012–
July 2016.
Over the years, numerous publications have emerged
attempting to determine the time-variable gravity
field using data from various techniques in addition
to GRACE data. Time-variable gravity field has been
directly determined using observations from SLR (Boninet al. 2018; Löcher and Kusche 2021; Sośnica et al. 2015)
and the GPS-based orbit of low-orbiting satellites, suchas the Swarm mission (Dahle et al. 2020). Combinationsof data from Swarm and GRACE have also been performed (Grombein et al. 2022), by setting the scale factor
based on differences in the sampling rate of the observed
data when stochastically combining other gravity solutions. Combination of data from different analysis centers is also successfully performing (Gauer et al. 2023;
Jäggi et al. 2020). Zhong et al. (2021) proposed monthly
gravity field models combined with SLR and high-low
satellite-to-satellite tracking (HLSST) data. Kang et al.
(2022) applied the combination of SLR gravity solution
up to degree and order 5 plus ­C61 and ­S61 and GRACE
solutions at the normal equation level. By combining thedata from different mission or analysis centers we canenhance our ability to model and predict future changes
in the Greenland and Antarctic ice sheet and its consequences for our planet. Matsuo et al. (2013) used SLRbased gravity field solutions to assess the long-term ice
mass changes in Greenland with the expansion up todegree and order 4, whereas Talpe et al. (2017) used the
SLR solutions up to degree 5 for the combination with
GRACE and DORIS data.
This analysis provides an assessment of the impact of­C20 and ­C30 data replacement versus SLR and GNSS combination for the determination of the time-variable gravity field by leveraging the advantages of two techniques:
the primary GRACE, and the supplementary but equallysignificant technique for low-degree coefficients, SLR.
The latter solutions are expanded up to degree and order
10 which gives in total 117 coefficients that are combined.
For the first time, the consequences of replacement two
SLR-based ­C20 and C­ 30 coefficients from different sources
are assessed for crucial regions of ice mass depletion andpost-glacial rebound, such as Greenland, Antarctica,

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Fennoscandia, and South Canada with a comparison toexternal and independent geophysical models.
Firstly, we focus on presenting the effects of replacing
spherical harmonics ­C20 and ­C30 in areas with significant ice loss, Greenland and Antarctica. We also show
the difference between standard approaches in GRACE
gravity field determination and we test the replacementsof ­C20 and C­ 30 based on own SLR solutions. Due to thefact that estimates of gravity field coefficients are notfully independent in both, SLR and GRACE solutions
because of correlations between coefficients, we propose
a combination instead of simple replacement. In Sect. 3,
we show a combination of SLR-based gravity field solutions up to degree and order 10 × 10 and GRACE up to
60 × 60 and its results, such as the contribution of eachtechnique to the final solution and the values of formal
error, slope, and amplitude for each spherical harmonic.
Subsequently, in Sect. 4, we provide results indicating theimpact of using data combinations on regions with significant ice mass changes. We test our solution by comparing with post-glacial rebound models in Fennoscandia
and South Canada, as well as mean errors over oceans
and we compare this results with hydrological models.
Finally, Sect. 5 contains the conclusion and summary.

2 Data and methodology
The International Association of Geodesy establishedthe International Combination Service for Time-variable
Gravity Fields (COST-G) dedicated to the combinationof monthly global gravity field models generated by individual GRACE analysis centers. Based on the combined
GRACE and GRACE Follow-On solutions provided by
the COST-G (Meyer et al. 2020a; Meyer et al. 2020; Peter

Page 3 of 20

et al. 2022) up to degree and order 60, we determined
the trend and amplitude of ice mass changes in Greenland and the Eastern and Western Antarctic regions. The
COST-G solution is a combination of monthly global
gravity field models provided by the NASA Jet Propulsion Laboratory (JPL 2019), University Of Texas Center
For Space Research (UTCSR 2018), the German ResearchCentre for Geosciences (GFZ, (Dahle et al. 2018)), andby several other analysis centers, e.g., the AstronomicalInstitute of the University of Bern (AIUB, (Lasser et al.
2020)), the Centre National d’Etudes Spatiales/Groupede Recherche de Géodésie Spatiale (GRGS, (Lemoine2023)), the Institute of Geodesy at Graz Universityof Technology (IfG, (Mayer-Gürr et al. 2018)), the Leibniz Universitxiät Hannover (LUH, (Koch et al. 2023)), theChinese Academy of Sciences (IGG, (Wang et al. 2015)),Huazhong University of Science and Technology (HUST,(Zhou et al. 2016)), Tongji University Shanghai (Chenet al. 2018), and Wuhan University (WHU, (Guo et al.
2017)). The COST-G workflow is described in Jäggi et al.
(2023).
We consider several possibilities for improvements of
GRACE-based solutions resulting from imperfections
associated with thermal effects and accelerometer issues
affecting the low-degree terms derived from GRACE andGRACE Follow-On. We investigate the impact of thesechanges on the trend and amplitude of annual signals in
Antarctica and Greenland, i.e., the specific regions withstrong ice mass loss. We replace the ­C20 and ­C30 coefficients according to the recommendations provided inTechnical Note 14 (TN14; Loomis et al. 2020) for theentire period of the ­C20 component and from March2012 for the C­ 30 component. TN14 provides a standard

Table 1 Impact of the replacement of ­C20 and/or ­C30 on trend (T) and amplitude (A) in ice mass change for the regions of Greenland,West Antarctica, and East Antarctica based on COST-G solutions
Greenland
Solution−1

T [Gt ­yr ]West Antarctica
A [Gt]

T [Gt ­yr ]−1

East Antarctica
A [Gt]

T [Gt ­yr−1]

A [Gt]

COST-G raw 60 × 60

− 179.8

70.2− 115.0

8.6

41.9

66.7

CSR raw 60 × 60

− 181.2

64.9− 117.7

17.8

35.4

94.7

C20 from TN14

− 177.7

68.2− 111.3

13.9

51.9

59.2

C20 from SLR S− 177.8

68.6− 111.4

12.1

51.7

59.7

C30 from TN14− 179.6

67.4− 115.3

3.1

41.2

81.3

C30 from SLR S− 180.9

54.1− 112.9

25.1

47.4

151.3
288.7

C30 from SLR S (full)− 180.6

28.0− 113.5

76.5

46.0

C20,C30 from TN14

− 177.6

65.3− 111.6

10.2

51.2

74.0

C20,C30 from SLR S− 178.8

53.1− 109.3

20.6

57.4

142.1

C20,C30 from SLR S (full)− 178.6

27.6− 109.9

72.3

55.7

278.2

The ­C30 replacement starts in March 2012 or is valid for the whole period (full)Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

procedure for all GRACE solutions. Additionally, we considered other replacements, which are listed in Table 1.
We replace the ­C20 and ­C30 components with monthly
models based on SLR derived using the methodology
described by Gałdyn et al. (2024). This research utilizes
SLR observations from two high-orbiting LAGEOS satellites and up to seven low-orbiting satellites, namely
Starlette, Stella, AJISAI, LARES, Larets, BLITS, and Beacon-C (Pearlman et al. 2019), to determine the Earth’s
gravity field coefficients. The authors generate 1-day
normal equations individually for each low-orbiting satellite, and 10-day normal equations for LAGEOS-1/2.
The methodology employed in this study follows theapproach described in Sośnica et al. (2015) and relies on
the updated version of the Bernese GNSS Software (Dachet al. 2015). The authors introduce a more stable methodby splitting the solutions into reduced degree expansions,eliminating unnecessary parameters, and combining the
split normal equations. This approach improves the stability of the results. Notably, by applying new method,
the researchers achieve a four-fold reduction in noiseover the Pacific Ocean compared to the primary results

Page 4 of 20

obtained in this study. The solution is expanded to degreeand order 10 using the normal equations from three consecutive monthly solutions.

3 The effects of spherical harmonics replacementon ice mass loss estimates
The Earth’s gravity field potential can be represented as aseries of spherical harmonic coefficients:
V (r, φ, ) =

GME
R

∞
n=0

R
r

n+1

n

Pnm (sin φ)
m=0

(1)

(Cnm cos m + Snm sin m)
where r, φ,  are the spherical coordinates in the reference frame, GME is the product of the gravitational
constant and the Earth’s mass, R is the semi-major axisof the Earth, Cnm and Snm are the Stokes coefficients ofspherical harmonics of degree n and order m, and Pnm
are the fully normalized Legendre polynomials. Thus, themost commonly used representation for the gravity fieldcomparison is the equivalent water height (EWH), which
represents the temporal (e.g., monthly) variations within

Fig. 1 Mass differences in Greenland and Antarctica due to various replacements of ­C20 and ­C30 values in the GRACE series. The COST-G GRACE
solution was subtracted from each of the models shown. The line represents the 3-month moving average. Note the scale difference for eachregionGałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Page 5 of 20

Fig. 2 Mass differences in Greenland and Antarctica due to various replacements of ­C20 and ­C30 values in the GRACE series. The COST-G GRACE
solution was subtracted from modified solution with ­C20 and ­C30 replaced from TN14, along with SLR S, and GRACE CSR solution. The line represents
the 3-month moving average. Note the scale difference for each region

the Earth’s hydrology. By determining this parameterfor each longitude and latitude, we can study its spatial
distribution. Comparisons are then made in areas characterized by their intense variations in this coefficient.
From the individual spherical harmonic values, the EWHparameter is calculated as:

EWH(φ,) =

∞

n

n=2

m=0

ME  2n + 1 Pnm
4πR2
1 + kn(sin φ)(�Cnm cos m + �Snm sin m)(2)where kn is the Love numbers of degree n, ME is the
Earth’s mass, and Cnm and Snm are temporal variationsof the Stokes coefficients. The  indicates that the values
are after subtracting the signal of the static gravity fieldpart.
Each of the resulting time series from spherical harmonics is transformed into EWH (Wahr et al. 1998)
and removes the static gravity field part by GOCO06s(Kvas et al. 2021). Table 1 compares two GRACE-only
solutions provided by CSR and the combined COST-G

solution, then, COST-G with replaced ­C20 or ­C30 orboth coefficients from TN14 and SLR solution (SLR S)
for the whole period and since March 2012.
We observe that the replacement of spherical harmonics according to TN14 leads to a change in Greenland from − 179.8 to − 177.6 Gt/year and 70.2 to 65.3
Gt in trend and amplitude from April 2002 to October
2021, respectively, without taking the leakage correction into account. However, in the case of replacing
components ­C20 and ­C30 for the same months with our
SLR S, the trend and amplitude are − 178.8 Gt/year and53.1 Gt, respectively, which implies a change in theamplitude of about 20% w.r.t. TN14.
More pronounced changes are observed in the Antarctic region. For East Antarctica, the trend changes from41.9, 51.2, to 57.4 Gt/year in the GRACE-only (COST-G),
TN14 case, and when replaced with SLR S data; which
represents a change of approximately 30% compared withthe original GRACE data. Replacing ­C30 only since March2012 does not change remarkably the trend, but changes
the amplitude of annual signal in Greenland and Antarctica and may lead to inconsistencies because part of theGałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Page 6 of 20

Fig. 3 C20 variation from GRACE COST-G (grey line), GRACE CSR (golden line), TN14 (red line), SLR S (green line)—the top figure; absolute differences
between the TN14 and SLR S—the bottom figure

Fig. 4 C30 variation from GRACE COST-G (grey line), GRACE CSR (golden line), TN14 (red line), SLR S (green line)—the top figure; absolute difference
between the TN14 and SLR S—the bottom figure

time series includes GRACE-only ­C30 whereas the otherpart SLR-only ­C30 values (Table 1). Large differences
between amplitudes in Antarctica occur between two
GRACE solutions, COST-G and CSR, however, COST-G
is characterized by a smaller noise as it combines seven
GRACE solutions, including CSR. Figures 1 and 2 illustrate the temporal changes resulting from the replacement of ­C20 and ­C30 separately and together, respectively.

Figures 1 and 2 show the differences between the raw
GRACE/GRACE-FO COST-G solution and individual
solutions utilizing various variants. We observe that the
introduced changes associated with C­ 20 and C­ 30 have asimilar impact on specific months in East and West Antarctica; however, the magnitude of changes is about 2.5
times larger in East Antarctica.
In each of the depicted regions, changes accumulate
over time, and their fluctuations become increasingly

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

significant, especially after 2014 when the accelerometer
issues became substantial.
Figure 3 shows a time series of C
­ 20 values for the
GRACE/GRACE-FO COST-G solution and CSR RL06,
comparing them with the values from TN14 and SLR S.
The absolute difference between the values from TN14
and SLR S is shown below the main graph in bars. The
root mean square error (RMS) of these differences for theentire period is 4.2 × 10–11, and the highest agreement
between the solutions is observed in the years 2012–2015, where the RMS is 3.7 × 10–11. Two SLR-based
solutions are quite consistent with each other, whereas
GRACE solutions show large discrepancies in the periodwith GRACE accelerometer issues in COST-G (2014–2018) and CSR for the accelerometer issues from GRACE
and GRACE-FO (2014–2022). Figure 4 compares the ­C30
values for the entire period. The colors used for different
solutions in Fig. 4 are the same as for Fig. 3. However,
SLR solutions for C­ 30 appear to be less consistent thanthose for the C­ 20 series, which is related to correlations
between all odd-degree zonal harmonics and a different
number of estimated zonal harmonics in each solution,in the case of this solution ­C30, ­C50, ­C70, and ­C90. We cansee that better convergence between SLR and the dataobtained from GRACE and GRACE-FO was achieved
only after the launch of the Laser Relativity Satellite(LARES). However, these data still differ from the values
in TN14, and the RMS of the differences for the entire
range is 7.7 × 10–11. We also observe that the period
from 2018 onward exhibits a lower RMS of differences of
5.7 × 10–11. SLR S include a larger number of estimatedcoefficients, including those, which are correlated with
­C30, i.e., C­ 70, and C­ 90, which are not estimated in TN14.
Moreover, once-per-revolution empirical orbit parameters in along-track are estimated for Starlette, Stella,and Ajisai in SLR S to compensate for the solar radiation

Page 7 of 20

pressure and albedo mismodeling effects. Hence, thesesmall modeling issues map into some spherical harmonic
differences. Therefore, SLR solutions are strongly vulnerable to the modeling applied.
Based on Figs. 3 and 4 and attempts at a direct replacement of C
­ 20, and ­C30, we can conclude that replacing
the entire time series of selected spherical harmonics
in models developed to different degrees may result inthe loss of information for certain components. Consequently, part of the signal potentially contained in the C­ 30
coefficient may be found in other coefficients correlatedwith ­C30. In addition, this approach of simple replacement of two coefficients completely neglects the formal
errors (standard deviations) of the spherical harmoniccoefficient determined in a given month from the leastsquares adjustment. To fully utilize the potential of thegravity field model and to avoid the commission andomission errors, all spherical harmonics determined up
to the least degree and order 10 from SLR should be usedbecause the SLR solutions are still sensitive to time-variable gravity for some coefficients of degree 10 (Sośnica
et al. 2015).

4 Results
4.1 Combination of GRACE COST‑G and SLR solution

The gravity field model combination can be done either
at the solution level or the normal equation (NEQ) level.
The NEQ level combination takes into account the covariances between the estimated parameters and is considered more robust. However, it relies on the assumption ofequivalent information content among all NEQs. In contrast, the solution level combination, particularly the stochastic combination method, focuses on combining thespherical harmonic coefficients of individual solutions
based on their respective standard deviations obtained
from the least-squares adjustment (Grombein et al.

Fig. 5 Median percentage contribution of SLR (left-hand side) and GRACE/GRACE-FO COST-G (right-hand side) for each spherical harmoniccoefficient in the combined SLR-GRACE solution from 2002.04 to 2021.10

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Page 8 of 20

Fig. 6 The average median of the GRACE-based spherical harmonic coefficient contribution for each month of the year up to degree and order 6Fig. 7 The comparison of parameters for GRACE/GRACE-FO (left-hand side), SLR (middle), and combination (right-hand side). The values presented
from top to bottom are formal error, slope, amplitude of the annual signal, and RMS

East Antarctica

8.6
41.9
66.7

Amplitude [Gt]

− 115.0

Trend [Gt y­ r−1]

Trend [Gt y­ r−1]

70.2

Amplitude [Gt]

Amplitude [Gt]

− 179.8

Trend [Gt y­ r−1]

GreenlandWest Antarctica

COST-G
unfiltered
60 × 60

Parameter

Area

94.7

35.4

17.8

− 117.7

64.9− 181.2

74.0

51.2

10.2

− 111.6

65.3− 177.6

CSR unfiltered C20,C30 from TN14
60 × 60

142.1

57.4

20.6

− 109.3

53.1− 178.8

C20,C30 from SLR

138.5

50.5

26.3

− 112.5

42.3− 177.9

COMB

126.1

51.9

18.1

− 112.0

50.2− 178.2

COMB 2

98.5

52.3

5.7

− 111.7

58.7− 178.3

COMB 4

84.3

52.2

1.4

− 111.6

62.4− 178.3

COMB 6

103.5

52.1

8.3

− 111.8

57.5− 178.2

COMB DOW

Table 2 Comparison of the trend and amplitude of the annual signal for the combined SLR-GRACE models in a given area for the period from April 2002 to October 2021
without removing the GIA effectGałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7
Page 9 of 20

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

2022). In the stochastic combination method, the spherical harmonic coefficients are weighted based on theirformal error. Solutions with lower standard deviations,indicating higher reliability, are given greater weight inthe combined result. This ensures that more reliable solutions contribute more significantly to the final estimateof the gravity field, while less reliable solutions have adiminished influence. The combination is defined by:x̂nm = N

1

N


i
i=1 wnm i=1

i
i
wnm
xnm
,

(3)

i is the weight of each spherical harmonic based
where wnm
i is theon the inverse of the square of formal error and xnm
value of each spherical harmonic. The NEQ level combinations require the full normal equations for COST-G
solutions, which are not publicly available, therefore, only
the combination at the solution level can be conducted.
Based on the above equation (Eq. 3), combined spherical harmonics have been derived. Their average percentage contribution over the period from April 2002 to
October 2021 is illustrated in Fig. 5. For the degree andorder of two, the median of the mean usage of SLR data is
91%. For the degrees and orders of four and five, it is 49%
and 39%, respectively. Analyzing the median of the meanvalues derived only from zonal spherical harmonics, weobserve a value of 24% for degree 10 and 44% for degree 6
from SLR. It is also noteworthy that the values of spherical harmonics C­ 44, ­S44, ­C43, and S­ 43 are highly utilized for
combination. Therefore, the SLR technique shall contribute not only to C
­ 20 and C­ 30 because the potential of SLR
goes beyond these coefficients.
Each solution over the entire period has a varied contribution from each technique. The models are definedas monthly datasets, which have different formal errors
depending on the time period. Based on Fig. 5, particularattention should be given to the spherical harmonics todegree and order 6. They have been analyzed month bymonth in Fig. 6. In Fig. 6, we can also observe a tendency
toward a higher contribution of observations based on
GRACE Follow-On data and early years of GRACE. This
indicates that the technology used in the GRACE-FO
mission, which involves additional laser interferometermeasurements, affects the formal errors and significantly
reduces them. Another evidence is the issues related tothe end stage of the GRACE mission connected withproblems with the power supply and accelerometer that
occurred in 2015–2017. During that period, the average
median is 52%, suggesting that, on average, around half of
the information in that period comes from SLR solutions.
However, over the whole period, the average GRACE
contribution is around 76% up to degree and order 6.

Page 10 of 20

The proposed combination takes full advantage of eachgravity recovery technique and reduces formal errors of
low degree and order spherical harmonics. This is illustrated in Fig. 7, where on the left we observe the parameters related to the GRACE/GRACE-FO COST-G data, inthe middle the parameters from the SLR S solution, andon the right the parameters associated with the combination. Figure 7 shows results for the entire data series from
April 2002 to October 2021.
From the top, Fig. 7 illustrates the median of the formal error, the slope indicating the secular trend in thespherical harmonics, the amplitude of the annual signal
corresponding to the seasonal changes described by coefficients, and last one the RMS value which represents thefitting error using the least squares method after removing the seasonal signals and a secular trend. The weakest part in the GRACE solutions, in terms of the formal
errors, are the degree-2 and degree-3 coefficients, as well
as some sectorial coefficients due to near-polar GRACE
orbits. The combination with SLR remarkably improvesformal errors of degree-2 coefficients and slightly
improves degree-3 and sectorial coefficients of degree-3
and 4.
By examining the slopes of individual spherical harmonics, a striking similarity between the GRACE/GRACE-FO and SLR data can be observed. Except
for the C­ 20 component, the other spherical harmonics
remain at a similar level in the combination. Regarding the amplitudes, an inverse relationship is evident. In
the case of SLR data, data exhibit lower amplitudes than
GRACE, which is due to the fact that the derived modelis based on data from three months. This provides aslight smoothing effect on individual components while
reducing noise originating from SLR observations andenabling better stability of SLR-only solution. The RMS
level is remarkably larger in SLR solutions due to a lower
stability of SLR solutions or some non-seasonal signal.
4.2 Impact of GRACE‑SLR combination on polar regions
using different scaling factors

In this section, we utilize the monthly GRACE-SLR combined gravity field solutions and compare them with
GRACE solutions with and without replacing low-degree­C20/C30 coefficients. Specifically, we compare it with the
GRACE/GRACE-FO solutions derived from the COST-Gand CSR RL06 approaches. According to the recommendations of TN14, we also replace the spherical harmoniccoefficients ­C20 and ­C30 in the COST-G data with thecorresponding and from SLR S. It is worth to mentioningthat the analyzed combination is performed using datafrom April 2002 onward when GRACE solutions became
available. The results regarding the basic parameters ofGałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Page 11 of 20

Fig. 8 The ice sheet loss from time-variable gravity field signal up to degree and order 60 of the Greenland (top), West Antarctica (middle) and East
Antarctica (bottom) recovered from GRACE/GRACE-FO COST-G (green), COST-G with TN14 (brown), COST-G with SLR (blue) and combination(orange). The delta sign indicates that the COST-G was subtracted from each solutions

the data series, such as trends and the amplitude of theannual signal, are shown in Table 2.
The SLR solution is based on stacking individual solutions from degree and order (d/o) 2–4 up to d/o 2–10
and generating 3-month solutions. Therefore, the observations are incorporated several time into the solution,
which results in keeping the formal error low. Moreover, the SLR solutions are characterized by increased
residual values in Fig. 7 when compared to GRACE
solutions. Therefore, we propose rescaling the SLR solutions by multiplying the formal errors by two, four, six,and degree-specific scaling factors, and examining whichof the combined solutions generates the lowest noiseover a very large area where the residual gravity signal is
not expected or expected to be exceptionally small. To

validate these combinations, we calculate the RMS overoceans and the trend over the Scandinavia and South
Canada regions, which are subject to the Glacial isostatic
adjustment (GIA) effect and thus, the gravity field solutions can be compared to external GIA models (Caronet al. 2018). Consequently, after reducing combined SLRGRACE models by GIA, the trend in EWH should beclose to zero.
For each of the analyzed regions with significant ice
mass loss, we observe that when the spherical harmoniccoefficients are replaced, the secular trend changes significantly. For example, for the unfiltered COST-G data in
the Western Antarctic region, the trend value is − 115.0
Gt/year, while the values are − 111.6 Gt/year and − 109.3
Gt/year after the replacement with ­C20 and C­ 30 values

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

from TN14 and SLR S, respectively. For the combined
solution (COMB), the value is − 112.5 Gt/year. Thus, thecombination is somehow closer to the GRACE-only solution than the GRACE solutions with data replaced by
TN14 or SLR S.
Large differences can be observed for the annual signal in Greenland. GRACE solutions provide amplitudes
of about 65–70 Gt, replacements lead to a reduction of
the signal to 53–65 Gt, whereas the COMB results inthe smallest amplitudes of annual signal of 42 Gt. An
opposite situation is observed for the West Antarctica,
where the combination COMB increases the amplitude
to 26 Gt from the level of 9 Gt in COST-G solutions.
For both Eastern and Western Antarctica, the amplitudeis approximately twice as large as in the case of replacing the ­C20 and ­C30 values from SLR S or combining in
COMB than the COST-G solution. These relationships
are also visible in Fig. 8. The figure is divided into regions,with each region showing the change in the average ice
mass over the period between April 2002 and OctoberPage 12 of 20

2021. Additionally, each region has a bar plot showing the differential changes between the initial GRACE/GRACE-FO COST-G solution and the final combined
solution COMB. We can observe that the differences predominantly exhibit seasonal variations throughout theanalyzed period. However, negative differential values
predominate in the early years of analysis, which has an
impact on the final trend in data.
Table 2 provides the values for combinations with
reduced weights for the SLR solution, labeled accordingly based on the scaling factor of the formal error,
COMB 2, 4, 6, and degree-specific downweighting(DOW). In the case of the Eastern part of Antarctica,
we obtain trend results at the levels of 51.9, 52.3, and52.2 Gt/year for COMB 2, 4, and 6, respectively. WhenDOW SLR contribution, the amplitudes tend to converge to the GRACE/GRACE-FO COST-G values, as
the contribution of GRACE/GRACE-FO solutions
increases for each region.

Fig. 9 The zonal C­ 20, ­C50, ­C60 and ­C70, sectorial ­C44 and tesseral S­ 42 spherical harmonics variations from SLR, GRACE/GRACE-FO COST-G,
and the GRACE-SLR combination (COMB)

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Page 13 of 20

Fig. 10 The difference between the GRACE/GRACE-FO COST-G and the combination with different weights of SLR solutions

In COMB DOW, the weights depend on the degreeand order of expansion and are determined by themethod used to generate the SLR S solution. Since
this solution is based on three-month data, the formal
errors of all spherical harmonics were multiplied by thesquare root of 3. Then, the stacking of observations andformal errors of spherical harmonics from d/o 4 were
multiplied by the square root of 4, from d/o 4 to 6 were
multiplied by the square root of 3, from d/o 6 to 8 by
the square root of 2, and from d/o 8 to 10, as they areused only once in the solution, were left unchanged.
Table 2 shows the results of this combination as COMB
DOW. For Greenland, COMB DOW trend is − 178.2
Gt/year with an amplitude of 57.5 Gt. For the Westernpart of Antarctica, the trend is − 111.8 Gt/year, and theamplitude is 52.1 Gt. Therefore, the COMB DOW combination yields results that fall between COMB 2 andCOMB 4 in terms of their outcomes.

Fig. 11 The difference between the trend from the GRACE/GRACE-FO COST-G with replaced C­ 20 and ­C30 from TN14 and the GRACE-SLR
combination COMB in EWH up to degree and order 60

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

4.3 Impact of GRACE‑SLR combinations on individual
spherical harmonic coefficients

Figure 9 shows examples of the time series of selected
spherical harmonics for COST-G, SLR S and the combined solution COMB. It reveals that the values of C
­ 20,
in most cases, are taken uniquely from the SLR solution,
confirming the information provided in Fig. 5 related tolarge formal errors in GRACE solutions. For C­ 50 obtained
from SLR, there is a significant underestimation of theamplitude compared to the COST-G solution until 2012
before the launch of LARES. Consequently, the combination attempts to suppress the high amplitudes relative to
SLR and the solution for that period. Figure 9 also showsthe example of C
­ 60, which is also well determined by SLR
but often fixed by analysis centers. Solutions provided
by CSR following (Cheng et al. 2011), which estimatesmodels complete up to degree and order 5 plus C­ 61 and­S61. A similar approach is used for the calculation of data
to TN14 by Loomis et al. (2019). Based on this example, it is possible to conclude that fixing the other coefficients of this order may result in a loss of informationand influence erroneous estimates of ice mass loss. In
addition to zonal coefficients, a strong agreement is also
observed between certain coefficients derived in SLR S
and GRACE/GRACE-FO COST-G. Figure 9 shows also
sectorial and tesseral ­C44 and ­S42 coefficients, which arestrongly influenced by seasonal variations and confirmsthat not only the zonal terms are well represented in the
SLR solutions but also other low-degree coefficients.
These components are necessary when GRACE solutions
are affected by high formal error, especially in the latest
months of GRACE mission. Moreover, the gaps between
GRACE and GRACE-FO can be replaced by SLR interms of the low-degree coefficients. Figure 9 shows theresults from the combination COMB with equal weights.
Successive reduction of the weight for the SLR solution brings each of the harmonics closer to the GRACE/GRACE-FO COST-G solution.

Table 3 Median RMS of EWH in the Pacific Ocean divided into
specific periods and in the full applying a 300 km Gaussian filter
Solution

I [cm]

II [cm]III [cm]

Full [cm]COST-G 60 × 60

2.5

2.8

2.9

2.7

C20,C30 from TN14

2.4

2.7

2.8

2.6

C20,C30 from SLR

2.4

2.7

2.9

2.6

COMB

2.6

3.2

2.9

2.9

COMB 2

2.4

2.9

2.8

2.7

COMB 4

2.3

2.7

2.8

2.6

COMB 6

2.4

2.7

2.8

2.6

COMB DOW

2.3

2.8

2.8

2.6

Page 14 of 20

Figure 10 depicts the sectorial and tesseral coefficients­ 33, ­S22, ­C22, and ­S32, considering combinations with difC
ferent formal error scaling factors, from 2012 onward. A
noticeable tendency can be observed in the means of aslight deviation from the seasonal pattern of SLR S compared to GRACE solutions when using formal errors
without a scaling factor. In such cases, the combined
solution excessively pulls the results toward the SLR solution, disturbing the inherent seasonality comprised in thedata. The appropriate scaling of formal errors becomes a
crucial aspect of proper representing seasonal signals andtrends included in the spherical harmonic coefficients.
Figure 11 depicts the differences between the trend
from the GRACE/GRACE-FO COST-G with replaced­C20 and ­C30 from the TN14 and the SLR-GRACE combination COMB derived using equal weights. Notably,
positive differences are observed near the north pole,indicating the influence of the C­ 20 and C­ 30 values, which
are predominantly derived from SLR data. This signifies
the Earth’s flattening effect and its impact on the trend inpolar regions. The non-zonal spherical harmonics, which
are relatively unrelated to large-scale phenomena andhave regional characteristics, affect the remaining areas,such as oceans and some land areas. Therefore, different
handling of low-degree gravity field coefficients may havean impact on the estimates of secular trends in sea level
rise or trends in land hydrology which are important forthe draughts and flood predictions. Reducing the weightof the SLR solutions reduces the differences between
GRACE-only and the GRACE-SLR combination, but theremaining differences are still visible in long-term trends
with a smaller magnitudes and slightly alter them compared to the solution adjusted according to TN14 ­C20 and­C30.
4.4 Validation of GRACE‑SLR combination using oceans
and GIA model

In the next step, we validate the gravity field model by
analyzing the median RMS of EWH in a west side of
the Pacific Ocean, covering an area of approximately
92 million square km which is half of the total area. As
shown by Chen et al. (2021), the variability observedin the oceans can be considered as noise, providing an
indication of the model’s overall quality, when excluding the secular changes related to eustatic sea level rise.
For the analysis of RMS values over the ocean, we divide
the data into three periods. The first period spans from
April 2002 to the end of 2009 with the best GRACE
data quality, the second period covers 2010 until theend of the GRACE mission, and the third period corresponds to the GRACE-FO mission until the end of SLR S,
which is October 2021. We also filter solutions by applying a 300 km Gaussian filter and then we calculate theGałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Page 15 of 20

Table 4 Trend and long-term standard deviations (STD) of EWH in the Scandinavia and South Canada regions calculated withremoving the GIA effect
Solution—GIA

Scandinavia−1

Trend [Gt ­yr ]South Canada
STD [Gt]

STD* [Gt]

Trend [Gt y­ r−1]

STD [Gt]

STD* [Gt]

COST-G 60 × 60

0.5

38.5

20.9

− 61.6

425.8

88.7

C20,C30 from TN14

1.2

37.8

20.4

− 57.8

398.4

70.9

C20,C30 from SLR

0.9

36.5

22.3

− 58.4

398.4

76.6

COMB

0.7

35.0

25.1

− 57.4

386.5

87.2

COMB 2

0.9

35.8

23.6

− 57.0

391.3

81.3

COMB 4

0.9

36.6

22.4

− 57.2

395.2

78.4

COMB 6

1.0

36.9

21.9

− 57.5

397.1

76.8

COMB DOW

0.9

36.4

22.8

− 57.2

394.5

78.2

The STD* was calculated with the removed hydrologically LSDM model, trend and the annual and semiannual signal

RMS value. Table 3 shows that the combinations with areduced weight of SLR by 4 times demonstrate comparable median noise levels over the oceans as the solution
with ­C20 and ­C30 replaced according to TN14. For theentire period, this median noise level amounts to 2.6 cm.
However, the noise in COMB 4 is slightly smaller in thefirst period when compared to COST-G and COST-G
with TN14 solutions. The combination COMB incorporating GRACE/GRACE-FO and SLR with the same
weight increases the noise levels over the oceans. However, the combination with a fourfold reduction in theweight of the SLR solution or COMB DOW with degreespecific weighting achieve similar or better results to thesolution based on TN14 data replacement. Therefore,
downweighting by a factor of four or degree-specificweighting is necessary to avoid the increased noise in the
combination.
For validation purposes, we also removed the effect of
GIA using a geophysical model (Caron et al. 2018) andcalculated the trends for the Scandinavian region andSouth Canada. For the combined solutions in Scandinavia, we obtained the following residual trends after
correcting for GIA: 0.7, 0.9, 0.9, 1.0, and 0.9 Gt/year for
COMB, COMB 2, COMB 4, COMB 6, and COMB DOW,
respectively. When replacing C­ 20 and C­ 30 with TN14, theresulting value was 1.2 Gt/year, suggesting an inferior
agreement with the GIA model when compared to direct
GRACE-SLR combinations. Values close to zero for the
combinations indicate that they are more consistent
with the geophysical model, whereas the smaller standard deviations (STD) of EWH indicate that combinations
do not introduce additional noise in smaller regions suchas Scandinavia, where the expected seasonal hydrological signals are small. A low value suggests that the data isless noisy and does not contain significant disturbancesfrom the solution. The trend value is lower for both the

Scandinavia region and South Canada in COMB 2 thanin TN14 solutions in terms of absolute values. Referring to the results obtained in Table 3, where COMB 2
yields a result of 2.7 cm and COMB 4 and COMB DOW
yield 2.6 cm, as well as the results from Table 4, we canconclude that solutions COMB 2, COMB 4, and COMB
DOW provide the best combination results. The advantage of COMB 2 is its long-term trend, which is proved
to be more consistent with GIA model than COMB 4,
whereas COMB 4 is characterized by a smaller noise overoceans than COMB 2. The residual trends and STD values are larger in South Canada than Scandinavia becauseof the impact of ice mass depletion and seasonal landhydrology variations in Canada. Moreover, the GIA effectis much better scrutinized in a smaller region of Scandinavia than in South Canada, where some leakage effects
from the north part of Canada can be observed. We alsodecide to subtract, from each of the models, the hydrological signal using the land and surface discharge model(LSDM; Dill 2008). We calculate the standard deviation(STD*) based on these results, after additionally subtracting the annual and semi-annual components from theremaining data signal. In this case, the results do not indicate that the combined model achieves the best outcomesbut rather strengthen the suggestion that the contribution of SLR data in the final solution should be reduced,
as COMB solutions exhibit much higher STD comparedto others, such as COMB 2 and COMB 4. The seculartrends from the COMB solution are shown in Fig. 12,
using a Gaussian filter with a radius of 300 km and with
the GIA effect removed. Figure 12 shows some residual
trends in the Canadian region, whereas the residual gravity field changes in Scandinavia are close to zero.

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Page 16 of 20

Fig. 12 The secular trend in EWH based on the GRACE-SLR combination COMB 2 up to degree and order 60, using a Gaussian filter of 300 km
and removing the GIA effect

4.5 Comparison of EWH changes with hydrological models

Now, we will focus on a detailed analysis of the Amazon River basin by calculating the averaged time series.
The Amazon Basin experiences the world’s most significant seasonal variations in water storage due to substantial precipitation. For this purpose, we will compare the
combination and replacement with the land and surface discharge model (LSDM, Dill 2008) and the global
land water storage dataset assimilating GRACE andGRACE-FO data (GLWS2.0; Gerdener et al. 2023). Figure 13 shows the data series from 2003 to 2020, which
is a consistent time frame for all models. The GLWS 2.0hydrological model, supported by GRACE/GRACE-FO

mission data due to its specificity and higher spatial resolution, exhibits lower amplitudes than the other models
being compared within the Amazon River basin. The
LSDM model, on the other hand, was developed independently of GRACE/GRACE-FO products and is basedon the existing Simplified Land Surface Scheme (SLS)and Hydrological Discharge Model (HDM). In LSDM theamplitude remains relatively consistent across this model,
the combined approach, and simple replacement. However, there is a slight tendency to overestimate the amplitude during the wet years of 2008–2010 in this region(Gloor et al. 2013). Additionally, a minor phase shift isobserved at the start of the analyzed period, spanningGałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Page 17 of 20

Fig. 13 Time series of EWH from COMB 2, COST-G with replaced ­C20 and ­C30 from TN14, land and surface discharge model (LSDM) and GLWS 2.0
for Amazon and Brahmaputra river basinsfrom 2003 to 2010. As we move further into the analysis,
these distinctions become less pronounced and eventually indiscernible. Similarly, we examine a time series
originating from the Brahmaputra River basin and conclude that these areas are not significantly burdened bydifferences in determined low-degree spherical harmonics. We do not observe significant distinctions between
our proposed combination and the replacement of ­C20
and ­C30 with TN14.

5 Discussion and conclusions
In this article, we analyzed the impact of replacing spherical harmonics C­ 20 and ­C30 on regions with significant ice
mass changes using data from TN14 and our own SLR
solution. We examined the effect of replacing both coefficients together as well as each coefficient individually.
We observed that by replacing ­C20 and ­C30 from SLR S in
the same periods as the TN14 data, we obtained similar
trends but with larger amplitudes of seasonal signals inGreenland and Antarctica. This was particularly evident
in the East Antarctic region. The influence of different­C30 values, which add significant seasonality to the data,
contributed to this effect. The secular trends and seasonal signal in low-degree gravity field coefficients are
similar in SLR and GRACE solutions. However, gravity
field coefficients are correlated in both technique-specificsolution, therefore, the gravity signal can be distributed
over different coefficients describing similar mass distribution. Therefore, we performed a combination basedon the formal errors of the spherical harmonics. For the
combination at the solution level, we used time-variablegravity field data from SLR S and GRACE/GRACE-FO
COST-G data.
For SLR solutions, we used the expansion up to degreeand order 10 × 10 obtained by splitting and re-stacking
the normal equations. The SLR S solutions were combined with GRACE solutions using formal errors for theweights in the combination. According to the approach,in months with high formal errors in the GRACE solution with issues with power supply or accelerometer
problems, SLR data were used to a greater extent. The
SLR data were used up to d/o 10, which corresponds tothe sensitivity of SLR solutions to time-variable gravity
field changes for which the SLR-only solutions are still
stable. The calculated average contribution showed thateven spherical harmonic C­ 90 can slightly benefit from
SLR data when combining them on a 1-to-1 basis. The
median contribution of SLR data was 91% for d/o 2, 49%
for d/o 4, and 39% for d/o 5. Monthly analysis showedthat during the period 2015–2017, SLR data accounted
for approximately 50% of the contribution used for combinations up to d/o 6.

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

We calculated trends in polar regions based on combined GRACE-SLR solutions. The unfiltered COST-G
data showed a secular trend value of − 115.0 Gt/yearin the Western Antarctic region, while the trend value
replaced with C­ 20 and ­C30 values from TN14 was − 111.6
Gt/year. For the combined solution COMB, the value wasin the middle equaling to − 112.5 Gt/year. To assess thequality of the solution, we calculated RMS values overparts of the Pacific Ocean, which indicate noise and thus
should be minimized. The computed values showed thatthe combination slightly increased the noise, with a valueof 2.9 cm for the combination and 2.6 cm for the replacement of ­C20 and ­C30 with TN14. Based on this knowledge, the weights of the SLR contribution were reduced
by increasing the scaling factor of SLR formal errors.
We applied scaling factors of 2, 4, and 6, and also created degree-specific weighting based on normal equation
stacking during the determination of SLR-based models,i.e., COMB DOW. The combination COMB DOW performed comparably to the combination with scaling factors of 2 and 4, which resulted in the superior statistics.
The noise in the oceans decreased in these cases andreached the level of replacing data based in TN14 or evenslightly better. The combinations also introduced changes
in areas related to ice melting. In the case of the easternpart of Antarctica, we obtained trend results of 51.9, 52.3,
and 52.2 Gt/year for COMB 2, 4, and 6, respectively.
Finally, we have determined the values of these changes in
the Scandinavian region and South Canada after subtracting the GIA effect. The expected residual gravity field signal should be close to zero provided that the GIA model iscorrect and no secular hydrological changes occur in theseregions, which is not fully true for Canada. For the combinations performed in the Scandinavian region, we obtained
the residual trend of 0.7, 0.9, 0.9, 1.0, and 0.9 Gt/year for
COMB, COMB 2, COMB 4, COMB 6, and COMB DOW,
respectively. When we replaced C­ 20 and ­C30 with TN14, thisvalue amounted to 1.2 Gt/year. The STD value, which represents long-term noise in Scandinavia and South Canada,
is the lowest for COMB 2. This suggests, based on both
Tables 3 and 4, that the COMB 2 solution provides stable
results and relies on a significant contribution from SLR
models for C­ 20 with non-negligible contribution to otherlow-degree gravity field coefficients from SLR. Therefore,
the combination of SLR and GRACE data provides moreconsistent results than replacing two coefficients in GRACE
series by SLR-derived values for arbitrarily selected periods.
Acknowledgements
We would like to acknowledge the International Laser Ranging Service (ILRS)(Pearlman et al. 2019) for providing SLR data, Astronomical Institute, Universityof Bern (AIUB) for providing SLR solution, and COST-G for providing GRACE
combined solutions (Meyer et al. 2023).

Page 18 of 20

Author contributions
KS and FG proposed the topic, conceived and designed the study. FG carriedout the experimental study. KS and FG writing review and editing manuscript.
KS acquired funding. All authors read and approved the final manuscript.
Funding
This work was supported by the National Science Centre, Poland (NCN), Grant
No. UMO-2021/42/E/ST10/00020. The article is part of a PhD dissertation titled
Time-variable Earth’s gravity field derived using SLR, GNSS, and GRACE data,prepared during Doctoral School at the Wrocław University of Environmental
and Life Sciences. The APC is financed by Wrocław University of Environmental
and Life Sciences.
Availability of data and materials
The combined GRACE-SLR solutions COMB 2 will be published by the IAG
Service International Center for Global Gravity Field Models (ICGEM).

Declarations
Competing interests
The authors declare that they have no competing interests.
Received: 25 July 2023 Accepted: 2 February 2024

References
Bonin JA, Chambers DP, Cheng M (2018) Using satellite laser ranging tomeasure ice mass change in Greenland and Antarctica. Cryosphere12(1):71–79. https://​doi.​org/​10.​5194/​tc-​12-​71-​2018
Caron L, Ivins ER, Larour E, Adhikari S, Nilsson J, Blewitt G (2018) GIA model
statistics for GRACE hydrology, cryosphere, and ocean science. Geophys
Res Lett 45(5):2203–2212. https://​doi.​org/​10.​1002/​2017G​L0766​44
Chen JL, Wilson CR, Seo K-W (2009) S2 tide aliasing in GRACE time-variablegravity solutions. J Geodesy 83(7):679–687. https://​doi.​org/​10.​1007/​
s00190-​008-​0282-1
Chen Q, Shen Y, Francis O, Chen W, Zhang X, Hsu H (2018) Tongji-Grace02s andTongji-Grace02k: high-precision static GRACE-only global earth’s gravity
field models derived by refined data processing strategies. J Geophys Res
Solid Earth 123(7):6111–6137. https://​doi.​org/​10.​1029/​2018J​B0156​41
Chen J, Tapley B, Tamisiea ME, Save H, Wilson C, Bettadpur S, Seo KW (2021)
Error assessment of GRACE and GRACE follow-on mass change. J Geophys Res Solid Earth 126(9):1–18. https://​doi.​org/​10.​1029/​2021J​B0221​24
Cheng M, Ries J (2023) C20 and C30 variations from SLR for GRACE/GRACE-FO
science applications. J Geophys Res Solid Earth 128(2):e2022JB025459.
https://​doi.​org/​10.​1029/​2022J​B0254​59
Cheng M, Ries JC, Tapley BD (2011) Variations of the Earth’s figure axis fromsatellite laser ranging and GRACE. J Geophys Res Solid Earth. https://​doi.​
org/​10.​1029/​2010J​B0008​50
Cheng M, Tapley BD, Ries JC (2013) Deceleration in the Earth’s oblateness. J
Geophys Res Solid Earth 118(2):740–747. https://​doi.​org/​10.​1002/​jgrb.​
50058
Dach R, Lutz S, Walser P, Fridez P (eds) (2015). Bernese GNSS Software Version
5.2. User manual. Astronomical Institute, University of Bern. https://​doi.​
org/​10.​7892/​boris
Dahle C, Arnold D, Bezdek A, Doornbos E, Ellmer M, Guo J, Dahle C, Zhang
Y, Jäggi A, Klokocník J, Krauss S, Mao X, Mayer-Gürr T, Meyer U, Sebera
J, Shum CK, Zhang C (2020) Description of the multi-approach gravity
field models from Swarm GPS data. Earth Syst Sci Data 12(2):1385–1417.
https://​doi.​org/​10.​5194/​essd-​12-​1385-​2020
Dahle C, Flechtner F, Murböck M, Michalak G, Neumayer H, Abrykosov O, Reinhold A, König R (2018) GRACE geopotential GSM coefficients GFZ RL06. V.
6.0. GFZ Data Services. https://​doi.​org/​10.​5880/​GFZ.​GRACE_​06_​GSM
Dill, R. (2008). Hydrological model LSDM for operational Earth rotation andgravity field variations. Sci Tech Rep. https://​doi.​org/​11.​2312/​GFZ.​
b103-​08095

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Gałdyn F, Sośnica K, Zajdel R, Meyer U, Jäggi A (2024) Long-term ice masschanges in Greenland and Antarctica derived from satellite laser ranging.
Remote Sens Environ. https://​doi.​org/​10.​1016/j.​rse.​2024.​113994
Gerdener, H., Kusche, J., Schulze, K. et al. The global land water storage data set
release 2 (GLWS2.0) derived via assimilating GRACE and GRACE-FO data
into a global hydrological model. J Geod 97, 73 (2023). https://​doi.​org/​10.​
1007/​s00190-​023-​01763-9
Gauer L-M, Chanard K, Fleitout L (2023) Data-driven gap filling and spatio-temporal filtering of the GRACE and GRACE-FO records. J Geophys Res Solid
Earth 128(5):e2022JB025561. https://​doi.​org/​10.​1029/​2022J​B0255​61
Gloor M, Brienen RJW, Galbraith D, Feldpausch TR, Schöngart J, Guyot J-L,
Espinoza JC, Lloyd J, Phillips OL (2013) Intensification of the Amazon hydrological cycle over the last two decades. Geophys Res Lett
40(9):1729–1733. https://​doi.​org/​10.​1002/​grl.​50377
Grombein T, Lasser M, Arnold D, Meyer U, Jäggi A (2022) Determination andcombination of monthly gravity field time series from kinematic orbitsof GRACE, GRACE-FO and Swarm BT—geodesy for a sustainable earth(Freymueller JT, Sánchez L (eds.), pp 191–201). Springer International
Publishing. https://​doi.​org/​10.​1007/​1345_​2022_​163
Guo X, Zhao Q, Ditmar P, Liu J (2017) A new time-series of GRACE monthly
gravity field solutions obtained by accounting for the colored noise in
the K-Band range-rate measurements. https://​doi.​org/​10.​5880/​icgem.​2017.​004
Intergovernmental Panel on Climate Change (IPCC) (2022) The ocean andcryosphere in a changing climate: special report of the intergovernmental panel on climate change. Cambridge University Press. https://​doi.​org/​
10.​1017/​97810​09157​964
Jäggi A, Meyer U, Lasser M, Jenny B, Lopez T, Flechtner F, Dahle C, Förste C,
Mayer-Gürr T, Kvas A, Lemoine J-M, Bourgogne S, Weigelt M, Groh A(2023) International combination service for time-variable gravity fields(COST-G). In: Freymueller JT, Sánchez L (eds) Beyond 100: the next century in geodesy. Springer International Publishing, pp 57–65
Jäggi A, Meyer U, Lasser M, Jenny B, Lopez T, Flechtner F, Dahle C, Förste C,
Mayer-Gürr T, Kvas A, Lemoine J-M, Bourgogne S, Weigelt M, Groh A(2020) International combination service for time-variable gravity fields(COST-G) (Freymueller JT, Sánchez L (eds), pp 57–65). Springer International Publishing
Kang Z, Ries J, Bettadpur S, Save H, Tapley B (2022) Combined gravity solution
from SLR and GRACE/GRACE-FO. In: International Association of Geodesy
Symposia, pp 1–7. Springer Berlin Heidelberg. https://​doi.​org/​10.​1007/​
1345_​2022_​148
Koch I, Flury J, Naeimi M, Shabanloui A (2023) LUH-GRACE2018: a new timeseries of monthly gravity field solutions from GRACE. In: Freymueller
JT, Sánchez L (eds) Beyond 100: the next century in geodesy. Springer
International Publishing, pp 67–75
Kvas A, Brockmann JM, Krauss S, Schubert T, Gruber T, Meyer U, Mayer-Gürr T,
Schuh W-D, Jäggi A, Pail R (2021) GOCO06s—a satellite-only global gravity field model. Earth Syst Sci Data 13(1):99–118. https://​doi.​org/​10.​5194/​
essd-​13-​99-​2021
Landerer FW, Flechtner FM, Save H, Webb FH, Bandikova T, Bertiger WI, Bettadpur SV, Byun SH, Dahle C, Dobslaw H, Fahnestock E, Harvey N, Kang
Z, Kruizinga GLH, Loomis BD, McCullough C, Murböck M, Nagel P, Paik M,
Yuan D-N (2020) Extending the global mass change data record: GRACE
follow-on instrument and science data performance. Geophys Res Lett
47(12):e2020GL088306. https://​doi.​org/​10.​1029/​2020G​L0883​06
Lasser M, Meyer U, Arnold D, Jäggi A (2020) AIUB-GRACE-FO-operational—
operational GRACE Follow-On monthly gravity field solutions. GFZ Data
Services. https://​doi.​org/​10.​5880/​icgem.​2020.​001
Lemoine J-M (2023) RL05 GSM monthly coefficients. https://​doi.​org/​10.​6096/​
5004
Löcher A, Kusche J (2021) A hybrid approach for recovering high-resolution
temporal gravity fields from satellite laser ranging. J Geodesy 95(1):1–15.
https://​doi.​org/​10.​1007/​s00190-​020-​01460-x
Loomis BD, Rachlin KE, Luthcke SB (2019) Improved earth oblateness rate
reveals increased ice sheet losses and mass-driven sea level rise. Geophys
Res Lett 46(12):6910–6917. https://​doi.​org/​10.​1029/​2019G​L0829​29
Loomis B, Rachlin KE, Wiese DN, Landerer FW, Luthcke SB (2020) Replacing
GRACE/GRACE-FO C30 with satellite laser ranging: impacts on antarctic
ice sheet mass change. Geophys Res Lett 47(3):1–7. https://​doi.​org/​10.​1029/​2019G​L0854​88

Page 19 of 20

Luthcke SB, Sabaka TJ, Loomis BD, Arendt AA, McCarthy JJ, Camp J (2013) Antarctica, Greenland and Gulf of Alaska land-ice evolution from an iterated
GRACE global mascon solution. J Glaciol 59(216):613–631. https://​doi.​
org/​10.​3189/​2013J​oG12J​147
Matsuo K, Chao BF, Otsubo T, Heki K (2013) Accelerated ice mass depletion
revealed by low-degree gravity field from satellite laser ranging: Greenland, 1991–2011. Geophys Res Lett 40(17):4662–4667. https://​doi.​org/​10.​
1002/​grl.​50900
Mayer-Gürr T, Behzadpour S, Ellmer M, Kvas A, Klinger B, Strasser S, Zehentner
N (2018) ITSG-Grace2018—monthly, daily and static gravity field solutions from GRACE. GFZ Data Services. https://​doi.​org/​10.​5880/​ICGEM.​
2018.​003
Meyer U, Lasser M, Jaeggi A, Dahle C, Flechtner F, Kvas A, Behzadpour S,
Mayer-Gürr T, Lemoine J-M, Koch I, Flury J, Bourgogne S, Bourgogne S(2020b) International combination service for time-variable gravity fields(COST-G) Monthly GRACE-FO Series. GFZ Data Serv. https://​doi.​org/​10.​5880/​ICGEM.​COST-G.​002
Meyer U, Lasser M, Dahle C, Förste C, Behzadpour S, Koch I, Jäggi A (2023)
Combined monthly GRACE-FO gravity fields for a global gravity-based
groundwater product. Geophys J Int. https://​doi.​org/​10.​1093/​gji/​ggad4​
37
Meyer U, Jaeggi A, Dahle C, Flechtner F, Kvas A, Behzadpour S, Mayer-Gürr T,
Bourgogne S (2020) International combination service for time-variablegravity fields (COST-G) monthly GRACE series. GFZ Data Services. https://​doi.​org/​10.​5880/​ICGEM.​COST-G.​001
NASA Jet Propulsion Laboratory (JPL) (2019) GRACE-FO monthly geopotential
spherical harmonics CSR release 6.0. NASA Physical Oceanography Distributed Active Archive Center. https://​doi.​org/​10.​5067/​GFL20-​MC060
Otosaka IN, Shepherd A, Ivins ER, Schlegel N-J, Amory C, van den Broeke MR,
Horwath M, Joughin I, King MD, Krinner G, Nowicki S, Payne AJ, Rignot E,
Scambos T, Simon KM, Smith BE, Sørensen LS, Velicogna I, Whitehouse PL,
Wouters B (2023) Mass balance of the Greenland and Antarctic ice sheetsfrom 1992 to 2020. Earth Syst Sci Data 15(4):1597–1616. https://​doi.​org/​
10.​5194/​essd-​15-​1597-​2023
Pearlman M, Arnold D, Davis M, Barlier F, Biancale R, Vasiliev V, Ciufolini I,
Paolozzi A, Pavlis E, Sośnica K, Blossfeld M (2019) Laser geodetic satellites:a high-accuracy scientific tool. J Geodesy 93(11):2181–2194. https://​doi.​
org/​10.​1007/​s00190-​019-​01228-y
Peter H, Meyer U, Lasser M, Jäggi A (2022) COST-G gravity field models for
precise orbit determination of low earth orbiting satellites. Adv Space Res
69(12):4155–4168. https://​doi.​org/​10.​1016/j.​asr.​2022.​04.​005
Sośnica K, Jäggi A, Meyer U, Thaller D, Beutler G, Arnold D, Dach R (2015) Timevariable earth’s gravity field from SLR satellites. J Geodesy 89(10):945–960.
https://​doi.​org/​10.​1007/​s00190-​015-​0825-1
Talpe MJ, Nerem RS, Forootan E, Schmidt M, Lemoine FG, Enderlin EM,
Landerer FW (2017) Ice mass change in Greenland and Antarctica
between 1993 and 2013 from satellite gravity measurements. J Geodesy
91(11):1283–1298. https://​doi.​org/​10.​1007/​s00190-​017-​1025-y
Tapley BD, Bettadpur S, Ries JC, Thompson PF, Watkins MM (2004) GRACE
measurements of mass variability in the earth system. Science
305(5683):503–505. https://​doi.​org/​10.​1126/​scien​ce.​10991​92
Tapley BD, Watkins MM, Flechtner F, Reigber C, Bettadpur S, Rodell M, Sasgen
I, Famiglietti JS, Landerer FW, Chambers DP, Reager JT, Gardner AS, Save
H, Ivins ER, Swenson SC, Boening C, Dahle C, Wiese DN, Dobslaw H,
Velicogna I (2019) Contributions of GRACE to understanding climatechange. Nat Climate Change 9(5):358–369. https://​doi.​org/​10.​1038/​
s41558-​019-​0456-2
University Of Texas Center For Space Research (UTCSR) (2018) Grace static fieldgeopotential coefficients csr release 6.0. NASA Physical Oceanography
Distributed Active Archive Center. https://​doi.​org/​10.​5067/​GRGSM-​20C06
Wahr J, Molenaar M, Bryan F (1998) Time variability of the Earth’s gravity field:hydrological and oceanic effects and their possible detection using
GRACE. J Geophys Res Solid Earth 103(12):30205–30229. https://​doi.​org/​
10.​1029/​98JB0​2844
Wang C, Xu H, Zhong M, Feng W (2015) Monthly gravity field recovery fromGRACE orbits and K-band measurements using variational equations
approach. Geodesy Geodyn 6(4):253–260. https://​doi.​org/​10.​1016/j.​geog.​
2015.​05.​010
Zhong L, Sośnica K, Weigelt M, Liu B, Zou X (2021) Time-variable gravity fieldfrom the combination of HLSST and SLR. Remote Sens 13:3491. https://​doi.​org/​10.​3390/​rs131​73491

Gałdyn and Sośnica ﻿Progress in Earth and Planetary Science(2024) 11:7

Zhou H, Luo Z, Zhou Z, Li Q, Zhong B, Hsu H (2016) A new time series of
GRACE monthly gravity field models: HUST-Grace2016. https://​doi.​org/​10.​5880/​ICGEM.​2016.​009

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Page 20 of 20