Impact of the combination and replacement of SLR-based low-degree gravity field coefficients in GRACE solutions

GRACE and GRACE Follow-On (FO) missions provide time-variable gravity field models of unprecedented quality that 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, C 20 and C 30 coefficients describing the Earth’s oblate-ness 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, depending on the source of SLR coefficients and period for which the coefficients are replaced. In SLR and GRACE solutions, all coefficients 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 a degree 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 trend of 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.


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 monthly data at a spatial resolution of 200-500 km enabling us to gain a better understanding of long-term changes in Earth's large ice sheets, land hydrology, and climate through time-variable gravity field and mass change 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 the Greenland 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 stability of 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 differed from 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 (Chen et 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 (Cheng et al. 2013;Cheng and Ries 2023;Loomis et al. 2020).
The GRACE and GRACE-FO missions are equipped with an onboard accelerometers to accurately measure non-gravitational forces such as atmospheric drag and solar radiation pressure and to distinguish them from gravitational perturbations.These accelerometers played a critical role in obtaining precise inter-satellite ranging measurements.Unfortunately, during the last few months of the GRACE mission, the accelerometer on one of 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 accelerometer measurements in the form of the data transplant.This algorithm took into account factors such as flight time, orientation, and satellite-specific events.The GRACE-FO 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 the launch 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 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 C 30 values in GRACE series should be replaced by SLR-based values since August 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 (Bonin et al. 2018;Löcher and Kusche 2021;Sośnica et al. 2015) and the GPS-based orbit of low-orbiting satellites, such as the Swarm mission (Dahle et al. 2020).Combinations of 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 C 61 and S 61 and GRACE solutions at the normal equation level.By combining the data from different mission or analysis centers we can enhance 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 to degree 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 C 20 and C 30 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 equally significant 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 C 20 and C 30 coefficients from different sources are assessed for crucial regions of ice mass depletion and post-glacial rebound, such as Greenland, Antarctica, Fennoscandia, and South Canada with a comparison to external and independent geophysical models.
Firstly, we focus on presenting the effects of replacing spherical harmonics C 20 and C 30 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 replacements of C 20 and C 30 based on own SLR solutions.Due to the fact that estimates of gravity field coefficients are not fully 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 each technique 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 the impact 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.

Data and methodology
The International Association of Geodesy established the International Combination Service for Time-variable Gravity Fields (COST-G) dedicated to the combination of 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 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 Research Centre for Geosciences (GFZ, (Dahle et al. 2018)), and by several other analysis centers, e.g., the Astronomical Institute of the University of Bern (AIUB, (Lasser et al. 2020)), the Centre National d'Etudes Spatiales/Groupe de Recherche de Géodésie Spatiale (GRGS, (Lemoine 2023)), the Institute of Geodesy at Graz University of Technology (IfG, (Mayer-Gürr et al. 2018)), the Leibniz Universitxiät Hannover (LUH, (Koch et al. 2023)), the Chinese Academy of Sciences (IGG, (Wang et al. 2015)), Huazhong University of Science and Technology (HUST, (Zhou et al. 2016)), Tongji University Shanghai (Chen et 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 and GRACE Follow-On.We investigate the impact of these changes on the trend and amplitude of annual signals in Antarctica and Greenland, i.e., the specific regions with strong ice mass loss.We replace the C 20 and C 30 coefficients according to the recommendations provided in Technical Note 14 (TN14; Loomis et al. 2020)

Greenland
West Antarctica East Antarctica procedure for all GRACE solutions.Additionally, we considered other replacements, which are listed in Table 1.

Solution T [Gt yr
We replace the C 20 and C 30 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 the approach described in Sośnica et al. (2015) and relies on the updated version of the Bernese GNSS Software (Dach et al. 2015).The authors introduce a more stable method by 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 noise over the Pacific Ocean compared to the primary results obtained in this study.The solution is expanded to degree and order 10 using the normal equations from three consecutive monthly solutions.

The effects of spherical harmonics replacement on ice mass loss estimates
The Earth's gravity field potential can be represented as a series of spherical harmonic coefficients: where r, φ, are the spherical coordinates in the ref- erence frame, GM E is the product of the gravitational constant and the Earth's mass, R is the semi-major axis of the Earth, C nm and S nm are the Stokes coefficients of spherical harmonics of degree n and order m, and P nm are the fully normalized Legendre polynomials.Thus, the most commonly used representation for the gravity field comparison is the equivalent water height (EWH), which represents the temporal (e.g., monthly) variations within (1) where k n is the Love numbers of degree n , M E is the Earth's mass, and C nm and S nm are temporal variations of the Stokes coefficients.The indicates that the values are after subtracting the signal of the static gravity field part.
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 (2)   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 C 20 and C 30 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 a similar 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 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 the entire 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 period with GRACE accelerometer issues in COST-G (2014COST-G ( -2018) ) and CSR for the accelerometer issues from GRACE and GRACE-FO (2014-2022).Figure 4 compares the C 30 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 than those 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 C 30 , C 50 , C 70 , and C 90 .We can see that better convergence between SLR and the data obtained 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 estimated coefficients, including those, which are correlated with C 30 , 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 pressure and albedo mismodeling effects.Hence, these small 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 C 30 , we can conclude that replacing the entire time series of selected spherical harmonics in models developed to different degrees may result in the loss of information for certain components.Consequently, part of the signal potentially contained in the C 30 coefficient may be found in other coefficients correlated with C 30 .In addition, this approach of simple replacement of two coefficients completely neglects the formal errors (standard deviations) of the spherical harmonic coefficient determined in a given month from the leastsquares adjustment.To fully utilize the potential of the gravity field model and to avoid the commission and omission errors, all spherical harmonics determined up to the least degree and order 10 from SLR should be used because the SLR solutions are still sensitive to time-variable gravity for some coefficients of degree 10 (Sośnica et al. 2015).

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 of equivalent information content among all NEQs.In contrast, the solution level combination, particularly the stochastic combination method, focuses on combining the spherical harmonic coefficients of individual solutions based on their respective standard deviations obtained from the least-squares adjustment (Grombein et al.  2022).In the stochastic combination method, the spherical harmonic coefficients are weighted based on their formal error.Solutions with lower standard deviations, indicating higher reliability, are given greater weight in the combined result.This ensures that more reliable solutions contribute more significantly to the final estimate of the gravity field, while less reliable solutions have a diminished influence.The combination is defined by: where w i nm is the weight of each spherical harmonic based on the inverse of the square of formal error and x i nm is the 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 and order 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 mean values derived only from zonal spherical harmonics, we observe 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 , S 44 , C 43 , 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 defined as monthly datasets, which have different formal errors depending on the time period.Based on Fig. 5, particular attention should be given to the spherical harmonics to degree and order 6.They have been analyzed month by month 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 interferometer measurements, affects the formal errors and significantly reduces them.Another evidence is the issues related to the end stage of the GRACE mission connected with problems 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. (3 The proposed combination takes full advantage of each gravity 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, in the middle the parameters from the SLR S solution, and on 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 the spherical harmonics, the amplitude of the annual signal corresponding to the seasonal changes described by coefficients, and last one the RMS value which represents the fitting 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 improves formal 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 model is based on data from three months.This provides a slight smoothing effect on individual components while reducing noise originating from SLR observations and enabling 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.

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 C 20 /C 30 coefficients.Specifically, we compare it with the GRACE/GRACE-FO solutions derived from the COST-G and CSR RL06 approaches.According to the recommendations of TN14, we also replace the spherical harmonic coefficients C 20 and C 30 in the COST-G data with the corresponding and from SLR S. It is worth to mentioning that the analyzed combination is performed using data from April 2002 onward when GRACE solutions became available.The results regarding the basic parameters of the data series, such as trends and the amplitude of the annual 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 which of the combined solutions generates the lowest noise over 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 over oceans 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 (Caron et al. 2018).Consequently, after reducing combined SLR-GRACE models by GIA, the trend in EWH should be close to zero.
For each of the analyzed regions with significant ice mass loss, we observe that when the spherical harmonic coefficients 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.3Gt/year after the replacement with C 20 and C 30 values 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 from TN14 and SLR S, respectively.For the combined solution (COMB), the value is − 112.5 Gt/year.Thus, the combination 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 in the 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 amplitude is approximately twice as large as in the case of replacing the C 20 and C 30 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 October 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 the analyzed 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, and 52.2 Gt/year for COMB 2, 4, and 6, respectively.When DOW 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.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 Western part of Antarctica, the trend is − 111.8 Gt/year, and the amplitude is 52.1 Gt.Therefore, the COMB DOW combination yields results that fall between COMB 2 and COMB 4 in terms of their outcomes.

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 to large formal errors in GRACE solutions.For C 50 obtained from SLR, there is a significant underestimation of the amplitude 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 shows the 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 estimates models complete up to degree and order 5 plus C 61 and S 61 .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 information and 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 C 44 and S 42 coefficients, which are strongly influenced by seasonal variations and confirms that 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 in terms of the low-degree coefficients.Figure 9 shows the results 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.
Figure 10 depicts the sectorial and tesseral coefficients C 33 , S 22 , C 22 , and S 32 , considering combinations with different formal error scaling factors, from 2012 onward.A noticeable tendency can be observed in the means of a slight 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 the data.The appropriate scaling of formal errors becomes a crucial aspect of proper representing seasonal signals and trends included in the spherical harmonic coefficients.
Figure 11 depicts the differences between the trend from the GRACE/GRACE-FO COST-G with replaced C 20 and C 30 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 in polar regions.The non-zonal spherical harmonics, which are relatively unrelated to large-scale phenomena and have regional characteristics, affect the remaining areas, such as oceans and some land areas.Therefore, different handling of low-degree gravity field coefficients may have an impact on the estimates of secular trends in sea level rise or trends in land hydrology which are important for the draughts and flood predictions.Reducing the weight of the SLR solutions reduces the differences between GRACE-only and the GRACE-SLR combination, but the remaining differences are still visible in long-term trends with a smaller magnitudes and slightly alter them compared to the solution adjusted according to TN14 C 20 and C 30 .

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 observed in 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 the end 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 RMS value.Table 3 shows that the combinations with a reduced weight of SLR by 4 times demonstrate comparable median noise levels over the oceans as the solution with C 20 and C 30 replaced according to TN14.For the entire period, this median noise level amounts to 2.6 cm.However, the noise in COMB 4 is slightly smaller in the first 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 the weight of the SLR solution or COMB DOW with degreespecific weighting achieve similar or better results to the solution based on TN14 data replacement.Therefore, downweighting by a factor of four or degree-specific weighting 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) and calculated the trends for the Scandinavian region and South 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, the resulting 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 such as Scandinavia, where the expected seasonal hydrological signals are small.A low value suggests that the data is less noisy and does not contain significant disturbances from the solution.The trend value is lower for both the Scandinavia region and South Canada in COMB 2 than in 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 can conclude 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 over oceans than COMB 2. The residual trends and STD values are larger in South Canada than Scandinavia because of the impact of ice mass depletion and seasonal land hydrology variations in Canada.Moreover, the GIA effect is 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 also decide 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 the remaining data signal.In this case, the results do not indicate that the combined model achieves the best outcomes but rather strengthen the suggestion that the contribution of SLR data in the final solution should be reduced, as COMB solutions exhibit much higher STD compared to others, such as COMB 2 and COMB 4. The secular trends 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.

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 sur 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 based on the existing Simplified Land Surface Scheme (SLS) and Hydrological Discharge Model (HDM).In LSDM the amplitude 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 is observed at the start of the analyzed period, spanning 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 from 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 by differences in determined low-degree spherical harmonics.We do not observe significant distinctions between our proposed combination and the replacement of C 20 and C 30 with TN14.

Discussion and conclusions
In this article, we analyzed the impact of replacing spherical harmonics C 20 and C 30 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 C 20 and C 30 from SLR S in the same periods as the TN14 data, we obtained similar trends but with larger amplitudes of seasonal signals in Greenland and Antarctica.This was particularly evident in the East Antarctic region.The influence of different C 30 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-specific solution, therefore, the gravity signal can be distributed over different coefficients describing similar mass distribution.Therefore, we performed a combination based on the formal errors of the spherical harmonics.For the combination at the solution level, we used time-variable gravity field data from SLR S and GRACE/GRACE-FO COST-G data.
For SLR solutions, we used the expansion up to degree and 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 the weights 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 to the sensitivity of SLR solutions to time-variable gravity field changes for which the SLR-only solutions are still stable.The calculated average contribution showed that even 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 showed that during the period 2015-2017, SLR data accounted for approximately 50% of the contribution used for combinations up to d/o 6.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/year in the Western Antarctic region, while the trend value replaced with C 20 and C 30 values from TN14 was − 111.6 Gt/year.For the combined solution COMB, the value was in the middle equaling to − 112.5 Gt/year.To assess the quality of the solution, we calculated RMS values over parts of the Pacific Ocean, which indicate noise and thus should be minimized.The computed values showed that the combination slightly increased the noise, with a value of 2.9 cm for the combination and 2.6 cm for the replacement of C 20 and C 30 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 and reached the level of replacing data based in TN14 or even slightly better.The combinations also introduced changes in areas related to ice melting.In the case of the eastern part 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 is correct and no secular hydrological changes occur in these regions, 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 C 30 with TN14, this value 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 other low-degree gravity field coefficients from SLR.Therefore, the combination of SLR and GRACE data provides more consistent results than replacing two coefficients in GRACE series by SLR-derived values for arbitrarily selected periods.

Fig. 1
Fig. 1 Mass differences in Greenland and Antarctica due to various replacements of C 20 and C 30 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 each region then, COST-G with replaced C 20 or C 30 or both 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 C 20 and C 30 for the same months with our SLR S, the trend and amplitude are − 178.8 Gt/year and 53.1 Gt, respectively, which implies a change in the amplitude of about 20% w.r.t.TN14.More pronounced changes are observed in the Antarctic region.For East Antarctica, the trend changes from 41.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 with the original GRACE data.Replacing C 30 only since March 2012 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

Fig. 2
Fig. 2 Mass differences in Greenland and Antarctica due to various replacements of C 20 and C 30 values in the GRACE series.The COST-G GRACE solution was subtracted from modified solution with C 20 and C 30 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

Fig. 3
Fig. 3 C 20 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. 5 Fig. 7
Fig. 5 Median percentage contribution of SLR (left-hand side) and GRACE/GRACE-FO COST-G (right-hand side) for each spherical harmonic coefficient in the combined SLR-GRACE solution from 2002.04 to 2021.10

Fig. 10
Fig. 10 The difference between the GRACE/GRACE-FO COST-G and the combination with different weights of SLR solutions face discharge model (LSDM, Dill 2008) and the global land water storage dataset assimilating GRACE and GRACE-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.0 hydrological model, supported by GRACE/GRACE-FO

Fig. 13
Fig. 13 Time series of EWH from COMB 2, COST-G with replaced C 20 and C 30 from TN14, land and surface discharge model (LSDM) and GLWS 2.0 for Amazon and Brahmaputra river basins

30 component. TN14 provides a standard Table 1 Impact of the replacement of C 20 and/or C 30 on trend (T) and amplitude (A) in ice mass change for the regions of Greenland
, West Antarctica, and East Antarctica based on COST-G solutionsThe C 30 replacement starts in March 2012 or is valid for the whole period (full)

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 effect

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

Table 4
Trend and long-term standard deviations (STD) of EWH in the Scandinavia and South Canada regions calculated with removing the GIA effectThe STD* was calculated with the removed hydrologically LSDM model, trend and the annual and semiannual signal