A comparative study on the parametrization of a time‑variable geopotential model from GRACE monthly solutions

The gravity field of the Earth is time-dependent due to several types of mass variations which take place on different time scales. Usually, the time-variability of the gravitational potential of the Earth is expressed by the monthly determination of a static geopotential model based on data from gravity field missions. In this paper, the variability of the potential is parameterized by a functional approach which contains a polynomial trend and periodic contributions. The respective parameters are estimated based on the monthly solutions derived from the GRACE and GRACE-FO gravity field mission up to a maximum degree of expansion n max = 96 . As a preliminary data analysis, a Fourier analysis is performed on selected potential coefficients from the available monthly solutions of the GFZ. The indicated frequency components are then used to formulate a time-dependent analytical approach to describe each Stokes coefficient’s temporal behaviour. Different approaches are presented that include both polynomial and periodic components. The respective parameters for modelling the temporal variability of the coefficients are estimated in a Gauss-Markov model and tested for significance by statistical methods. Extensive comparative numerical studies are carried out between the newly generated model variants and the existing monthly GRACE, GRACE-FO and the existing time dependent EIGEN-6S4 solutions. The numerical comparisons make it clear that estimated models based on all available monthly solutions describe the essential periods very well, but such monthly events that deviate strongly from the mean behaviour of the signal show less precision in the space domain. Models that are estimated based on fourteen consecutive monthly solutions, covering one selected year, represent the amplitudes much more precise. The statements made apply to four initial data used, which are filtered to varying degrees. In particular, DDK2, DDK5 and DDK8, as well as unfiltered coefficients were used. For all the model approaches used, it can be seen that the potential coefficients contain up to about n ≈ 40 in case of DDK5 periodically signals with annual, semi-annual or quarterly, as well as Luna nodal periods and do not vary significantly beyond that degree. Only an offset can be estimated significantly for all Stokes coefficients.


Introduction
The determination of the Earth's external gravity field with its temporal and spatial variations is one of the main tasks of physical geodesy.It contributes to the exploration of the Earth system (Flury et al. 2017).Based on the deep knowledge about the behaviour of the global gravity field, many geodynamic processes can be measured and understood.By recording its temporal variations, it is possible to infer mass displacements, fluctuations in sea level, deep ocean currents, run-off and groundwater storage on landmasses, interactions between ice sheets or glaciers and the oceans, the rate of ice melt, changes in groundwater storage and the extent of droughts.Reliable statements about non-linear processes are also possible.Furthermore, the geoid as a selected equipotential surface of the gravity field presents a basis for geodetic reference systems to determine vertical movements and unify height systems (Sánchez et al. 2021;Zhang et al. 2022;Grombein et al. 2016Grombein et al. , 2017)).The gravity field missions GRACE (Gravity Recovery And Climate Experiment) and its follow-up mission GRACE-FO provide observations with global homogeneous precision (Tapley et al. 2019;Ince et al. 2019), which can be used to derive functionals of the gravity field of the Earth including its behaviour with respect to time.This is realised in terms of time-dependent spherical harmonic coefficients (Stokes coefficients).Each coefficient is represented by a time-varying function whose parameters are estimated from the monthly GRACE solutions.
The EIGEN (European Improved Gravity model of the Earth by New techniques) models by Förste et al. (2016) introduce a series of time-variable geopotential models from the data of satellite gravity field missions.For the satellite-only EIGEN-6S4 model, longterm observations from the LAGEOS (LAser GEOdynamics Satellite) as well as data from the GRACE and GOCE missions are used to determine time-dependent parameters for the Stokes coefficients (potential coefficients) up to degree and order n = 80.
The presented work deals with the determination of a time-dependent geopotential model from monthly solutions of the GRACE and GRACE-FO gravity field mission.Each Stokes coefficient is described by a time-dependent functional approach.To specify the functional modelling, which is in the focus of this investigation, the GRACE solutions are first subjected to a Fourier analysis.This analysis qualitatively confirms the expected drift and periods in the Stokes coefficients.Within the framework of a Gauss-Markov model, different sets of time-dependent parameters for the Stokes coefficients are determined.This contains periodic and polynomial parts.The estimated parameters are tested for their significance.These parameters are then used to calculate the gravitational potential and geoid undulations at specific epochs.To assess the quality of these newly evaluated time-dependent models and the adopted functional approach, they are compared in the space and frequency domain with existing monthly GRACE and GRACE-FO solutions (ICGEM 2023) and with the time-dependent EIGEN-6S4(v2) model.
In further studies beyond the scope of this article, the generated time-varying models can be used to bridge the data gap between GRACE and GRACE-FO.AI (artificial intelligence) methods could also be explored to perform this task.
The structure of the paper is as follows: The gravitational potential as a time dependent function is described in Sect. 2. The approaches developed for this purpose are described here in detail.The least-squares estimation in the Gauss-Markov model and applied statistical tests regarding the significance of the respective parameters are compiled in Appendix 2. A short timeline of the satellite gravity field missions GRACE and GRACE-FO is drawn in Sect.3. In addition, the products derived from the observations are specified.The geopotential models of the GFZ (Deutsches Geoforschungszentrum Potsdam) which are used as input data in order to derive a time-dependent geopotential model are explained.Section 4 presents the evaluation of the estimated time-dependent parameters up to degree and order n max = 96 .The parameters describing the time-variability are tested with respect to their significance.Several estimated time-dependent sets of Stokes coefficients are discussed in Sect. 5. Finally, Sect.6 contains a summary and provides an outlook on future (Hartmann andWenzel 1994, 1995;Wenzel 1998) and the relation to a dedicated tidal system (Ekman 1988(Ekman , 1989)).The atmospheric effects can be reduced approximately by the use of weather models (ECMWF 2022; Kumar et al. 2009).For the global observation of these variations, satellite missions on near-Earth and polar orbits are suitable.
It is obvious from Table 1 that the Stokes coefficients will variate with different frequencies caused by seasonal reasons.Even longer periods can be resolved if the responsible forces can by identified, like the nutation of the rotational axes of the Earth.Otherwise additional parameters can be introduced (quadratic and cubic terms) to describe (mathematically) those effects which is common practise in geodesy.This topic will be discussed in more detail in Sects.2.2 and 4. Time dependent effects on Z(r, , t) , such as those listed in Table 1, can be taken into account by using appropriate models.This is accompanied by the establishment of a geodetic datum (ITRF, Altamimi et al. 2011Altamimi et al. , 2016)).Their periods are not investigated further in this paper.
The static gravitational field is commonly represented in spherical harmonics (Heiskanen and Moritz 1967): The geocentric gravitational constant is the product of Newton's gravitational constant G and the mass M of the Earth.The radius of the reference sphere is denoted by R. The fully normalized Stokes coefficients are denoted by C nm and S nm .The Laplace's spherical harmonics Y nm , which are also fully normalised, create an orthonormal base system (Heis- kanen and Moritz 1967): The P nm are the fully normalized Legendre functions of the first kind of degree n and order m with 0 ≤ m ≤ n .Based on the orthogonality relations of Laplace's spherical harmonics (2), the power spectrum or absolute degree variances can be defined: Analogously, absolute error degree variances 2 n c n can be determined: The variances of the Stokes coefficients C nm and S nm are given by C nm and S nm , respectively.The signal-to-noise ratio of the gravitational potential can be represented in the frequency domain by comparing (3) and ( 4).

Time-dependent functional models
The representation ( 1) is used to model the time-dependent gravitational potential by introducing fully normalized Stokes coefficients C nm (t) and S nm (t) which are timedependent.It was decided to model the time-dependent coefficients f(t), the effects on the gravitational potential, by a trend and periodically approach: A trend is represented in terms of a polynomial up to degree q max .The periodic behav- iour of the time dependent potential coefficients is modeled by periods P i .In Sect. 5 it is brought forward the argument that q max = 3 is recommended to model the trend with a cubic polynomial.A spectral analyses by FFT (Fast Fourier Transformation) supports to select the periods of the oscillations to P 1 = 1 y, P 2 = 1 ∕2 y, P i = 1 ∕i y.Due to the sampling rate of one month, p is limited to p ≤ 6 (Jekeli 2017).
In addition, considering the analyzes in Sect. 4 revealed that modeling the lunar nodal cycle of P L = 18.6 y (Roy 2004) appears to make sense.This period appears in the nutation of the Earth's rotational axes and also in the solid Earth tides (Benjamin et al. 2006;Botai and Combrinck 2012).
Five representative analytical approaches are taken into account to formulate and find the optimal functional model.The functional relationships are based on the one used by Förste et al. (2016) and Grombein et al. (2021) and on those supported by the results of the Fourier analysis of the monthly GRACE solutions.A summary of the different approaches, which are subsets of the general approach in Eq. ( 5), can be found in Table 2.
In the first functional model f 1 , for the time series of each coefficient, the trend is modeled by an offset c 0 , a secular term c 1 (t − t 0 ) , and the amplitudes a 1 and b 1 of annual periods according to the functional relation: (5) The f 2 corresponds to the EIGEN-like approach (Förste et al. 2016) Coefficients for P L with the yearly periods of the oscillations P 1 = 1 y have to be estimated.
In the second functional model f 2 , the model f 1 is extended by the semi-annual peri- ods P 2 = 1 ∕2 y with the amplitudes a 2 and b 2 : to be estimated.The function f 2 corresponds to the modelling of EIGEN-6S4(v2) (Förste et al. 2016).
In the third functional model f 3 , for the time series of each coefficient, an offset c 0 , a linear trend c 1 , and the amplitudes of annual a 1 and b 1 , semi-annual a 2 and b 2 and quarterly a 3 and b 3 periods are calculated according to the functional relation: with the quarterly periods of the oscillations P 3 = 1 ∕4 y to be estimated.This takes into account temporal variations in the gravitational field that are due to periodic or long-term mass variations such as seasonal changes in the continental water balance and long-term sea level changes (Cazenave et al. 2014) caused by variations in the mass balance between ice and seawater due to global warming.The functional relations f 1 and f 2 are subsets of f 3 .They are considered in order to evaluate the individual periods of the oscillations.
The fourth model f 4 , in addition to f 3 , also takes into account quadratic and cubic trend terms for a better local fit to the time series, i.e. it contains a third-degree polynomial: The fifth model f 5 , in addition to f 3 , includes the long period of P L = 18.6 y , the lunar nodal cycle (nutation): Since it is expected that the long-period lunar nodal effect of 18.6 y is completely modeled by a L and b L , neither a quadratic nor cubic trend term is estimated in the functional model f 5 in order to avoid over-parameterization.
Based on the input values listed in Table 11 (see Appendix 2) the parameters are estimated for each of the k Stokes coefficient according to the functional model f j .For each coefficient, the estimated parameters listed in Table 12 (see Appendix 2) result.The representation of a oscillation with two periods of amplitudes a and b can be transformed into the phase-modulation with the phase angle and amplitude A:

Determination and significance test of the time-varying parameters
In the framework of the Gauss-Markov approach the model parameters which are introduced in Sect.2.2 are now determined.The estimated parameters, as well as the assumed functional and stochastic models are then tested for significance (see Appendix 2).

Monitoring the time-variable gravity field of the Earth
Through GRACE (GFZ 2021c), it was possible to determine the Earth's gravitational field and its changes over time with high precision (Dahle et al. 2018a, b).The mission ended in October 2017 after about 15 years.Since May 2018, the two satellites of the GRACE-FO mission have been in orbit to continue the objectives of GRACE.The GRACE mission utilised electromagnetic radiation in the microwave range (K-band).GRACE-FO also uses such a sensor to collect the primary observation data.In addition on GRACE-FO a laser interferometer is installed, to carry out a technical experiment for future missions.The aim of the laser interferometer is to achieve higher precision of the range measurements (GFZ 2021b).The Intergovernmental Panel on Climate Change underlines the important contribution of GRACE to the observation of the temporal mass variation in the Earth system (Portner et al. 2022).

GRACE and GRACE-FO products
The products of the GRACE and GRACE-FO missions can be obtained from three analysis centres, namely the Jet Propulsion Laboratory of the California Institute of Technology (JPL, Caltech), the Center for Space Research of the University of Texas at Austin (CSR, UT) and the German Research Centre for Geosciences in Potsdam (GFZ 2021a).
Level-2 products comprise sets of fully normalised potential coefficients C nm (t k ) and S nm (t k ) for each month k to describe the Earth's external gravitational potential according to (1).The Level-2 data products are generated by the three analysis centres through the use of different processing techniques.The data can be downloaded from the ICGEM website with (uncorrelated) variance information (ICGEM 2023).They are available with different levels of smoothing as listed in Table 3.The applied filters are explained in Kusche et al. (2009).Consult also Wahr et al. (1998) and Swenson and Wahr (2006) for the topic of filtering the Stokes coefficients itself, instead of filtering in the space domain.

Utilized database for the calculations
The data basis for the calculations within the scope of this work are fully normalised Stokes coefficients up to degree and order n max = 96 , determined by GFZ from GRACE monthly observations.The latest Release 06 (R06) is used where different degree-dependent correlated noise filter DDK (Kusche et al. 2009;Guo et al. 2017) are applied for smoothing (Chen et al. 2021).If nothing is specified, the investigations in this article refer to DDK5 filtered coefficients.Studies were also carried out with the DDK2, DDK8 and unfiltered coefficients, which are also available on the ICGEM web-page.
The variations within the monthly solutions of the Earth's gravity field reflects the changes in surface and deep ocean currents, run-off and groundwater storage on landmasses, interactions between ice sheets or glaciers and the oceans, and mass displacement within the Earth (Angermann et al. 2022).This so called level-2 products can be downloaded from the International Centre for Global Earth Models (ICGEM) website with (uncorrelated) accuracy informations (ICGEM 2023).
The solutions at the beginning of the GRACE mission have some gaps of several days, which is why a total of 161 out of the 163 available monthly solutions for the period from 2002 to 2017 are included in this research study.Appendix 1 contains an overview of the available monthly solutions of the GRACE and GRACE-FO missions in Table 10.
In addition to the monthly solutions, the time-variable geopotential model EIGEN-6S4(v2) (Förste et al. 2016), published by the GFZ in cooperation with the Groupe de Recherche de Géodésie Spatiale (GRGS), is used for the comparison to assess the quality of the newly developed time-dependent GRACE model in the space and frequency domain.These EIGEN models are based on long-term observations from the gravity field missions GOCE, GRACE, CHAMP and/or LAGEOS.The parametrization of the time-dependent coefficients up to degree and order n max = 80 are described in Rudenko et al. (2014).An offset, a trend and the amplitudes of an annual oscillation and a semi-annual oscillation are estimated for each year from 2002 to 2014 (Abd-Elmotaal 2020;Rudenko et al. 2014).This corresponds to the determination of p max = 2 and q max = 1 in the general model approach given in Eq. ( 5).In the present work this corresponds to the model approach f 2 .The result- ing time-dependent behaviour of the coefficients is continuous; at instances of strong earthquakes, for example in Chile 2010 (Tong et al. 2010) or in Japan 2011 (Sato et al. 2011), an offset can be detected.

Calculations and analysis
Within the framework of a least squares estimation, time-dependent parameters are estimated in the Gauss-Markov model in order to represent the fully normalised Stokes coefficients up to degree and order n max = 96.
Standard deviations for these parameters are computed via variance propagation.Additionally, the parameters are tested for significance.To assess the quality of the new timevariable geopotential model, various comparisons are carried out in the space and frequency domains.

Preliminary investigation for the frequencies to be expected
The discrete Fourier transformation (DFT) (Jekeli 2017) is used to analyse the spectral properties of selected coefficient time series.The aim is to show which frequencies are to be expected in the signal and whether the model approaches explained in Sect.4.2 are justified.This is introduced as prior information for the formulation of the functional model ( 5) for the estimation calculation.To do this, the offset and trend in the data are removed via a least-squares fit.The input signal is not smoothed, but interpolated to a sampling rate dt = 1 d as a regular measure of time, also in order to close the gaps towards the end of the GRACE time series x n .The resulting time series is N = 5413 days long.A Hanning win- dow is applied to the time series x n in order to reduce leakage effect.A Hanning window is used because the Fourier transform of a Hann window results in a Hann window again.Thus, the effect of the rectangular window in the time domain, whose FFT is the sinc function, is weakened (Testa et al. 2004).From the discrete FFT, the amplitudes of the corresponding frequencies of the oscillations occurring in the time series are determined and plotted (cf.Fig. 2) as amplitude spectra of the signal.

Determination of the time-variable geopotential model
According to Appendix 2 the time-variable coefficients are determined.The quality of the estimation results is checked using hypothesis testing, which is also provided in Appendix 2.
All available 161 monthly solutions of the GRACE-mission filtered with DDK5 (Kusche 2007), complete to degree and order n max = 96 and the constants listed in Table 4 are entered in the processing of the time-variable gravitational field.
The total number of coefficients is listed in Table 4. From C nm and S nm , the observation vectors C nm and S nm are set up for each coefficient and gathered in matrices.In order to avoid numerical instability, the value 0 = 10 −12 is factored from the standard deviations for the stochastic model.The a priori variance factor is thus 2 0 = 10 −24 .Since there is no information on correlations between the observations available on ICGEM (2023), and  5) and ( 19).For this purpose, the reference date of the first observation is selected as the reference epoch t 0 .The Julian date is introduced as time parameter.The calculation in days leads to unfavourably scaled normal equation matrices, the inverse of the condition number (N) is numerically close to zero with 1∕ ≈ 10 −24 .For this reason, calculations are car- ried out in years for numerical stabilization.
By introducing the Stokes coefficients C nm and S nm as uncorrelated (pseudo-) observations C nm and S nm , because no correlations are available, they can be analyzed independently of each other.The five considered functional representations of the time variable gravitational potential are presented in Sect.2.2.
To verify the stochastic and functional model, a test according to Eq. ( 27) with a significance level of = 5% is performed.It should be remembered that the statement whether H o or H A is correct obviously depends on the choice of .It would also be conceivable to determine the value at which H 0 is just accepted.If the null hypothesis is accepted, the variances of the estimated parameters are calculated using the a priori variance factor From these variances, the statistical test quantities T prio or T post for the parameter tests are calculated according to Eq. ( 29) and the limits for defining the rejection or acceptance of the null hypothesis are determined depending on the result of the global test.The procedure for the two-dimensional parameter tests according to Eq. ( 30) is analogous.
For verification by comparison with existing monthly solutions, the respective month is removed from the observation vectors.This reduces the number of observations n to n V = 160 .Then, the result provides a measure of the precision of the applied functional model.

Verification in frequency and space domain
The starting point for the comparisons of the time-dependent geopotential models according to the functional relationship f 1 (6), f 2 (7), f 3 (8), f 4 (9) and f 5 (10) are the estimated parameters xj and their variance-covariance matrix C xx j .The derived time-dependent GRACE models are compared with an existing monthly solution.For this purpose the time-dependent model is created in terms of fully normalized Stokes coefficients on the basis of n V = 160 monthly solutions at a certain interpolation time t.This returns fully normalized potential coefficients C nm 1 and S nm 1 using four estimated parameters according to Eq. ( 6), six estimated parameters according to Eq. ( 7), eight estimated parameters according to Eq. ( 8), and ten estimated parameters according to Eq. ( 9).The variances of the new coefficients result from variance propagation.According to (Niemeier 2008, p. 56) the row vector A i corresponds to a row of the design matrix A .The variance 2 of a coefficient K nm can be written as: In the next step, the differences ΔC nm and ΔS nm between the calculated coefficients C nm 1 and S nm 1 and the actual coefficients C nm 2 and S nm 2 of the GRACE model of the considered months are created: For the comparison in the space domain, the potential is synthesized on a global grid for the difference model.The synthesis (see Eq. ( 1)) of the difference ΔV( , ) is carried out according to the parameters listed in Table 13 which is presented in Appendix 3. The effect on geoid undulation is calculated from Brun's theorem (Torge 2003): with normal gravity (Moritz 1980).
In addition to a comparison in the space domain, the differences in the Stokes coefficients are also considered in the frequency domain.For this, the absolute degree variances are calculated according to Eq. ( 3).The error degree variances for the difference model result from the linear variance propagation law: The comparison of the developed time-dependent geopotential model from GRACE data with the EIGEN-6S4(v2) is performed also in the space and frequency domain.The timedependent parameters of the EIGEN-6S4(v2) with its standard deviations are used.Since these parameters are determined in one-year intervals, the corresponding interval from a total of 14 epochs is selected to calculate a new model.The differences ΔC nm and ΔS nm between the calculated GRACE coefficients C nm 1 and S nm 1 and the calculated EIGEN- 6S4(v2) coefficients C nm E and S nm E are defined as: Several statistical parameters such as the minimum, maximum and mean value, as well as the standard deviation (std) are evaluated for the presentation and comparison of the sum of squared residuals and the effect on the geoid undulation N in Sect. 5.
All evaluations and analyses are carried out with software realisations in Matlab ® .

Results
At first, the individual monthly solutions are compared in the frequency domain and then the results of the Fourier analysis are presented.The results of the least-squares estimation in the Gauss-Markov model are presented together with the results of the global test and the parameter tests.To assess the external precision, a comparison is made in the space and frequency domain with a GRACE monthly solution (ICGEM 2023), GRACE-FO monthly solutions and finally the new time-variable GRACE model is compared with the EIGEN-6S4(v2) model (Förste et al. 2016). (13) (16)

Investigations of the time series of given monthly solutions in the frequency domain
For the preliminary examination of the monthly solutions, the absolute degree variances of the individual models are considered.They coincide graphically as can be seen in Fig. 1a.
Residual models are therefore created for the individual months in respect to the arithmetic mean values of the time series from 2002 to 2017.Their respective degree variances show some differences as expected (Fig. 1b).To show the noise content of the fully normalized Stokes coefficients, the error degree variances are presented in Fig. 1c.From the Fig. 1b, c it can be seen that the majority of the degree variances and error degree variances of the monthly solutions lie in one coherent frequency band.Models outside of this band are ) reports for January 2015 an unfavourable ground track coverage due to short-period repeat orbit patterns.This caused problems in computing the data products and also influenced the solution for February 2015.The impact can be observed in Fig. 1.The same conditions occurred in September and October 2004 (GFZ 2021a).
In addition to the consideration of degree variances and error degree variances, a Fourier analysis was carried out to examine the signal contained in the time series in order to select the appropriate functional model (cf.Eqs. ( 6)-( 10)).The GRACE zonal coefficient C 3,0 is presented as example.Figure 2a shows the time series for the coefficient C 3,0 in its original monthly basis, the time-series that was interpolated to daily resolution as well as the daily time series after applying a Hanning window.The width of the Hanning window corresponds to the length of the time series.This makes the effect of the window visible: using a Hanning window, the observations at the beginning and end of the time series are weighted less than the observations in the middle.Since the window is normalized, the energy content of the signal is preserved.The values fluctuate around zero because the mean and trend have been removed via a least-squares fit.The associated amplitude spectrum is shown in Fig. 2b.From this, the effect of the window in the frequency domain becomes obvious: the amplitudes of the secondary maxima are slightly higher, the broadening of the main maxima at approx.0.003 d −1 and approx.0.006 d −1 however, cannot be avoided and is also slightly wider due to the windowing.The first maximum is at the frequency f = 0.00018 d −1 .This signal corresponds to the frequency of the applied Hanning window and is therefore not included in the frequency spectrum of the signal without the applied window.With a frequency of around f = 0.00037 d −1 , a frequency component that corresponds to 7.5 years appears, which at least gives an indication of long-period components in the signal.The maximum at f = 0.0028 d −1 corresponds to a period of 361 days and thus indicates an annual oscillation.The maximum at f = 0.0055 d −1 corresponds to an oscillation with a period of about half a year.At approx.f = 0.0115 d −1 a weak signal with a period of about three months can be seen.These findings are used to set up the functional models in the estimation.
The FFT was also carried out for C 2,0 and showed a comparable result but with a higher noise level.This is to be expected since the dynamic form factor is difficult to determine from GRACE observations.C 2,0 can be determined with higher accuracy from Satellite Laser Ranging observations (Cheng and Ries 2017;Dahle et al. 2019b).Furthermore, the Fourier analysis is performed for other coefficients, as an example for the time series of the zonal coefficient C 90,0 .The significant signal components for coefficients of a higher order are much more difficult to identify because they are associated with higher levels of noise.Nevertheless, for C 90,0 there is also a maximum at f = 0.0052 d −1 (period P ≈ 192 days ) and f = 0.0105 d −1 (period P ≈ 95 days ).Both the error degree variances in Fig. 1c and the Fourier analysis show that the higher degree coefficients are associated with more noise.The maxima in the amplitude spectrum are less clearly identifiable.

Presentation and evaluation of the estimation results
The time dependent functions f j are fitted to the GRACE time series of each Stokes coef- ficient.The time series of the coefficient C 2,0 with the estimated observations is shown in a joint figure and the associated residuals are presented (Figs. 3,4,5,6,and 7).Temporal changes in the C 2,0 coefficient are directly related to variations in the flattening of the Earth.On the one hand, the time series shows annual oscillations due to seasonal variations in the ice cover at the poles and the amount of water stored in the large river systems.On the other hand it shows long-term changes in the form of a negative trend for the years from 2002 to 2017 due to the accelerated melting of ice masses caused by climate change, which superimposes the signal of the post Ice Age land uplift.Therefore, there has been a decrease in flattening since 1998 (Förste 2013).
The representations in Figs. 3, 4, 5, 6, and 7 show that the observations in the middle of the time series have lower residuals and therefore fit better to the functional model.The regular temporal variations can be modelled well in this period, since GRACE models are available for each month.Towards the end of the time series, gaps in the observation data sets occur more frequently, which is why there are several months for which no solution exists (cf.Table 10).The residuals are therefore also larger than in previous years.This is evident for all applied functional relationships f j , 1 ≤ j ≤ 5 presented.A The sum of squared residuals for the observations of the C nm and S nm for the different analytical approaches is determined as a quality indicator for the estimation in the framework of the applied global test.For better comparability, statistical parameters such as the arithmetic mean, the standard deviation as a representation of scatter and the maximum and minimum weighted sum of squared residuals are determined using the estimation results of all coefficients (Table 5).From this it can be deduced which approach best suits the data.Fitting the time series of each Stokes coefficient with the functional relationships f 4 and f 5 give in average for all coefficients the smallest weighted sum of squared residuals.The scatter and thus the interval between maximum and minimum are also smaller than with f 2 and f 3 .There are minor differences between f 2 and f 3 .This indicates that the quarterly oscillation, included in f 3 , may not be significant for all coefficients.
The respective minimum parameters for f j are highlighted in yellow.The spectral com- ponents of the sum of squared residuals for the tested functional approaches f j are evalu- ated.Table 5 summarize their statistical parameters.In a functional context, f 1 takes into account an offset, a linear trend and an oscillation with a period of one year (Table 2) and on average leads to the highest of the weighted sum of squared residuals.The std values for f 2 and f 3 are nearly the same, as highlighted in blue.This shows that considering only an annual oscillation leads to less adapted functionals compared to those including the semiannual or quarterly oscillation.The additional consideration of the periods of six and three months leads to a smaller sum of squared residuals (Table 5, f 3 ), since random signal components can also be reduced with it.
In addition to the weighted sum of squared residuals, the results of the global test can also be used to assess the fit to the analytical approaches (Figs. 9, 10).The results for the The values highlighted in bold are similar for f 2 and f 3 .The statistical values for f 4 and f 5 are highlighted in italics and show the smallest, very similar numbers different functional relationships are very similar: The test ( H 0 -hypothesis) is rejected for the coefficients up to degree around n ≈ 40 and accepted for the coefficients of a higher degree and order.A rejection means that the a priori and the a posteriori variance factor from the estimation differ significantly ( = 5% ).The individual test statistic T G for each Stokes coefficient and the five different analytical approaches ( f 1 to f 5 ) used to determine the time-dependent geopotential model are shown in Fig. 8.The test statistic is computed according to Eq. ( 27), which involves the ratio of the a priori and a posteriori variance  factor weighted by the redundancy.Consequently, the test statistic T G is unit-less.For the zonal and near-zonal tesseral coefficients the value of the test statistic is highly significant (dark green).For coefficients with degrees larger than n ≈ 40 , the value of the test statistic falls below the critical value 2 r which is in the order of ≈ 180 (cf.Table 6) and lead to acceptance of the H 0 hypothesis.The coefficients with multiples of order m = 15 are more significant due to the resonant orbit perturbations, which can be nicely observed.The result of the global test is shown in Fig. 9.
The reasons for this are either inconsistencies in the stochastic model or in the functional model-compared to the observation data.Since this pattern is roughly the same for all approaches, it can be concluded that the functional model is unlikely to be the reason for the rejection of the test.However, it is possible that there are discrepancies in the stochastic model for the observations of the coefficients.The higher degree coefficients are subject to higher variances than the lower degree coefficients (Fig. 1c).It is therefore possible that the uncertainties for the observations obtained from a variance propagation from the uncertainties of the satellite data are too optimistic and do not fit to the actual prevailing conditions, which can lead to a rejection of the global test.Correlations between the data that are not taken into account can also affect the test result.
These statements relating to DDK5 filtered solutions can be supplemented with the following observation: As expected, the significance limit also depends on the filter strength  (Table 3), as can be clearly seen in Fig. 10.While in DDK2 the global tests for the coefficients below n = 30 are significantly, in DDK8 the significance limit is n = 60 .For the unfiltered data, 0 < m < n − 60 with n > 60 is no longer significant.for f 4 based on the DDK8 data are presented in Fig. 13.The results of the other approaches are very similar and are therefore not shown.If the parameter test is rejected (black) for the parameters related to f 4 , the parameter is significantly different from zero.If, on the other hand, the null hypothesis of the test is accepted (yellow), no significant difference from zero can be determined for the value of the parameter.With = 5% the null hypothesis is rejected in 5% of the cases, although it is true.This shows that the offset is determined significantly for almost all coefficients.The linear term, the quadratic term and the cubic term of the trend, as well as the amplitudes of the annual and semi-annual oscillations differ significantly from zero for most coefficients up to about degree n ≈ 40 (DDK5).From this, it can be deduced that the Stokes coefficients of higher degree and order are subject to less strong temporal variations or have shorter periods, or are less precise, which is why the time-dependent parameters are determined to be insignificant.The amplitudes of the quarterly oscillation differ from zero for only a few moderate-degree coefficients at a significance level of 5% .The bands in Fig. 11i at orders m = 15 and m = 30 could be due to resonant orbital perturbations.The amplitudes of the disturbances are superimposed in a resonance-like manner due to a strict or approximate commensurability between the mean motion of the satellite and the rotation of the Earth.In these cases, the ratio of the angular velocities is an integer (Kaula 2000).
The results of the parameter tests strengthen the assumption that the functional models are quite suitable for describing the temporal variation in the low-degree coefficients.Despite these results, all parameters are taken into account for the synthesis of the new models, even for the higher degree coefficients.This is done because the Fourier analysis has shown that periods of three and 6 months are also contained in the data and in order to additionally model part of the random signal contained.
The relationships f 2 and f 3 , which differ in the quarterly period P 3 , give similar results.The relationships f 4 and f 5 , where the quadratic and cubic terms are replaced by the nutation period P L , lead to the lowest weighted sum of squared residuals on average (cf.Table 5).Therefore, f 3 and f 4 are used for further comparisons.
In the present work, the authors did not have a fully occupied variance-covariance matrix at their disposal.Hence, the respective weight matrices were generated as diagonal matrices from the available accuracy information for the Stokes coefficients download from (ICGEM 2023).The formal errors described in this way are, as is well known (Dahle et al. 2019a, b), too optimistic because the temporal correlations are not taken into account.Thus, using the formal errors yields an increased significance compared with more realistic errors.The used (too optimistic) accuracies lead to a higher probability of error and thus to a higher probability of an error of first kind.This means that H 0 is rejected even though the estimated parameter is actually not significant.In other words: as seen from ( 27) the test statistic T G will become larger if 0 is to optimistic, and this leads to a earlier significance.The question of realistic a priori accuracies is a very important research question.

Verification by comparison with monthly solutions
As explained at the end of the Sect.5.3, approaches f 3 and f 4 are now compared with each other.They differ in the additional quadratic and cubic terms in the polynomial approach.Therefore, as a measure of the external precision, a comparison with existing monthly solutions is carried out with these approaches.Figure 14 shows the effect on geoid undulation for October 2010 with the approaches f 3 and f 4 .For the determina- tion of the time-dependent parameters within the framework of the estimation denoted by n V = 160 the GRACE model for October 2010 did not enter the estimation process.
The associated statistical parameters are listed in Table 7.This shows that f 4 leads to a smaller deviation on average, but has a larger minimum than f 3 .The largest differences are found in South America around the equator in the Amazon Basin.This region is subject to large seasonal fluctuations due to seasonal change in the total water storage (TWS) (Tapley et al. 2004(Tapley et al. , 2019)).The year 2010 saw extreme droughts in the Amazon Basin compared to previous years, peaking in October.The reasons for this were the occurrence of the El Niño phenomenon in connection with warming of the tropical Atlantic (Marengo and Espinoza 2016).Since these phenomena do not occur at regular intervals, they cannot be modeled with f 3 and f 4 .Therefore they are visible as poten- tial differences transformed to geoid undulations applying Eq. ( 14).The use of wavelets (Freeden and Schneider 1998;Freeden and Schreiner 2022) could possibly show that the frequency spectrum is stable, but its amplitudes vary over time.
The comparison in the frequency domain using the degree variances (Fig. 15) shows only small differences between f 3 and f 4 for the low degree coefficients.The error degree variances are shown in Fig. 16.They are also very similar for the different approaches, but differ from the error degree variances of the GRACE monthly solution from October 2010: The coefficients up to about degree n = 22 are associated with a greater uncertainty for the model, determined according to f 3 and f 4 than the coefficients of the GRACE monthly solution.On the other hand, the Stokes coefficients with higher degree are associated with less noise for the calculated model.There is also a kink in the error degree variances at around degree n ≈ 40 (cf.global test, Appendix 2).
In addition, the differences in geoid undulations for January 2010 are on average lower than for October 2010 and vary symmetrically around zero.As with October 2010, the models determined using approach f 4 lead to a smaller mean difference, but show a higher maximum than f 3 , which is also located in the Amazon Basin.

Assesment of the developed models against GRACE-FO
In addition to the comparison with GRACE monthly solutions to assess the quality of the calculated time-dependent parameters, a comparison with a GRACE-FO month is carried out.The GRACE comparison shows how good the agreement is within the observation  period.Observations from GRACE-FO can also be used to check how the time-dependent parameters behave outside the period of the GRACE mission (2002-2017).
In the following, a model for January 2019 determined from approaches f 3 and f 4 is checked for consistency with the corresponding monthly solution of the GRACE-FO satellite mission.Figure 17 shows the effects on geoid undulation N for the difference model.For f 3 the agreement is slightly worse than in October 2010 (Table 7, GRACE), which is probably due to the fact that it is an extrapolation to 2019.Significant differences can given in the statistical values (cf.Table 8 when calculating a new model based on f 4 with a polynomial approach.Figure 17b and the associated degree variances in Fig. 18b show that there are large differences, especially in the low-degree coefficients.These Stokes coefficients are subject to greater temporal variations than coefficients of higher degree.The extrapolation leads to larger deviations when using the approach f 4 due to the third-degree polynomial, whereas it shows a good fit for interpolations.The reason for this lies in the properties of third-degree polynomials: The additional consideration of a polynomial for the periodic processes means a better local adjustment of the data to the functional context within the period from 2002 to 2017.In this way, for example, a slightly longer periodic signal can be intercepted.However, the polynomial is unsuitable for the extrapolation, as this tends towards infinity over time and is therefore not suitable as a model for describing the displacement of water masses.The error degree variances of the difference model are also larger for f 4 (Fig. 19b) than for f 3 (Fig. 19a).

Comparison to EIGEN-6S4(v2)
To assess the agreement between the time-dependent geopotential model from GRACE monthly solutions and the EIGEN-6S4(v2), models for the month of January 2010 are generated using the time-dependent parameters.Since the parameters of the EIGEN model are determined for one year in each case, in addition to the adjustment over all n = 161 obser- vations, adjustments are also carried out in windows of one year.The monthly solutions from one year and the solutions from December of the previous year and January of the following year are included in the adjustment ( n E = 14 ).In the EIGEN model, six param- eters are given according to Eq. ( 7) to model each time dependent coefficients up to degree and order n max = 80 .Therefore, in addition to approaches f 3 and f 4 , a comparison is also carried out using the EIGEN-like f 2 .

Adjustment in annual windows
First, the results of the comparison with n E = 14 observations are presented.The Fig. 20a, b show the geoid undulation for the respective difference model depending on the functional approach f 2 and f 4 .Above all, finely structured differences can be seen: The EIGEN model takes into account time-dependent parameters up to degree and order n = 80 and then uses mean values to describe the Stokes coefficients.The GRACE model, on the other hand, additionally models the time dependence for the coefficients up to degree n max = 96 .This results in deviations between the models, e.g.around 120 • east longitude and 60 • south latitude.There are only small differences between the individual analytical approaches, with f 2 leading to the smallest minimum and maximum deviation.This model- ling with six parameters (cf.Table 2) corresponds to the modelling of the EIGEN-6S4(v2) model.The statistical parameters for the differences in geoid undulation listed in Table 9 show that f 4 leads to the largest deviations.
The degree variances of the two models (Fig. 21a, b) show the differences in the definition of the coefficients C 1,0 , C 1,1 and S 1,1 .In all GRACE monthly solutions the potential coefficients of degree n = 1 are set to zero.This means that the origin of the coordinate system is linked to the centre of mass of the Earth.At GRACE, all developments are performed relative to the Earth's centre of gravity.The EIGEN model, on the other hand, also takes into account time-dependent parameters for the coefficients of degree n = 1 .Changes in this term indicate scaled variations of the centre of mass with respect to the origin due to mass shifts (Peters 2007).The degree variances of the difference model decrease up to degree n = 80 and then increase again due to the time-dependent parameters up to degree n max = 96 for the GRACE model.The error degree variances (Fig. 22a, b) are up to degree n = 20 higher for the EIGEN model than for the GRACE model.The higher degree coef- ficients, on the other hand, are subject to lower inaccuracies, and the error degree variances are lower than for the GRACE model.

Adjustment over all observations
The panels in Figs.23 and 24 show the results for the comparison of the EIGEN model with the time-dependent GRACE model in the space and frequency domain for January 2010, respectively.The time-dependent GRACE parameters are thereby calculated from n = 161 observations.For the EIGEN model, the parameters from epoch January 2009 to February 2010 are used.In contrast to the results for n E = 14 where f 2 gives the best agree- ment, f 4 shows the smallest std from n = 161 (cf.Table 9).The scatter is also less than for f 2 and f 3 , although the range of the residuals is larger.On average, the deviations are some- what larger than for n E = 14 observations.For n = 161 observed monthly solutions, the differences in geoid undulations show sim- ilar structures, such as in the Pacific near Indonesia, and deviations, such as in Greenland, in North and South America, which are shown in the panels of Fig. 20.These deviations in Greenland, North and South America can also be seen in a comparison with the GRACE monthly solution from January 2010 in the space domain (Fig. 26) and in the frequency domain (Figs. 27,28).This shows that there are more changes in these areas than in previous years, which is why they are modeled poorly.This can be better intercepted by determining the time-dependent parameters for each year (Fig. 20a, b).
The comparison in the frequency domain shows that the difference between the determination of the parameters over a year ( n E = 14 ) and over 15 years ( n = 161 ) lies mainly in the coefficients up to degree n = 50 , almost regardless of the functional model ( f 2 , f 3 or f 4 ).The Fig. 25a, b also show the changes for the accuracies: The error degree variances of the GRACE model are lower compared to the error degree variances in the Fig. 21a, b, because with n = 161 observations the redundancy is greater and the time-dependent parameters for the GRACE model can be determined with higher accuracy.The accuracies for the coefficients from degree n = 80 on are comparable for the EIGEN and GRACE models.

Summary and outlook
Based on monthly GRACE solutions, time-dependent geopotential models are determined.For the potential coefficients up to degree and order n max = 96 , time-dependent parameters are determined using five different analytical approaches in the framework of the Gauss-Markov model.
Long-term changes in the gravitational field are taken into account using a linear trend or a third-degree polynomial.The numerical investigations have shown that the latter can be modeled very well using the f 5 approach due to the long-period effect of lunar nuta- tion with a period of 18.6 years.Periodic variations are described mathematically by periodical functions with annual, semi-annual and a quarterly periods, as suggested by a initial FFT for a few selected coefficients of the monthly potential models.To check the quality of the models, global tests and parameter tests are carried out for the estimation results.In addition, model variants are calculated for selected months and are carefully compared to GRACE, GRACE-FO and EIGEN-6S4 monthly solutions in the space and frequency domain.
The statistically significant periods found in the analysis are in accordance with the expected physically implied effects (cf.Table 1).
It should be noted that two types of basis functions (polynomial and trigonometric) were used to model the time variations of the potential coefficients, which are not orthogonal to each other.This can mean that a very long-period component can also be included in the determination of the polynomial trend parameters and vice versa.The correlation between the estimated parameters are studied for a few Stokes coefficients.The estimated correlation coefficients show very high correlation (nearly ±1 ) for the trend parameters.This might be an effect of the choice of t 0 .Only a weak correlation was found between the parameter sets of the trend and periodically parameters.
The models according to approach f 3 with eight time-dependent parameters (linear trend, six amplitudes, the later ones can be transformed to three amplitudes and three phase parameters) and approach f 4 with ten parameters (third degree polynomial, six amplitudes) lead to the lowest sum of squared residuals and are therefore used for further comparisons.The monthly GRACE comparison for the month of October 2010, when there was extreme drought in the Amazon Basin, leads to mean differences in geoid undulations of N 3 = −0.24mm for f 3 and N 4 = −0.11mm for f 4 .Another monthly comparison without any special occurrences (January 2010) shows significantly smaller mean deviations, which are slightly smaller for f 4 than for f 3 .Thus, a third degree polynomial can fit the trend within the time series of the coefficients slightly better than applying a linear trend alone.Of course, adding additional parameters, here the cubic parameters, leads to an increase in the degree of freedom of the functional model.As a result, there can be a reduction in the residuals, which is desirable.The parameter test prevents any model expansion because its significance is examined.In the analyzes carried out it was significantly shown that the cubic model approach is equivalent to the approach f 5 , which includes the lunar period instead of the cubic terms.
As mentioned in Göttl et al. (2019) the DDK5 filter with 400 km Gaussian smoothing (Kusche 2007;Kusche et al. 2009) does not only smooth or reduce the striping effects in the GRACE date, or noise, but of course also the high frequency signal.This might be a reason why the Stokes coefficients are not significantly time dependent for higher degree than n ≈ 40.
In the context of determining time-dependent parameters for spherical Stokes coefficients, there are some possibilities for future research fields.If, in addition to the currently available accuracy information for the Stokes coefficients, their correlations are also made available, a study on their influence on the estimation of the time dependent model parameters would be of great interest.
The determination of the parameters for describing the time variability of the coefficients on the basis of all monthly solutions has shown that a worse fitting can occur for individual years.Further investigations based on the estimation of the parameters for only one year would be appropriate.A direct comparison with EIGEN-6S4(v2) could then be made.It can be expected that the estimation in windows of one year will improve the consistency with GRACE or GRACE-FO months, and also the EIGEN-6S4 geopotential model.This allows to observe how the parameters and their precision change over time.The use of wavelets could possibly show that the frequency spectrum is stable, but its amplitudes vary over time.As an alternative model approach to describe the time-varying gravity field, a decomposition using the independent component analysis (ICA) could also be investigated.The ICA was successfully carried out by Forootan et al. (2020) on the functional of the time-varying part of the gravity field in the form of the time-varying global terrestrial water storage changes.
In addition, it was revealed that the time-dependent parameters are not significant for every coefficient.Therefore, coefficients of different degree could be better described via different functional models.As a next step, a Fourier analysis can be applied to the residuals of the estimations.With the result, the functional approach can be adjusted until the Fourier analysis only shows white noise.
Further investigations based on the unfiltered monthly solutions could be carried out in the frequency domain separately for each degree.This could be an advantage over general smoothing of the entire signal.
Because strong earthquakes cause abrupt changes in the gravity field, incorporating an offset, as in the EIGEN models, could also result in better agreement with existing models.
with E(⋅) denoting the expectation vector and v the residual vector.The stochastic model is given by the (regular) variance-covariance matrix of the observations.It contains information on the uncertainty of the observations in terms of variances and covariances.If an a priori variance factor 2 0 is extracted, the cofactor matrix includes the uncertainty relations: The weight matrix P equals the inverse cofactor matrix occurring in Eq. ( 20) (Koch 1999): From Eqs. ( 18) and ( 19), the unique solution for the estimated parameters is given by: provided that the inverse of the quadratic normal equation matrix N = A T PA does exist.

The vector of estimated residuals is
The cofactor matrix for the estimated parameters Q xx results from variance propagation: For the statistical hypothesis tests, the a posteriori variance factor σ2 0 is computed as an estimate of the a priori variance factor 2 0 .The a posteriori variance factor is composed of (20) Table 10 Availability of the GRACE and GRACE-FO monthly solutions by GFZ: "-" no solution available; "(✓ )" solution is not used, because the data used to create the solution of this month do not cover the whole month The entries refer to a query in the ICGEM (2023) database from October 2023 the weighted sum of squared residuals Ω and the redundancy r of the estimation problem (Koch 1999;Niemeier 2008): Inserting ( 20) and ( 21) into (25) it is obvious that the weighted sum of the squared residuals does not depend on the choice of 2 0 .

Statistical hypothesis tests to assess estimation results
In the context of the estimation problem, statistical tests are applied as validity check and quality assurance of the results.Two types of statistical hypothesis tests are used: a test for consistency of model and observations (global test) and a test of estimated parameters.
For testing any hypothesis, a procedure according to Koch (1999) is recommended: in the first step, the null hypothesis H 0 and the alternative hypothesis H A are set up.The sig- nificance level is defined.In the next step, a suitable test statistic T is calculated.According to the distribution of the test statistic T and the significance level, the test's range of acceptance or rejection for H 0 is defined.Finally, the test decision is made.Based on this scheme for statistically testing hypotheses, a test is applied to check the (global) consistency of the model and the observations.H 0 implies that the a posteriori variance factor σ2 0 as estimated value is not significantly different from the a priori variance factor 2 0 as theoretical value.The alternative hypothesis states that the σ2 0 differs significantly from 2 0 .The test is performed one-sided: According to Eq. ( 20) an a priori variance factor 2 0 can be extracted from the variancecovariance matrix.The chosen factor acts simply as a scaling and can help to stabilize the numerical behavior of the weighting matrix P .The test statistic T G contains the ratio between σ2 0 and 2 0 (see Eq. ( 27)) and is therefore independent on the choice of 2 0 .The test statistic T G is (under the H 0 hypothesis) chi-square distributed with f = r degrees of freedom: T G ∼ 2 f .The number of degrees of freedom f corresponds to the redundancy r of the estimation problem defined in (17) (Koch 1999).Taking the relation (25) into account one finds: The null hypothesis H 0 is rejected with a significance level of if T G is greater than the corresponding quantile of the chi-square distribution.According to Niemeier (2008), there are two possible reasons for rejecting the null hypothesis: firstly, it is possible that 2 0 does not fit the existing measurement conditions, and secondly, a rejection can be due to inconsistencies in the functional model or in the stochastic model.This test will be denoted as global test throughout the paper.
In contrast to this test, which applies to the functional and stochastic model, parameter tests are used to test null hypotheses concerning the estimated parameters xj .This is a one- sided test with null hypothesis H 0 and alternative hypothesis H A , which examines whether ( the expected value of the estimated parameter is significantly different from a value x 0 , (Neyman and Pearson 1933): The parameters x 0,j ∋ {c q , a i , b i } from the functional approach ( 19) are tested on their sig- nificance.Therefore the parameters x 0,j were tested to determine if they significantly differ from the expected value x 0,j = 0 .However, it is important to note that in general x 0 is not always set to zero.For example, when testing a scale factor, x 0 could be set to 1.
The test statistics depend on the parameter xj and the associated a priori standard deviations  xj and a posteriori σx j , respectively.The square of the resulting test statistics is Fisher distributed, with the degree of freedom ∞ and r, respectively.The limits for the ranges of acceptance or rejection of the null hypothesis are taken from the corresponding quantile of the F-distribution: Analogously to the one-dimensional parameter test, an m-dimensional test can be performed: The model expectation vector = (x 0,1 , x 0,2 , … , x 0,u ) T has the dimension u = m.
The used in-and output parameters in the framework of the estimation of the time dependent Stokes coefficients C nm (t) and S nm (t) and the applied tests on their signifi- cance are briefly listed in Tables 11 and 12.

Fig. 1
Fig. 1 Comparison of the 161 GRACE monthly solutions (2002-2017) in the frequency domain based on degree and error degree variances.In panels b and c a few conspicuous monthly solutions are visible

Fig. 8
Fig. 8 Test statistic T G for the different analytical approaches ( f 1 to f 5 ) used to determine the time-dependent geopotential model from GRACE monthly solutions (DDK5) according to Table 2.The values of the test statistic are unit-less.Note the logarithmic scale

Fig. 9
Fig. 9 Results of the global tests for the different analytical approaches ( f 1 to f 5 ) used to determine the time-dependent geopotential model from GRACE monthly solutions (DDK5) according to Table 2 (black: H 0 rejected → significant, yellow: H 0 accepted → not significant).The global test for all f i look quite similar and indicate that H 0 is rejected for n <≈ 40

Fig. 10 Fig. 11
Fig. 10 Results of the global tests for the analytical approach ( f 4 ) used to determine the time-dependent geopotential model from GRACE monthly solutions based on DDK2, DDK5, DDK8 and unfiltered data.(black:H 0 rejected → significant, yellow: H 0 accepted → not significant)

Fig. 12 Fig. 13
Fig. 12 Results of the parametric tests for f 5 based on the DDK5 filtered data.(black: H 0 rejected → significant, yellow: H 0 accepted → not significant) Fig. 14 Comparison of f 3 and f 4 with the GRACE monthly solution in the space domain by geoid undulation differences (October 2010).Units are (mm) Fig. 15 Comparison of f 3 and f 4 with the GRACE monthly solution in the frequency domain by degree variances (October 2010).Units are (m 4 s −4 ) Fig. 16 Comparison of f 3 and f 4 with the GRACE monthly solution in the frequency domain by error degree variances (October 2010).Units are (m 4 s −4 )

Fig. 17 Fig. 18
Fig. 17 Comparison of f 3 and f 4 with the GRACE-FO monthly solution in the space domain by geoid undulation differences (January 2019).Units are (mm)

Fig. 19
Fig. 19 Comparison of f 3 and f 4 with the GRACE-FO monthly solution in the frequency domain by error degree variances (January 2019).Units are (m 4 s −4 )

Fig. 26 Fig. 27
Fig. 26 Comparison of f 3 and f 4 with the GRACE monthly solution in the space domain by geoid undulation differences (January 2010).Units are (mm)

Table 2
Compilation

Table 3
Monthly solutions are available on ICGEM (2023) with different levels of smoothingThe resolution in the space domain is restricted by the degree of smoothing

Table 4
Constants and parameters entering the determination of the time-variable description of the gravitational potential The Stokes coefficients and their standard deviations are available each month (epoch t) with an maximum degree and order n max = 96 .The used constants , R and n max are taken from the header of the coefficients files (Dahle et al. 2018b) max + 1)(n max + 2)∕2 Number of Stokes coefficients ( k = 4753 for n max = 96) Total number of adjustments required are set up as diagonal matrices.The functional model is formulated based on Eqs. (

Table 6
Critical values of the test statistics based on DDK5

Table 2
(black: H 0 rejected → significant, yellow: H 0 accepted → not significant).The global test for all f i look quite similar and indicate that H 0 is rejected for n <≈ 40

Table 7
Statistical values:Comparison with the GRACE monthly solution in the space domain by differences in geoid undulations

Table 8
Statistical values:Comparison with the GRACE-FO monthly solution in space domain by differences in geoid undulations for January 2019

Table 9
Design matrix of dimension n × u j for the j− th functional approach with j ∈ {1, … , 5} Cofactor matrix of the observations of dimension n × n × k j Vector of observations of dimension n × k 0 A priori standard deviation