GCL-Mascon2024: a novel satellite gravimetry mascon solution using the short-arc approach

This paper reports on an innovative mass concentration (mascon) solution obtained with the short-arc approach, named “GCL-Mascon2024”, for estimating spatially enhanced mass variations on the Earth's surface by analyzing K- and Ka-band ranging satellite-to-satellite tracking data collected by the Gravity Recovery And Climate Experiment (GRACE) mission. Compared to contemporary GRACE mascon solutions, this contribution has three notable and distinct features: first, this solution recovery process incorporates frequency-dependent data-weighting techniques to reduce the influence of low-frequency noise in observations. Second, this solution uses variably shaped mascon geometry with physical constraints such as coastline and basin boundary geometries to more accurately capture temporal gravity signals while minimizing signal leakage. Finally, we employ a solution regularization scheme that integrates climate factors and cryospheric elevation models to alleviate the ill-posed nature of the GRACE mascon inversion problem. Our research has led to the following conclusions: (a) GCL-Mascon2024 mass anomaly estimates from GRACE data show strong agreement with the (Release) RL06 versions of mascon solutions (GSFC, CSR, JPL) in both spatial and temporal domains; (b) in Greenland and global hydrologic basins, the correlation coefficients of estimated mass changes between GCL-Mascon2024 and other RL06 mascon solutions exceed 95.0 %, with comparable amplitudes, and, especially over non-humid river basins, the GCL-Mascon2024 suppresses random noise by 27.8 % compared to contemporary mascon products; and (c) in desert regions, the analysis of residuals calculated after removing the climatological components from the mass variations indicates that the GCL-Mascon2024 solution achieves noise reductions of over 29.3 % as compared to the GSFC and CSR RL06 mascon solutions. The GCL-Mascon2024 gravity field solution (Yan and Ran, 2025) is available at https://doi.org/10.5281/zenodo.15525467.

1 Introduction

Comprehending the Earth as a dynamic system relies heavily on our knowledge of its gravity field, mass variations induced by fluid layers, and geophysical or climatic processes (e.g., Wahr et al., 1998; Pail et al., 2015). Over the last 2 decades, significant achievements have been made through the availability of observations collected by satellite gravimetry missions, such as the Gravity Recovery And Climate Experiment (GRACE) (Tapley et al., 2004; Tapley et al., 2019) and its successor, GRACE Follow-On (GRACE-FO) (Flechtner et al., 2016; Landerer et al., 2020). These satellite gravity missions have not only enhanced our understanding of temporal variations in the Earth's gravity field but also played a crucial role in advancing various disciplines, including glaciology, hydrology, geophysics, oceanography, and atmosphere and climate science (e.g., Han et al., 2006a; Chen et al., 2009; Rignot et al., 2011; Jacob et al., 2012; Rodell et al., 2018).

Gravity field variations expressed in spherical harmonics have been extensively and widely employed in satellite geodesy for decades (Chen et al., 2022). However, certain limitations persist in the application of spherical harmonic solutions from GRACE and GRACE-FO data, including the presence of north–south “stripes” (Swenson and Wahr, 2006), as well as signal leakage (Kusche et al., 2009), particularly in regions adjacent to the land–sea boundary. The main reasons for the above problems are temporal aliasing (Wiese et al., 2011b) and the design of the satellite orbits and tracking systems (Wiese et al., 2011a), including inclination, altitude, inter-satellite distance, and a co-planar low–low satellite-to-satellite tracking system. Conventional approaches involve the removal of stripes through empirical smoothing (e.g., Wahr et al., 1998), de-striping (e.g., Swenson and Wahr, 2006), or regularization techniques (e.g., Save et al., 2012). It is important to note that, although these methods are largely effective in preserving signals and suppressing noise, the elimination of stripes also results in a reduction in the genuine geophysical signals (e.g., Han et al., 2005; Yi and Sneeuw, 2022; Zhou et al., 2023). Moreover, the efficacy of destriping is highly dependent on the characteristics of the signals, including their size, shape, and orientation (Watkins et al., 2015). It is worth mentioning that the impact of aliasing errors can be mitigated by combining gravity satellite formations within optimal constellation configurations (Yan et al., 2024) or by recovering the temporal gravity field at a higher temporal resolution (Yan et al., 2023).

Alternatively, mass concentration (mascon) solutions can be utilized to model the temporal gravity field. This technique was initially introduced by Muller and Sjogren (1968) in their efforts to develop a model for the static gravity field of the moon. Thereafter, mascon solutions utilizing GRACE level-1B data were initially conducted in a regional context (e.g., Rowlands et al., 2005; Luthcke et al., 2006) and were subsequently extended to encompass diverse global parameterizations (e.g., Luthcke et al., 2013; Watkins et al., 2015; Save et al., 2016; Allgeyer et al., 2022). Besides this, some attempts have been made to enhance mascon solutions' spatial (e.g., Loomis et al., 2021) or temporal resolutions (e.g., Croteau et al., 2020). Additionally, to mitigate the computational complexity, alternative variants of the mascon approach have been put forward, which utilize monthly sets of spherical harmonic coefficients (SHCs, i.e., level-2 data) as input (e.g., Forsberg and Reeh, 2006; Baur and Sneeuw, 2011; Schrama and Wouters, 2011). Numerous recent publications have used mascon solutions released by responsible agencies, including the NASA (National Aeronautics and Space Administration) Goddard Space Flight Center (GSFC), the NASA Jet Propulsion Laboratory (JPL), and the University of Texas at Austin Center for Space Research (CSR) in the United States. The mascon solution released by JPL (JPL RL06 mascon) utilizes explicit partial derivatives with analytical expressions for the mascons to establish the relationship between inter-satellite range rate measurements and individual mascons (Wiese et al., 2018), whereas the latest variants of GSFC mascon solutions (GSFC RL06 mascon) and CSR mascon solutions (CSR RL06 mascon) are characterized by a finite series of spherical harmonic functions, with the corresponding partial derivatives being computed using the chain rule (Loomis et al., 2019; Save, 2020). These GRACE and GRACE-FO gravimetry data-processing centers also offer visualization tools for their mascon products, facilitating analysis and comparison of the latest mascon solutions, as well as generating time series data for specific regions.

Various methods, including the dynamic approach (e.g., Kvas et al., 2019), the short-arc approach (e.g., Mayer-Gürr, 2008), the celestial mechanics approach (e.g., Beutler et al., 2010), the energy balance approach (e.g., Han et al., 2006b), and the acceleration approach (e.g., Ditmar and van der Sluijs, 2004), play a vital role in modeling the temporal gravity field from level-1B satellite gravimetry data. To date, most publicly available global mascon products based on level-1B data commonly rely on longer arcs (e.g., 24 hr ones). This includes the mascon solutions recovered using the dynamic approach by GSFC (Loomis et al., 2019), CSR (Save, 2020), and JPL (Watkins et al., 2015), as well as the mascon solution by the Australian National University (ANU), which utilizes the celestial mechanics approach (Allgeyer et al., 2022; Tregoning et al., 2022; McGirr et al., 2023). This study represents the first application of the short-arc approach to recover the global mascon solution. A distinguishing feature of this methodology compared to other conventional approaches lies in its substantially reduced arc length integration interval (Mayer-Gürr et al., 2005; Mayer-Gürr, 2008). The temporal gravity field based on the short-arc approach exhibits enhanced stability and superior accuracy owing to the substantially reduced condition number of the normal equation system (Chen et al., 2015).

Frequency-dependent noise in GRACE measurements significantly limits GRACE from reaching the pre-launch baseline accuracy; thus, modeling this noise is a critical aspect for improving the accuracy of temporal gravity field recovery. In the context of spherical harmonic coefficient solutions, the impact of frequency-dependent noise in observations is typically accounted for by introducing empirical parameters (Liu et al., 2010; Zhao et al., 2011) to absorb errors or by using frequency-dependent data-weighting (FDDW) techniques (Klees et al., 2003; Ditmar et al., 2007). However, the potential of suppressing frequency-dependent errors in mascon modeling with the FDDW technique remains largely unexplored.

The Geodesy and Cryosphere Laboratory (GCL) from the Southern University of Science and Technology has released a new series of mascon solutions (hereafter referred to as GCL-Mascon2024) using the short-arc approach and FDDW, as well as advanced regularization schemes. These mascon solutions incorporate pertinent physical constraints to estimate global mass variations directly from inter-satellite range rate measurements. To alleviate the effects of errors introduced by signal leakage, the GCL-Mascon2024 solution employs a strategy that involves segmenting the mascon shape based on land–sea boundaries and the boundaries of distinct hydrologic basins. Subsequently, this paper aims to investigate the impact of selecting arc length and accelerometer calibration parameters in the short-arc approach on the mascon solutions while also providing a quantitative evaluation of the GCL-Mascon2024 solution.

The article is organized as follows. Section 2 describes the methodology for recovering global mascon solutions with the short-arc approach. Section 3 discusses the parameter determination in global mascon solutions using the short-arc approach. Section 4 evaluates the scientific results of real data processing with the proposed approach. Finally, Sect. 5 provides the main conclusions. Section 6 provides detailed information and links for accessing the dataset utilized in this study, along with the GCL-Mascon2024 solution released in this work.

2 Methodology

Building upon the earlier studies by Ran et al. (2018) and Ran et al. (2021), we propose a new mascon approach recovered from GRACE level-1B tracking data based on the short-arc approach. The primary distinction between GCL-Mascon2024 and the aforementioned mascon solutions lies in the type of exploited input data (i.e., level 1B vs. level 2). The mascon solutions released by Ran et al. (2018, 2021) are based on spherical harmonic coefficients and cover only mass anomalies over Greenland. The GCL-Mascon2024 solution is a series of global mascons with analytical partial derivatives. In other words, we establish a direct relationship between the mass variations of mascons and the inter-satellite measurements. Section 2.1 elaborates on the utilized functional model, which links GRACE level-1B data to mascon solutions. Section 2.2 outlines the strategy for defining mascon geometry during the data inversion process. Section 2.3 describes the background force models and input data employed to recover the GCL-Mascon2024 solution. Section 2.4 explains the suppression of frequency-dependent errors by using the FDDW technique. Finally, the advanced spatial constraints exploited in the inversion procedure are presented in Sect. 2.5.

2.1 Mathematical formulation

A satellite in orbit around the Earth is subject to gravitational forces, which are governed by Newton's law of universal gravitation. The temporal gravity field can be modeled as a series of N mascons, with the surface mass density (mass per unit area) of mascon M_i being represented by ${\mathit{\rho }}_i(i=\mathrm{1},\mathrm{2},\mathrm{\dots },N)$. When the satellite is at measurement point p, the gravitational forces f_p exerted on the satellite by the mass variations of the Earth's surface can be expressed as

$$\begin{array}{}\text{(1)}& {\mathit{f}}_p=\mathrm{G}\sum _{i=\mathrm{1}}^N{\mathit{\rho }}_i\underset{M_i}{\int }{\displaystyle \frac{{\widehat{\mathit{d}}}_p\cdot \mathrm{d}s}{{\left(l_p\right)}^{\mathrm{2}}}}=G\sum _{i=\mathrm{1}}^N{\mathit{\rho }}_i\cdot {\widehat{\mathbf{I}}}_{i,p}.\end{array}$$

Here, G is the universal gravitational constant; ${\widehat{\mathit{d}}}_p$ is the unit vector directed from the satellite measurement point toward the surface mass; l_p is the distance between the satellite measurement point p and an integration point on the mascon; and ${\widehat{I}}_{i,p}$ is a vector pointing from the satellite measurement point p to the given mascon M_i, which is calculated using numerical integration. To that end, we utilize a composed Newton–Cotes formula (Gonzalez, 2010) applied to the Fibonacci nodes, i.e., the Fibonacci nodes as integration points, as mentioned above. By defining the surface area and the number of the Fibonacci nodes of mascon M_i as S_i and K_i, we can calculate ${\widehat{I}}_{i,p}$ as

$$\begin{array}{}\text{(2)}& {\widehat{\mathbf{I}}}_{i,p}\approx \sum _{j=\mathrm{1}}^{K_i}{\displaystyle \frac{S_i}{K_i\cdot {\left(l_{i,j,p}\right)}^{\mathrm{2}}}}\cdot {\widehat{\mathit{d}}}_{i,j,p},\end{array}$$

where $l_{i,j,p}$ represents the distance between a Fibonacci point j located in the mascon M_i and the satellite measurement point p, and ${\widehat{d}}_{i,j,p}$ is a unit vector pointing from the satellite measurement point p to a Fibonacci point j located in the mascon M_i.

Then,

$$\begin{array}{}\text{(3)}& {\mathit{f}}_p=\underset{\mathit{x}}{\underbrace{\sum _{i=\mathrm{1}}^N}}{\mathit{\rho }}_i\cdot \underset{G_p}{\underbrace{\mathrm{G}{\displaystyle \frac{S_i}{K_i}}\cdot \sum _{j=\mathrm{1}}^{K_i}{\displaystyle \frac{{\widehat{\mathit{d}}}_{i,j,p}}{{\left(l_{i,j,p}\right)}^{\mathrm{2}}}}}}.\end{array}$$

Combining G_p over multiple positions or epochs within an arc yields the matrix G, which is used in the observation model (Mayer-Gürr, 2008) with orbit and range rate measurements as observation types.

2.2 Parameterization

The choice of an appropriate mascon partitioning strategy is crucial for mitigating noise amplification during the data inversion process (Ran et al., 2018). In this study, the selection of mascon geometry is based on incorporating pertinent physical constraints, such as the geometry of the coastlines and basin boundaries. The definitions of these basin boundaries are derived from Scanlon et al. (2018). Regarding the aforementioned parameterization, the primary assumption is that there is no signal correlation between mascons located in different basin systems (Ran et al., 2021), meaning that basins do not share mascons with their neighboring basins to reduce signal leakage between the corresponding basins.

In the GCL-Mascon2024 processing scenario, the estimated monthly mascon solution has a spatial resolution of about 300 × 300 km and 400 × 400 km on land and in the ocean, respectively. The total number of mascons is 4069, with 1852 terrestrial mascons and 2217 ocean mascons. Figure 1 provides the mascon partitioning of GCL-Mascon2024. It is important to note that the mascons located within the basins and coastal regions are defined in accordance with the boundary geometry. The numerical integration points, as discussed in Sect. 2.1, are distributed on a Fibonacci grid with an average spacing of 10 km, requiring the generation of approximately 5.1 million Fibonacci grid points for global coverage. Parallel message-passing-interface (MPI) computing is used to increase computational efficiency.

Figure 1 Mascon partitioning of GCL-Mascon2024 solution.

2.3 Background force models and input data

Using the aforementioned methodology and mascon partitioning strategy, we have produced a time series of mascon solutions (GCL-Mascon2024) from the GRACE level-1B data covering the time period from January 2003 to December 2015. Here, we concisely introduce the background force models and input tracking data.

Table 1 provides an overview of the background force models, which encompass various components, including Earth's static gravity field, third-body attractions, solid Earth (pole) tides, ocean (pole) tides, atmospheric tides, atmospheric and oceanic dealiasing effects, and general relativistic correction. In addition to the background force models mentioned above, these additional force models and corrections are discussed in detail below, specifically (i) the elastic response of the solid Earth to mass transport at the Earth's surface, (ii) glacial isostatic adjustments (GIAs), (iii) Earth ellipsoidal corrections, (iv) low-degree term corrections, and (v) GAD corrections. Following a standardized processing workflow (Watkins et al., 2015; Save et al., 2016; Loomis et al., 2019; Tregoning et al., 2022), the uncorrected mascon solutions (i.e., MASCON_uncorrected, which we will return to in Sect. 2.5) are systematically integrated with the aforementioned corrected components to generate corrected mascon grids. The formula to generate the corrected mascon grid is

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathrm{MASCON}}_{\mathrm{corrected}}={\mathrm{MASCON}}_{\mathrm{uncorrected}}\\ & -{\mathrm{MASCON}}_{C_{\mathrm{20}}}+{\mathrm{SLR}}_{C_{\mathrm{20}}}+\mathrm{DEG}\mathrm{1}-\mathrm{GIA}+\mathrm{GAD}.\end{array}\end{array}$$

Table 1 also lists the input used in mascon recovery, including nongravitational accelerations, satellite attitudes, reduced-dynamic orbits, kinematic orbits, and K-band range rate measurements. The level-1B data used in the mascon recovery are mainly from JPL, e.g., ACC1B, SCA1B, GNV1B, and KBR1B. Additionally, the kinematic orbit product released by the Graz University of Technology (Strasser et al., 2019) was used in the GCL-Mascon2024 recovery framework.

Table 1 Summary of background force models and data used in GRACE mascon recovery.

* Data missing, indicating that the data or strategy are unavailable.

2.3.1 Earth's elastic response

The solid Earth is not perfectly rigid but exhibits some elastic response to surface loads (Boy and Chao, 2005). Here, we estimate the effect of surface load or surface mass changes based on the elastic loading theory of a spherical Maxwell Earth, as formulated by Wahr et al. (1998), who used load Love numbers (represented as k_l) to quantify Earth's elastic deformation.

In this study, the temporal gravity field model released by the Institute of Geodesy of the Graz University of Technology (ITSG-Grace2018 (Kvas et al., 2019)) is used as the signal source to compute the Earth's elastic deformation. Because this model is represented in terms of unfiltered spherical harmonic coefficients, there exist north–south stripes and high-frequency noise in the spatial domain. Thus, post-processing in the form of the DDK4 filter (Kusche et al., 2009) is used to mitigate these issues. The elastic deformations induced by the filtered ITSG-Grace2018 solutions are incorporated into the GCL-Mascon2024 recovery framework as an additional background force model.

2.3.2 Glacial isostatic adjustments

We apply GIA corrections in the GCL-Mascon2024 recovery process as another background force model. The official mascon products (i.e., CSR RL06 mascon, JPL RL06 mascon, and GSFC RL06 mascon) represent the surface mass deviation relative to the 2004.0-2009.999 time mean baseline. Subsequently, we model the GIA signals relative to the middle epoch of 2007.000, utilizing the GIA model ICE-6G, which was developed by Stuhne and Peltier (2015).

2.3.3 Earth ellipsoidal corrections

Temporal Stokes coefficients derived from GRACE satellite data are typically converted into mass anomalies at the Earth's surface using spherical harmonic synthesis, as formulated by Wahr et al. (1998). However, the results obtained using this approach reflect mass transport at a spherical surface with a fixed radius of 6378 km, which can introduce inaccuracies. Ditmar (2018) demonstrated that such a conversion may lack sufficient precision and proposed a revised formulation for converting Stokes coefficients into mass anomalies. This updated approach assumes that (i) mass transport occurs at the reference ellipsoid, and (ii) at each point of interest, the ellipsoidal surface is approximated by a sphere with a radius equal to the local radial distance from the Earth's center (the “locally spherical approximation”). In this study, we adopt the spherical harmonic synthesis method proposed by Ditmar (2018) to account for the effects of the Earth's oblateness and improve the accuracy of mass anomaly estimation.

2.3.4 Low-degree term corrections

Given the inherent limitations of the GRACE twin-satellite tracking, it is not feasible to determine the effects of geo-center motion, which can be represented in terms of time-varying degree-1 coefficients. Consequently, we utilize the coefficients derived by combining GRACE data with geophysical models (Sun et al., 2016). Furthermore, we incorporate the C₂₀ (degree-2 order-0) coefficients derived from satellite laser ranging (SLR) measurements (Chen et al., 2005; Cheng et al., 2013) to enhance accuracy. To this end and in line with previous studies (Watkins et al., 2015), the mascon grid solutions are first converted to the spherical harmonic coefficients by using spherical harmonic analysis. Then, we replace the low-degree terms (i.e., degree-1 and C₂₀) and utilize the spherical harmonic synthesis proposed by Ditmar (2018); the coefficients are converted back into mascon grid solutions to correct the implied low-degree term component of GCL-Mascon2024, considering the influence of the Earth's oblateness as detailed in Sect. 2.3.3.

2.3.5 GAD corrections

To explicitly contain seafloor pressure anomalies in the corrected mascon solutions, the AOD1B RL06 GAD product (Dobslaw et al., 2017) is reintegrated into the mascon calibration framework.

2.4 Frequency-dependent data weighting

The concept of FDDW originates from the fast collocation technique (Bottoni and Barzaghi, 1993), which assumes stationary measurement noise, thereby imparting a Toeplitz structure onto the noise covariance matrix. Subsequently, Ditmar et al. (2007) provided a detailed discussion of the FDDW concept and employed the technique to estimate the static Earth gravity field from the kinematic orbital acceleration of the CHAllenging Minisatellite Payload (CHAMP) satellite (Ditmar et al., 2006). The FDDW technique was later adapted for solving the temporal gravity field model using the GRACE inter-satellite acceleration (Liu et al., 2010). Afterward, Guo et al. (2018) utilized the FDDW technique to account for KBR (K-band ranging) frequency-dependent noise in the classical dynamic approach, leading to the development of the WHU RL01 model. Chen et al. (2019) further extended the application of the FDDW technique by incorporating both orbit and KBR frequency-dependent noise into the optimized short-arc approach and released the temporal gravity model named the Tongji-Grace2018 solution.

As indicated in numerous previous studies (e.g., Guo et al., 2018; Chen et al., 2019), the inter-satellite range rate measurements are affected by frequency-dependent noise. Before applying the FDDW technique, it is essential to build a stochastic noise model using, e.g., post-fit residuals from the GRACE measurements. As an example, we select the post-fit residuals from June 2009, calculated using the preliminary mascon solution for that month. As shown in Fig. 2a and c, the time series of post-fit residuals from orbit and range rate measurements on 5 June 2009 exhibit a clear dependence on frequency. This is further illustrated by the power spectral densities (PSDs) displayed in Fig. 2b and d, which indicate that both orbit and range rate measurements, particularly the former ones, are contaminated by low-frequency noise. The frequency-dependent noise in GRACE observations is largely attributed to errors in the GRACE orbits (Ditmar et al., 2012). This type of noise, in the essence of perfect orbital, instrumental, and/or other models, is typically addressed by either estimating (once or twice per orbital revolution) periodic parameters to account for unmodeled accelerations or by incorporating variance–covariance matrices to mitigate these errors (Zhou et al., 2024). In this study, noise-whitening filters W, constructed based on post-fit residuals derived from orbit and range rate measurements using the autoregressive (AR) noise model implemented in the ARMASA toolbox (Broersen and Wensink, 1998; Broersen, 2000), are applied to transform frequency noise ε into Gaussian white noise ε₀. Following the methodology of Chen et al. (2019), the variance–covariance matrix Σ can be constructed using the law of variance–covariance propagation:

$$\begin{array}{}\text{(5)}& \mathbf{\Sigma }={\mathit{W}}^{-\mathrm{1}}\cdot \mathrm{diag}\left[{\left(\mathit{W}\cdot \mathit{\epsilon }\right)}^{\mathrm{2}}/{\left({\mathit{\epsilon }}_{\mathrm{0}}\right)}^{\mathrm{2}}\right]\cdot {\left({\mathit{W}}^T\right)}^{-\mathrm{1}}.\end{array}$$

Figure 2 Time series and power spectrum densities (PSDs) of post-fit residuals from the orbit and KBR range rate.

2.5 Advanced spatial constraints

The linear system that connects satellite range rate observations to the mass anomalies within each mascon for estimation is rank-deficient. To stabilize the rank-deficient system of equations in mascon recovery, we employ Tikhonov regularization techniques (Tikhonov, 1963). Herein, we estimate the mascon elements using the following equation:

$$\begin{array}{}\text{(6)}& \widehat{\mathit{x}}={\left({\mathbf{A}}^T\mathbf{PA}+\mathit{\mu }{\mathbf{C}}_{\mathrm{M}}\right)}^{-\mathrm{1}}\cdot {\mathbf{A}}^T\mathbf{PL},\end{array}$$

where $\widehat{x}$ represents the estimated mascons without any corrections; A is the design matrix of partial derivatives; L is the residual vector, which is obtained by subtracting the kinematic orbit or KBR measurements from the reference orbit positions or KBR data; P is the weight matrix derived from the inverse of the variance–covariance matrix Σ (refer to Sect. 2.4); μ is the regularization factor; C_M is a diagonal constraint (or regularization) matrix of size n×n, named the mass variation regularization constraint normalized (MVRCN) matrix; and n is the number of the mascons to be estimated.

For the advanced spatial constraints, we construct the MVRCN matrix, which primarily comprises two components: one derived from the continental region aridity–wetness index, which is defined as the ratio of mean annual precipitation to mean annual reference evapotranspiration (Trabucco and Zomer, 2018), and the other from the ETOPO Global Relief Model of ice sheet regions (i.e., Greenland and Antarctica), an ice surface version that portrays the topography of the top layer of the polar ice sheets (MacFerrin et al., 2025). The fundamental premise is that humid basins on the continent require looser constraints for recovering higher temporal gravity signals, while arid basins require tighter constraints. Similarly, on polar ice sheets, areas at lower elevations necessitate looser constraints to recover mass variations, whereas regions at higher elevations require tighter constraints. Figure 3 shows the spatial distribution of the mascon-size MVRCN matrix. We employ the L-curve method to determine the appropriate regularization factor μ, employing monthly varying factor values to ensure that the resulting regularization matrix is sufficiently tight to suppress noise yet loose enough to allow the mascons to be adjusted to their optimal values.

Figure 3 The mass variation regularization constraint normalized (MVRCN) matrix used in the GCL-Mascon2024 recovery framework.

3 Short-arc approach for gravity field inversion

The short-arc approach, initially introduced by Schneider (1968), is a commonly utilized method for satellite gravity data inversion. Mayer-Gürr et al. (2005) introduced the short-arc approach to determine a CHAMP gravity field model. Mayer-Gürr (2008) further proposed a gradient correction algorithm to enhance the accuracy of the short-arc approach and applied it to real GRACE data inversion. Since then, the short-arc approach has been employed in processing GRACE data (e.g., Ran et al., 2014; Chen et al., 2019), demonstrating its effectiveness and efficiency in recovering temporal gravity field models. Section 3.1 is devoted to the optimal choice of the arc length. Next, Sect. 3.2 discusses the design of calibration parameter estimation during the gravity inversion process.

3.1 Arc length determination

Longer arcs (e.g., 24 h ones) are usually utilized in the dynamic approach to the temporal gravity solution recovery, whether it be in the form of the mascon solution (e.g., Watkins et al., 2015) or spherical harmonic solutions (e.g., Mayer-Gürr et al., 2018). Regarding the short-arc approach, the tendency is to select shorter arc lengths, such as 1 h arcs for Bonn University's ITG-GRACE2010 (Mayer-Guerr et al., 2010) and 6 h arcs for Tongji University's Tongji-Grace2018 (Chen et al., 2019). However, as the arc length decreases, the number of parameters per day increases. Given that the total number of observations remains constant, this increases the condition number of the estimation process in the temporal gravity field recovery. In the mathematical sense, the smaller the condition number of the normal matrix, the more stable the resulting estimate of the gravity field (Chen et al., 2019).

To determine the appropriate arc length for GCL-Mascon2024, we conducted computations of a monthly mascon model using different arc lengths to compare the stability of the resulting estimates. Table 2 presents the condition numbers of the unconstrained normal matrices and the corresponding computational time needed for different arc lengths. From this standpoint, the 2 h arc length corresponds to the most stable arc length in the GCL-Mascon2024 recovery. Figure 4 illustrates that increasing the arc length beyond 2 h in the short-arc approach leads to a significant increase in noise in gravity field estimates as the normal equations become more ill-conditioned. This observation aligns closely with what we conclude from Table 2. Therefore, an arc length of 2 h is determined to be the most suitable for the short-arc approach employed in this work. Additionally, we incorporate the gradient correction algorithm proposed by Mayer-Gürr (2008) to consider the influence of the kinematic orbit errors.

Figure 4 Geoid height differences per degree with regard to GOCO06s from mascon solutions of different arc lengths.

Table 2 Condition numbers (log10) of normal matrixes and inversion time cost in the GCL-Mascon2024 recovery framework with different arc lengths.

3.2 Calibration parameter estimation

The accelerometer represents a significant source of errors in the GRACE mission (Kim, 2000), necessitating the implementation of robust strategies to manage and mitigate accelerometer errors effectively. Simultaneously, in the analysis of GRACE observations, it is necessary to estimate not only the gravity field parameters but also arc-related parameters, such as the two boundary position vectors of each arc (Mayer-Gürr, 2008). That is, the error occurring at the boundaries of each arc is also of non-negligible magnitude. A commonly used strategy in temporal gravity field recovery is the incorporation of calibration parameters to mitigate the impact of the aforementioned errors.

The GRACE raw accelerometer measurements exhibit systematic errors, including bias, scale error, and drifts (Han et al., 2006b) along three axes (i.e., along, cross, and radial) for both satellites. The findings of Meyer et al. (2016) demonstrate that the scale calibration of accelerometer data at daily intervals significantly reduces the impact of solar activity on the derived gravity field models. To this end, we conduct the daily estimation of accelerometer scales along three axes for both satellites in this study. In addition, bias is a frequently employed parameter for estimating the local parameters of accelerometers (Kim, 2000). Based on prior studies (e.g., Han et al., 2006b; Bettadpur, 2007), we also incorporate the estimation of drift parameters into the recovery of the mascon solution. Combining the biases, drifts, and scales, the calibration formula for the accelerometer data can be constructed as

$$\begin{array}{}\text{(7)}& {\mathit{f}}_{\mathrm{new}}=\mathrm{bias}+\mathrm{scale}\times {\mathit{f}}_{\mathrm{ori}}+\mathrm{drift}\times t,\end{array}$$

where f_ori and f_new denote the nongravitational accelerations prior to and after calibration, respectively; bias, scale, and drift are the estimated local parameters of the accelerometers; and t represents the period during which the drift of nongravitational accelerations is calibrated. Figure 5 illustrates the geoid height differences per degree with respect to the GOCO06s of the mascon solutions, with accelerometer calibration parameters (i.e., bias, drift, and scale) co-estimated over different periods.

Table 3 Estimation periods of accelerometer calibration parameters (unit: minutes).

Table 3 provides a detailed definition of each considered case, characterized by three pre-defined periods for accelerometer calibration parameters: bias, drift, and scale. One can see from Fig. 5 that the inversion performs optimally when bias and drift are co-eliminated per arc, as well as with scale elimination on a per-day basis, with the premise of estimating the boundary position parameters per arc. After generating the normal equation for each arc, the calibration parameters of the boundary position can be eliminated immediately. Then, once the normal equations for a specific period are generated, the corresponding accelerometer calibration parameters are eliminated as well. Last, by combining all of the reduced daily normal equations, we obtain the final monthly normal equation, which is solved for the mascon coefficients.

As mentioned above, a 2 h arc is selected for the GCL-Mascon2024 computation. The calibration parameters for accelerometer observations include biases and drifts estimated per arc, as well as scales estimated per day, for the twin satellites along three axes.

Figure 5 Geoid height differences per degree with regard to GOCO06s from mascon solutions under different scenarios. Each scenario corresponds to a distinct set of parameters, reflecting variations in the estimation periods for accelerometer calibration parameters (i.e., bias, drift, and scale). Refer to Table 3 for detailed information on parameter settings.

4 Analysis of scientific results

To evaluate and validate the GCL-Mascon2024 solution, we compare the estimates of mass variation globally and over specific regions with the RL06 mascon solutions released by GSFC, CSR, and JPL. Here, annual amplitudes, monthly mass variations, basin hydrological signals, and polar region mass balances are utilized to assess the performance of temporal-signal retrieval. At the same time, continental random-noise and desert residuals are used to evaluate temporal noise.

4.1 Global comparisons

We first analyze the global mass change signals in GCL-Mascon2024 and in the RL06 mascon solutions provided by GSFC, CSR, and JPL. To emphasize the differences in the four mascon solutions, the results are presented as anomalies in relation to the baseline defined as the time mean during the period from January 2004 to December 2009. In Fig. 6, we specifically present the long-term trends in temporal gravity signals for the time span ranging from January 2003 to December 2015. Upon observing Fig. 6, we can discern a high level of consistency in the global mass change signals across all four models.

The annual amplitudes in mass change are depicted in Fig. 7 for the time span ranging from January 2003 to December 2015. It is evident that the spatial distribution of the four monthly mascon solutions exhibits a substantial level of concurrence. Regions characterized by a more pronounced annual fluctuation in total water storage are predominantly concentrated in specific areas, namely the Amazon basin in South America, the Niger basin in West Africa, the Zambezi basin in South Africa, and the Ganges and Mekong basins in Southeast Asia.

Figure 6 Long-term trends from January 2003 to December 2015 (in equivalent water height or EWH).

Figure 7 Annual amplitudes from January 2003 to December 2015 (in equivalent water height or EWH).

Ditmar (2022) proposed a technique to combine and regularize GRACE-based mass anomaly time series and, at the same time, to quantify the standard deviation (SD) of random noise in each time series. The latter is estimated using variance component estimation (VCE) as adapted from Koch and Kusche (2002). Figure 8 illustrates the spatial distribution of the random-noise SD estimated for various mascon solutions. The noise SD of the mass anomaly time series over the globe obtained for the mascon solutions from GSFC, CSR, JPL, and GCL-Mascon2024 are 3.5, 4.1, 3.9, and 3.7 cm, respectively. In northern Africa, the Arabian Peninsula, and eastern Asia (the border region between China and Mongolia), GCL-Mascon2024 and JPL mascon solutions exhibit similar spatial distributions, with smaller SDs in terms of random noise compared to GSFC and CSR solutions. Given the predominant desert coverage in these regions, it is reasonable to expect lower standard deviations in terms of random noise. Further quantitative analyses of random noise over specific local regions, including river basins, Greenland, and desert areas, are provided in the following section.

Figure 8 Spatial distribution characteristics of random noise of GSFC RL06 mascon, CSR RL06 mascon, JPL RL06 mascon, and GCL-Mascon2024 (in equivalent water height or EWH), with the standard deviation computed according to Ditmar (2022).

4.2 Regional comparisons

For a more comprehensive comparative analysis of signal magnitudes across various mascon solutions, this study selects distinct river basins, Greenland drainage systems, and typical deserts. These specific selections allow us to discern temporal signals associated with hydrological processes, ice melting dynamics, and temporal noise, respectively.

4.2.1 Hydrology

Continental water storage is a pivotal constituent within both terrestrial and global hydrological cycles, exerting a significant degree of control over intricate processes involving water, energy, and biogeochemical exchanges (Famiglietti, 2004). As such, it plays a paramount role in shaping and influencing the Earth's climate system (Chen et al., 2010). Of significant importance in terrestrial basins, the comprehensive analysis of total water storage (TWS) aids in understanding the intricate dynamics of water distribution and availability (Long et al., 2013). TWS refers to the summation of all water present within a given region, accounting for its various forms, such as surface water, groundwater, soil moisture, and snowpack. The GRACE mission can accurately capture the total mass variation caused by terrestrial water storage change (e.g., Ramillien et al., 2008; Rodell et al., 2018).

Given the potential divergence in the temporal signals of mass variations across river basins characterized by distinct sizes and climate classifications, we have statistically analyzed the temporal signals within the 42 largest basins (area > 5 × 10⁵ km²) in the world, which encompass different climate types. This selection intends to showcase the performance of the temporal-signal recovery by the different mascon solutions. The basic definitions of the aforementioned river basins are all taken from and credited to Scanlon et al. (2018).

Figure 9 illustrates the time series of basin mass variations derived from the WaterGAP Global Hydrology Model (WGHM); GCL-Mascon2024; and the mascon solutions from GSFC, CSR, and JPL. The WaterGAP model (Müller Schmied et al., 2021, 2024), primarily developed at the universities of Kassel and Frankfurt, simulates water flows, storage, withdrawals, and consumptive use globally, serving as a tool to evaluate the human–freshwater system under the influence of global change. As shown in Fig. 9, GCL-Mascon2024 exhibits a high level of agreement with the other models in terms of mass anomalies across all analyzed river basins. Using WGHM-based mass variations as control data, the time series derived from GCL-Mascon2024 for the 42 largest basins demonstrates a reduction of approximately 5.3 % in error compared to the other three mascon solutions released by GSFC, CSR, and JPL. Notably, in the Murray Basin, which exhibits the sub-arid climate type, the GCL-Mascon2024 time series shows a 44.6 %–58.0 % reduction in error compared to the other mascon solutions. As shown in Table 4, the correlation coefficients for mass variations within the selected regions between GCL-Mascon2024 and the other mascon solutions exceed 95.0 %.

Figure 9 Time series of mass anomalies over typical river basins from the hydrology model WaterGAP (outlined by the gray zone) and mascon solutions recovered by GSFC, CSR, JPL, and GCL (yellow, green, blue, and red lines, respectively). The base map illustrates the 42 largest basins (area > 5 × 10⁵ km²) extracted from the Total Runoff Integrating Pathway database, as from Scanlon et al. (2018).

Table 4 Correlation coefficients between mass anomaly time series over the representative river basins from GCL-Mascon2024 and official RL06 mascon solutions.

According to Table 5, the noise SD of the mass anomaly time series over the aforementioned river basins for the mascon solutions from GSFC, CSR, JPL, and GCL-Mascon2024 are 4.2, 4.6, 5.0, and 4.1 cm, respectively. It is important to highlight that the ability of the GCL-Mascon2024 solution to suppress random noise is optimal in most non-humid (i.e., subhumid, semiarid, and arid) basins. This indicates that the noise reduction of the GCL-Mascon2024 solution is 21.6 %, 29.2 %, and 32.6 %, respectively, compared to the GSFC, CSR, and JPL RL06 mascon solutions. Those improvements provided by the GCL-Mascon2024 solution may benefit from incorporating advanced spatial constraints derived from the aridity–wetness index of continental regions.

The results presented in Fig. 9 and Tables 4 and 5 demonstrate strong evidence that GCL-Mascon2024 is equally sensitive to hydrological signals compared to the official mascon solutions despite employing a shorter arc length (i.e., 2 h) and exhibiting a superior capacity for random-noise suppression.

Table 5 The root mean square of random noise over the 42 largest basins (area > 5 × 10⁵ km²) from the mascon solutions recovered by GSFC, CSR, JPL, and GCL. The definitions of these basin boundaries are derived from Scanlon et al. (2018). The bolded value indicates the lowest rms of random noise (unit: centimeters).

4.2.2 Cryosphere

The Greenland Ice Sheet (GrIS) is home to one of the largest freshwater reserves on our planet. Due to its substantial accumulation rate and considerable meltwater runoff, the GrIS is a highly dynamic system (Chen et al., 2006). Rapid transformations within the GrIS have the potential to raise the mean sea level substantially (Ran et al., 2024) and could significantly impact the North Atlantic thermocline circulation, thereby affecting the global climate (Velicogna and Wahr, 2005). One of the primary means for monitoring mass variation in the GrIS is the GRACE satellite mission (e.g., Schlegel et al., 2016; Velicogna et al., 2020).

Figure 10 Time series of detrended mass anomalies for individual drainage systems and the entirety of Greenland based on the mass balance from the input–output method, i.e., surface mass balance – ice discharge (outlined by the gray zone) and mascon solutions recovered by GSFC, CSR, JPL, and GCL (yellow, green, blue, and red lines, respectively). The middle panel presents the schematic illustration of the mascon division, and its base map portrays the topography of the Greenland Ice Sheet. In this study, Greenland is partitioned into 21 mascons and 7 individual drainage systems: north (NO), northeast (NE), northwest (NW), central east (CE), central west (CW), southeast (SE), and southwest (SW).

In Greenland, it is critical to emphasize that the mascon geometry of GCL-Mascon2024 is delineated based on the boundaries of the Greenland drainage system and the coastline. Greenland is partitioned into 21 mascons and 7 individual drainage systems: north (NO), northeast (NE), northwest (NW), central east (CE), central west (CW), southeast (SE), and southwest (SW). The various mascon solutions over different drainage systems of Greenland are validated using the input–output method (IOM) as control data, i.e., mass balance = surface mass balance −  ice discharge. Mass variations caused by surface mass balance (SMB) processes are derived from the MARv3.14.0 polar regional climate model run at a resolution of 10 km over the whole of the GrIS and with a 6 hourly forcing by the ERA5 reanalysis at its lateral boundaries and over the ocean (Fettweis et al., 2017). The middle panel of Fig. 10 presents the schematic illustration of the mascon division and the topography of the ice surface on Greenland. The other subfigures of Fig. 10 illustrate the time series of the detrended mass anomaly based on the mass balance from the IOM outlined by the gray zone and the different mascon solutions integrated over seven drainage systems, as well as over the entirety of Greenland. As indicated in Fig. 10, the time series of mass changes over Greenland is generally consistent across the four different mascon solutions, with all models effectively capturing the overall mass change in Greenland. The correlation coefficients of mass changes across the seven different drainage systems between GCL-Mascon2024 and the other three RL06 mascon solutions exceed 98.6 %. Furthermore, the correlation coefficient for capturing the total mass change of Greenland across all four models is as high as 99.8 %. Particularly in the northeastern drainage system of Greenland, where the mass variation is minimal, the time series for this region, extracted from GCL-Mascon2024, demonstrates a reduction in error of about 20 % compared to the other three mascon solutions from GSFC, CSR, and JPL. By extracting the noise SD of the mass anomaly time series within the Greenland drainage system from various mascon solutions (Table 6), we find that the noise SD for the GCL-Mascon2024 and GSFC RL06 mascon solutions is 8.7 and 9.7 cm, respectively, whereas it is 14.4 and 13.7 cm for the CSR and JPL RL06 mascon solutions. This indicates that the GCL-Mascon2024 solution achieves a random-noise reduction of 32.6 % and 29.2 % compared to the CSR and JPL RL06 mascon solutions. The observed discrepancies and the improvement offered by our mascon solution could be attributed to differences in the definition of mascon geometry and the processing methodology.

4.2.3 Desert

The preceding two sections have delved into the signal characteristics exhibited by the GCL-Mascon2024 solution over river basins and Greenland. In this section, we aim to evaluate the uncertainties of our mascon solutions over deserts and compare them with those of the other mascon solutions. Our impetus stems from an understanding that precipitation within desert regions is limited. It is critical to emphasize that aridity cannot be equated with negligible temporal mass variations (e.g., Scanlon et al., 2022). Conversely, low precipitation may stimulate an extensive consumption of groundwater. To that end, the residuals, calculated after removing the climatological components (i.e., bias, trend, and amplitude) from the mass variations, can be regarded as mis-modeling signals or temporal noise signals that persist in the temporal gravity fields (e.g., Zhou et al., 2024). Consequently, we analyze the error characteristics inherent to the mascon models over typical deserts, such as the Sahara Desert in Africa, the Taklamakan Desert in Asia, and the Atacama Desert in South America.

Deserts are territories characterized by low precipitation. They can be classified into several categories based on their respective geographical locations and prevailing weather patterns, which include trade wind deserts, rain shadow deserts, and coastal deserts (Whitford and Duval, 2019). Trade wind deserts are typically found on both sides of the horse latitudes, between ±30° and ±35°. These regions are characterized by subtropical anticyclones and the large-scale descent of dry air masses (Glennie, 1987). The Sahara Desert, the largest hot desert in the world, is an example of this type. By extracting the mass variation residuals of the Sahara Desert from varying mascon solutions, the residual of the GCL-Mascon2024 solution and the JPL RL06 mascon solution is 62.7 and 50.5 Gt, but it is 84.5 and 70.2 Gt for the GSFC and CSR RL06 mascon solutions, respectively. This indicates that the noise reduction of the GCL-Mascon2024 solution is 25.8 % and 10.7 %, respectively, when compared to the GSFC and CSR RL06 mascon solutions. Rain shadow deserts are formed by the rain shadow effect. Orographic lift forces air masses to rise over mountains, cooling and losing moisture on the windward slopes. As the air descends on the leeward side, it warms, increasing its moisture capacity and creating a drier region with reduced precipitation (Sun et al., 2008). The Taklamakan Desert, the largest in China, located in the rain shadow of the Himalayas, exemplifies this phenomenon. The residuals of mass variations in this region are estimated to be 3.3, 2.0, 1.2, and 1.7 Gt according to the GSFC RL06, CSR RL06, JPL RL06, and GCL-Mascon2024 mascon solutions, respectively. The Atacama Desert, a prime example of a coastal desert, is one of the driest regions on Earth, characterized by an almost complete absence of life due to its extreme aridity. This hyperarid climate is primarily caused by the orographic effects of the Andes Mountains to the east and the Chilean Coast Range to the west, which prevent the desert from receiving significant precipitation. Additionally, the cold Humboldt Current and the persistent Pacific anticyclone play critical roles in maintaining the region's dryness (Westbeld et al., 2009). The root mean square (rms) of mass variations over the Atacama Desert, as derived from the mascon solutions by GSFC, CSR, JPL, and GCL-Mascon2024 are 3.8, 1.9, 1.4, and 1.2 Gt, respectively, indicating that the GCL-Mascon2024 solution has the smallest error.

Figure 11 Time series of mass change residuals over deserts derived from the RL06 mascon solutions from GSFC, CSR, JPL, and GCL-Mascon2024. The residuals indicate that the climatological components (i.e., bias, trend, and amplitude) have been removed from the mass variation. The deserts chosen are the Sahara Desert, Sechura Desert, Atacama Desert, Kyzylkum Desert, Gobi Desert, Taklamakan Desert, and Arabian Desert.

Table 6 The root mean square of random noise over individual drainage systems of Greenland and desert regions from the mascon solutions recovered by GSFC, CSR, JPL, and GCL. The bolded value indicates the lowest rms of random noise (unit: centimeters).

Figure 11 illustrates the mass variations and the rms of residuals of typical deserts. The deserts selected for this study include the Sahara Desert, the Sechura Desert, the Atacama Desert, the Kyzylkum Desert, the Gobi Desert, the Taklamakan Desert, and the Arabian Desert. The GCL-Mascon2024 incorporates well-defined physical constraints, such as coastal and basin boundaries, along with advanced spatial constraints based on the MVRCN matrix, enabling it to reduce errors in desert regions by approximately 29.3 % compared to the GSFC and CSR RL06 mascon solutions. Meanwhile, the JPL RL06 mascon demonstrates a slightly superior error suppression capability compared to the GCL-Mascon2024 solution in the aforementioned deserts. Notably, especially in the Atacama Desert, which is a long and narrow coastal desert from north to south and the driest desert in the world, GCL-Mascon2024 can achieve noise suppression ranging from 35.5 % to 68.2 % compared to the mascon solutions provided by GSFC and CSR. As shown in Table 6, the noise SD of the mass anomaly time series over the selected desert regions for the GSFC, CSR, JPL, and GCL-Mascon2024 mascon solutions are 2.1, 1.5, 1.5, and 0.9 cm, respectively. This translates to a random-noise reduction ranging from 40.0 % to 57.1 % compared to the GSFC, CSR, and JPL RL06 mascon solutions.

4.2.4 Lake and ocean

The utilization of mass variations in large lakes (e.g., the Caspian Sea) to assess noise levels in GRACE solutions is a well-established approach (e.g., Loomis and Luthcke, 2017; Ditmar, 2022). Herein, we choose the largest lake on Earth, the Caspian Sea, as an example for verification. We follow the approach proposed by Ditmar (2022), wherein the mass anomaly time series derived from GRACE is compared with the water level time series obtained from satellite altimetry observations. The latter time series is empirically rescaled (with a scaling factor of 0.687 for the Caspian Sea, as provided by Ditmar, 2022) to account for signal damping in the GRACE solution. Figure 12 presents the mass anomaly time series over the Caspian Sea derived from various mascon solutions and satellite altimetry data. As illustrated, the GCL-Mascon2024 solution shows strong consistency with the other models in capturing mass variations in this region. Using satellite-altimetry-derived mass variations, scaled by a factor of 0.687, as the reference, the noise SD values for the GSFC, CSR, JPL, and GCL-Mascon2024 mascon solutions are 5.7, 5.8, 5.6, and 5.2 cm, respectively.

Figure 12 Comparison of GRACE-derived mass anomaly time series (expressed in equivalent water height, EWH) from different mascon solutions with satellite-altimetry-based water level variations over the Caspian Sea. The time series derived from satellite altimetry has been downscaled using a scale factor of 0.687 to account for signal attenuation (Ditmar, 2022).

GRACE satellite gravity measurements over oceanic regions directly correspond to ocean bottom pressure variations at spatial scales of ∼ 300 km (Watkins et al., 2015). Figure 13 illustrates the time series of basin mass variations derived from different mascon solutions. To assess the quality of our solutions for ocean signals, we compute the correlation coefficients between GCL-Mascon2024 and the RL06 mascon solutions released by GSFC, CSR, and JPL. The resulting correlations are 95.7 %, 98.0 %, and 98.2 %, respectively, indicating a high level of consistency between our products and official mascon products.

Figure 13 Comparison of GRACE-derived mass anomaly time series (expressed in equivalent water height, EWH) over the global sea from different mascon solutions.

5 Data availability

All datasets used in this study were last accessed on 20 August 2025. The specific data repositories include GRACE level-1B data, downloaded from JPL at https://podaac.jpl.nasa.gov (Case et al., 2010), and kinematic orbits available from the Graz University of Technology at ftp://ftp.tugraz.at (Suesser-Rechberger et al., 2022). The ITSG-Grace2018 monthly solutions can be accessed via https://doi.org/10.5880/ICGEM.2018.003 (Mayer-Gürr et al., 2018). The RL06 mascon solutions released by JPL, CSR, and GSFC are available at, respectively, https://doi.org/10.5067/TEMSC-3JC63 (Wiese et al., 2023), https://doi.org/10.15781/cgq9-nh24 (Save, 2020), and https://earth.gsfc.nasa.gov/geo/data/grace-mascons (Loomis et al., 2019). The WaterGAP Global Hydrology Model, for comparisons, can be downloaded from https://doi.org/10.1594/PANGAEA.948461 (Müller Schmied et al., 2022). The MAR (version 3.14) model used in this study can be downloaded from http://ftp.climato.be:80/fettweis/MARv3.14/Greenland/ (Fettweis et al., 2017). The time series of water level variations in the Caspian Sea is derived from satellite altimetry data provided by the USDA/NASA G-REALM program (https://ipad.fas.usda.gov/cropexplorer/global_reservoir/, Birkett, 1995). The GCL-Mascon2024 model is available at https://doi.org/10.5281/zenodo.15525467 (Yan and Ran, 2025).

6 Conclusions

Mascon solutions of Earth's temporal gravity field can be considered to be more “user-friendly” compared to spherical harmonic solutions as they remove the need to apply empirical post-processing filters to eliminate errors in the unconstrained spherical harmonic solutions. Given this significant advantage, mascon solutions have been garnering increased interest from the GRACE applications community. Herein, the Geodesy and Cryosphere Laboratory from the Southern University of Science and Technology presents a novel mascon solution named GCL-Mascon2024, derived utilizing the short-arc approach and the level-1B data from GRACE. GCL-Mascon2024 features uniquely variably shaped mascon geometries integrated with relevant physical constraints such as coastline and basin boundary geometry, which ensures an accurate representation of temporal gravity signals while minimizing signal leakage. Meanwhile, this series of mascon recovery processes incorporates frequency-dependent data-weighting techniques to reduce the influence of low-frequency noise in observations. GCL-Mascon2024 utilizes advanced spatial constraints based on the MVRCN matrix, which is constructed by integrating a priori basin climate factors and cryosphere elevation models. The MVRCN matrix is carefully incorporated into the inversion process as a regularization matrix to minimize errors, ensuring the improvement of the signal-to-noise ratio in the GCL-Mascon2024 recovery framework.

To evaluate the quality of GCL-Mascon2024, we analyze the signal and error levels across continental regions globally, assess signal strengths over selected river basins and Greenland, and examine noise levels in representative desert areas. Based on these analyses, we draw the following conclusions:

1. Over the continental regions, GCL-Mascon2024 mass anomaly estimates from GRACE data show strong agreement with the RL06 mascon solutions (GSFC, CSR, JPL) in both spatial and temporal domains. For global ocean signals, the correlation between GCL-Mascon2024 and the RL06 mascon products from GSFC, CSR, and JPL exceeds 95.7 %.

2. The long-term trends and amplitudes from GCL-Mascon2024 over river basins and Greenland exhibit strong consistency with the RL06 mascon solutions from GSFC, CSR, and JPL. In particular, within non-humid river basins, the GCL-Mascon2024 suppresses random noise in the range of 21.6 % to 32.6 % compared to contemporary mascon products. With SMB-based mass balance as the benchmark, GCL-Mascon2024 achieves about 20 % error reduction compared to the other three mascon solutions in the northeastern drainage system of Greenland, where mass variation is minimal.

3. Mass variations in deserts, regions characterized by low precipitation, are typically minimal, offering an ideal basis for assessing the temporal errors of different mascon models. Building on this premise, the work investigates the error characteristics across diverse desert types, including the Sahara Desert (trade wind type), the Taklamakan Desert (rain shadow type), and the Atacama Desert (coastal type), along with other deserts. The GCL-Mascon2024 reduces temporal errors in these desert regions by approximately 29.3 % compared to the RL06 mascon solutions from GSFC and CSR. Meanwhile, GCL-Mascon2024 achieves a random-noise suppression ranging from 40.0 % to 57.1 % compared to the other three mascon solutions.