Introducing an Improved GRACE Global Point-Mass Solution—A Case Study in Antarctica

: In the so-called point-mass modeling, surface densities are represented by point masses, providing only an approximated solution of the surface integral for the gravitational potential. Here, we propose a reﬁnement for the point-mass modeling based on Taylor series expansion in which the zeroth-order approximation is equivalent to the point-mass solution. Simulations show that adding higher-order terms neglected in the point-mass modeling reduces the error of inverted mass changes of up to 90% on global and Antarctica scales. The method provides an alternative to the processing of the Level-2 data from the Gravity Recovery and Climate Experiment (GRACE) mission. While the evaluation of the surface densities based on improved point-mass modeling using ITSG-Grace2018 Level-2 data as observations reveals noise level of approximately 5.77 mm, this ﬁgure is 5.02, 6.05, and 5.81 mm for Center for Space Research (CSR), Goddard Space Flight Center (GSFC), and Jet Propulsion Laboratory (JPL) mascon solutions, respectively. Statistical tests demonstrate that the four solutions are not signiﬁcant different (95% conﬁdence) over Antarctica Ice Sheet (AIS), despite the slight differences seen in the noises. Therefore, the estimated noise level for the four solutions indicates the quality of GRACE mass changes over AIS. Overall, AIS shows a mass loss of − 7.58 mm/year during 2003–2015 based on the improved point-mass solution, which agrees with the values derived from mascon solutions.


Introduction
The observations acquired by the time-resolved satellite gravimetry namely Gravity Recovery and Climate Experiment (GRACE) have provided essential information on the global mass changes since April 2002 [1]. For instance, the monthly fields of mass changes inverted from GRACE observations have been used in monitoring ice sheet and glacier mass balance as well as for understanding their contributions to sea level changes (e.g., Reference [2]). GRACE monthly solutions of the time-varying Earth's gravity field are usually expressed in terms of spherical harmonic coefficients called Level-2 data. The Level-2 products are freely distributed by the centers identified in the GRACE Science Data System (SDS). These centers are the University of Texas Center for Space Research (CSR, [3]), the NASA Jet Propulsion Laboratory (JPL, [4]), and the German Research Centre for Geosciences (GFZ, [5]). Likewise, other institutions have provided independent solutions for Level-2 data based on alternative approaches (see, e.g., References [6][7][8][9][10][11][12]). Moreover, the recently established International Combination Service for Time-variable Gravity Fields (COST-G) generates combined monthly Level-2 data using

The Improved Point-Mass (IPM) Solution
Considering that the static gravity field of the Earth has been removed, the time variable residual gravitational potential as sensed by the satellite gravimetry can be seen as a column from the satellite's altitude to the center of the Earth due to large-scale mass variations. If hydrological mass variations are in the focus of interest, the other contributions from atmosphere and ocean mass variations, polar motion, solid and oceanic tides, and solid Earth variations (e.g., due to glacial isostatic adjustment-GIA-and mass loadings), can be reduced using background models and assumptions of a deformable Earth, respectively. Thus, Newton's integral for the residual gravitational potential (δV) due to surface densities, which are assumed to be concentrated on a virtual layer of water thickness ∆r (equivalent water height-hereafter abbreviated as e.w.h.) on Earth's surface, is given as: (r + 3r cos ψ) + r 2 (3cos 2 ψ − 1) which was discretized for each panel k. In Equation (1), G is the Newtonian constant of gravitation (G = 6.67408 × 10 −11 m 3 kg −1 s −2 , cf. Reference [41]), w is the density of fresh water (position-dependent, which can be considered constant ∼ 1, 000 kg/m 3 ). The Euclidean distance ( ) between the (attracted) computation point P(r, ϕ, λ) and the running integration (attracting) point Q(r , ϕ , λ ) is given by the al-Kāshi theorem: with ψ being the spherical distance as: cos ψ = sin ϕ sin ϕ + cos ϕ cos ϕ cos λ − λ .
Furthermore, ∆r = r 2 − r 1 is the e.w.h. in meters, which is assumed to generate the residual gravitational potential δV at the satellite gravimetry. The integral in Equation (1) is an elliptic integral in terms of the geocentric spherical coordinates (ϕ , λ ) and, as such, it is not integrable in elementary terms. Taking the first-order approximation of the kernel in Equation (1) into account, that is, neglecting O(∆r 2 ), which gives 2∆rR 2 / (see Equation (A10) in Appendix A), the resulting integral becomes an elliptic integral of the first kind. In Equation (4), is the Euclidean distance given by (A9), and the Landau symbol O(∆r 2 ) indicates that the second-order terms in ∆r are omitted. Given the first-order approximation of the kernel appearing in the surface integral in Equation (4), and considering the surface density µ given by w · ∆r, the remain integral becomes equivalent to the potential of a single layer (or density layer, simple layer). The so-called mascon solutions represent the time-varying gravity field as the potential of a single layer [21]. Such solutions assume that the mass changes take place at the Earth's surface [42]. This assumption is somewhat valid since the time-varying signals mainly reflect the superficial processes such as the water cycle given that the non-surficial sources such as GIA have been removed. Consequently, the challenges due to the non-uniqueness of the inversion of the volume integral to find a unique distribution of masses enclosed by a body given the external potential as observations are overcome when considering the surface integral for the Earth's surface, of which the inversion is unique [42].
Since elliptic integrals occur in Equation (4), it cannot be solved by elementary integration. Therefore, it is proposed to perform a Taylor series expansion of the integrand in Equation (4), that is, the kernel function, and its subsequent integration. Noteworthy, this has already been suggested for the volume integral by, for example, Heck & Seitz [40], Wild-Pfeiffer [43], and Grombein et al. [44]. The solution presented by these authors has been proposed in the context of the forward modeling of the gravity field due to the topographic masses, and it would become a nonlinear problem in the present context. Nevertheless, the kernel of the integral given in Equation (4) (that is, second-order and higher terms in ∆r were not considered) can be formulated in Taylor series as: where the coefficients K i,j are the partial derivatives of the integral kernel K defined by The coefficients K i,j in Equation (6) depend on the relative positions of the computation point P(r, ϕ, λ) at satellite altitude and the geometric center Q 0 (R, ϕ 0 , λ 0 ) of the spherical panel (i.e., the Taylor point defined by ϕ 0 = (ϕ 1 + ϕ 2 )/2 and λ 0 = (λ 1 + λ 2 )/2 ). In this case, the zeroth-order coefficient K 0,0 results in where 0 is given by Equation (A9) at which ψ = ψ 0 (center of the panel). The second-order coefficients K 2,0 and K 0,2 , computed by means of Equation (6), are given as and respectively. Note that in the context presented here, where R < r, no singularities can occur in Equations (7)-(9) when ψ = 0 • .
Considering a Taylor series expansion up to degree 3, the gravitational potential according to the discretization presented in Equation (4) can be expressed by where ∆ϕ = ϕ 2 − ϕ 1 and ∆λ = λ 2 − λ 1 denote the horizontal dimensions of a spherical panel k.
In Equation (10), O(∆ 4 ) indicates that the fourth-order terms in ∆ϕ and ∆λ are omitted. That is, Equation (10) is approximate (accurate) up to third-order terms. Due to the specific choice of the Taylor point Q 0 , only the terms with even powers of i and j in the series expansion remain after integration. In the case of ϕ , the same holds for λ , the remaining even power can be verified by inserting Equation (5) into Equation (4), which provides: where Γ is the gamma function, which for an integer i the function relates with factorial as Γ(i + 1) = i!.
Hence, the odd-order terms in i or j in the Taylor series in Equation (5) cancel out after subsequent integration as per Equation (4) and only terms with even-order remain. The zeroth-order approximation of Equation (10), that is, neglecting first-order terms and higher, provides the so-called point-mass solution (hereafter referred to as "original point-mass" solution and abbreviated as "OPM"). Therefore, considering higher-order terms in Equation (10) represents the refinement of the original point-mass solution (henceforward referred to as "improved point-mass" solution and abbreviated as "IPM"), which was specified as objective (ii) in Section 1. Thus, the OPM solution only provides an approximate solution of the surface integral appearing in Equation (4), where the linear approximation of its integrand, which provides 1/ , is considered to be constant over each panel. Furthermore, it is the functional model for the inverse modeling, which relates the observations (residual gravitational potential) and the parameters (surface densities). It considers terms up to third-order and is the main contribution of the present work with regard to the solution used in previous studies (cf. References [34][35][36][37][38][39]). Equation (10) has an advantage w.r.t. the formulation presented by Heck & Seitz [40] since there are no nonlinear terms in the mathematical model of type G(m) = d. For comparison with a 3-D tesseroid approach, please refer to the contribution of Grombein et al. [45].

The Surface Integral Considering the Exact Kernel
In Section 2.1.1, it was presented the zeroth-order (i.e., OPM solution) and third-order (i.e., IPM solution) approximations of Equation (1), and both require assessment. Therefore, it is necessary to find a solution for the gravitational potential that is superior and independent to both, IPM O(∆ 4 ) and OPM O(∆ 1 ) solutions. The kernel in the integral in Equation (1) can be rewritten as I(ϕ , λ , r ) = cos ϕ (r + 3r cos ψ) + r 2 (3cos 2 ψ − 1) ln + r − r cos ψ r =R+∆r which, when evaluated at the Taylor point Q 0 , gives: where 0 is the same as defined for Equation (1) at which ψ = ψ 0 . The partial derivatives are: with and I 0,2 =r cos ϕ 0 6r(K 2 λ + K λλ cos ϕ 0 ) ln(µ 0 ) with Thus, considering the above developments, the potential in Equation (1) can be written as where the terms with odd orders and the mixed partial derivatives vanished after integration as shown in Equation (11). Within the solutions provided by Equation (20), singularities occur when ψ = 0 and r ≥ r . For satellite geodesy applications, as in the current case, r > r . Therefore, when ψ = 0, the argument of the natural logarithm (µ 0 ) is zero. Equation (20) is superior to Equation (10) since the radial integration was analytically performed, see the kernels of the integrals in Equation (1) and Equation (4), which differ by O(∆r 2 ). Equation (20) is used to evaluate the inverse modeling employing the OPM and IPM solutions as the functional models for recovering the surface densities at the Earth's surface.

GRACE Level-2 Data
The ITSG-Grace2018 monthly gravity field solutions generated by the Institute of Geodesy of the Graz University of Technology (see References [8,46] for further details) were used up to degree and order (d/o) 120 to synthesize the gravitational potential (residual). The ITSG-Grace2018 data is available at www.tugraz.at/institutes/ifg/home/. The degree-1 coefficients (C 1,0 , C 1,1 , and S 1,1 ), which represent the changes in the geocenter, and the degree-2 coefficient (C 2,0 ), which is associated with the oblate shape of the geopotential, were replaced with those computed from the methodologies of Swenson et al. [47] and Cheng & Tapley [48], respectively. They are also based on Release 06 (RL06), which uses new background models and were retrieved from https://grace.jpl.nasa.gov/. The GIA contribution to the secular changes of GRACE inferred surface densities was removed using the model provided by A et al. [49], which is available at https://grace.jpl.nasa.gov/. However, signals due to large earthquakes have not been considered; hence, the analysis of water mass gain/loss must be treated with caution in regions near large earthquakes.

Experimental Design
At this stage, it is necessary to investigate the OPM solution, as wished-for in objective (i), and to validate the IPM solution, as per objective (ii). Furthermore, an application of GRACE Level-2 data to assess the mass changes over Antarctica, as proposed in objective (iii), is carried out. Firstly, an experiment is advised to validate both original and improved point-mass solutions in closed-loop simulation. Secondly, as to apply the proposed improved method in the context of GRACE Level-2 data, it is performed the inversion of global residuals of the gravitational potential into mass changes expressed as e.w.h. Figure 1 presents the main steps of the experimental procedure.
The first experiment is based on inverse modeling in a closed-loop scheme ( Figure 1). However, in this experiment, it is necessary to know the residual gravitational potential, which has to be superior and independent of Equation (10). Hence, the monthly residuals of the gravitational potential are computed employing Equation (20) at 500 km altitude (approximately the GRACE's altitude) using the respectively monthly fields of e.w.h. from the GSFC-M (Section 2.2.2). Next, inverse modeling is conducted using the forwarded gravitational potentials as the observations and the zeroth-and third-order approximations of Equation (10), that is, the OMP and IPM solutions, respectively, as the functional models. The parameters are a set of equal-area panels covering the whole Earth. The choice of spatial sampling using panels with a constant area equivalent to the area of 1-by-1 arc-degree quadrangle at the equator (approximately 12,390 km 2 ) is to characterize better the units (e.g., basin and regions) and the coastlines. The partitioning of the sphere into equal-area panels used here is precisely the one used for GSFC-M, as shown in Figure 2, in which further details can be found in Luthcke et al. [24]. Albeit the inversions using GSFC-M as synthetic data are carried out at a global scale, the analyses of the different inversion results are also considered for Antarctica Ice Sheet (AIS). The outcomes from these experiments are presented in Section 3.1.
GSFC-M monhly fields of equivalent water height (e.w.h.) ITSG-Grace2018 ("take data as they are") synthesized potentials at 500 km altitude, which are taken as the observations inverted e.w.h. fields using Eq. (9) in its third-order approx. as functional model inverted e.w.h. fields using Eq. (9) in its zeroth-and third-order approximations forwarded potentials at 500 km altitude using Eq. (19) error assessments of the inversions of Eq. (9) in its zeroth-and third-order approx.
do higher order terms in Eq. (9) provide any improvement w.r.t.
the zeroth-order approx.?   (10) as well as a study case over Antarctica using GRACE Level-2 data as given by ITSG-Grace2018. Finally, a case study over the AIS is carried out using global monthly synthesized residual gravitational potentials at 500 km altitude. The residual gravitational potentials are computed from the ITSG-Grace2018 Level-2 data up to d/o 120, and they are then used as observations in the inversion, which the third-order approximation of Equation (10) is taken as the functional model. The mass change time series are averaged over AIS considering its limits given by the grounding line ( Figure 2). The AIS' delimitations are those according to Rignot et al. [53] as provided by NASA National Snow and Ice Data Center [54]. The outcomes from these experiments are presented in Section 3.2.

Assessment of OPM and IPM by Means of a Closed-Loop Simulation
One of the goals of this study is to validate the original point-mass (OPM) modeling since it has been identified as a viable methodology for processing GRACE Level-2 data (cf. Reference [36]) and applied in many relevant studies (e.g., References [34][35][36][37]39]). In this context, a closed-loop simulation is performed to assess the point-mass solution. In this section it is presented the evaluation of the zeroth-and third-order approximations of Equation (10)   The monthly e.w.h. fields were used with Equation (20) to forward the respective gravitational potentials at an altitude of 500 km with a spatial resolution of 1 • in ∆ϕ and varying ∆λ per latitude-band, totalizing 41,168 points. Figure 3b summarizes the respective gravitational potentials in terms of the RMS values. The use of the integral in Equation (20) to forward the gravitational potential at GRACE's altitude is assumed to be an independent and superior solution for the gravitational potential regarding Equation (10). This forward modeling, which aims to find the gravitational potential at 500 km altitude given the water masses on the Earth's surface, is essentially an upward continuation. The upward continuation works as a low-pass filter, as can be seen in Figure 3b, where the gravity field weakens with altitude, and short wavelengths are more attenuated than longer ones. Similar to Figure 3a, the spatial patterns of the gravitational potential shown in Figure 3b illustrate the long-wavelength features of the water masses, where the changes over South America and the tropical wet and dry regions of Africa remain evident. Noteworthy, the values of the gravitational potential over Greenland and West AIS are shown higher than central South America, specifically at Amazon basin. These values over West AIS are due to the ice melting rather than the seasonality due to the accumulation and oblation variations over AIS. (The simulated mass changes represented by the GSFC-M were not de-trended.) Inverse modeling was performed to estimate the surface densities as the parameters (m), expressed as e.w.h. ∆r, on the spherically approximated Earth's surface with radius 6378 km with the gravitational potential fields at 500 km altitude taken as observations. This was performed in two different ways. First, the gravitational potential values at 41,168 bins with altitude equals 500 km were computed through Equation (20) due to the 11,930 mass change panels from GSFC-M covering only the continents (i.e., the oceans were masked out). Second, the gravitational potential values at 41,168 bins at an altitude of 500 km were also computed through Equation (20) due to the 41,168 mass change panels from GSFC-M covering whole Earth. Noteworthy, the use of GSFC-M mass changes and Equation (20) is to ensure independent observations (d) in a closed-loop scenario to evaluate Equation (10) in its zeroth-(OPM) and third-order (IPM) approximations. Equation (10) become a linear matrix equation as d = Gm, where the coefficient matrix G ∈ R n×m , has a dimension of 41,168-by-11,930 considering only continents, and 41,168-by-41,168 given the entire Earth. Such a system is rank-deficient, which produces an infinite number of solutions. Consequently, the inverse problem was solved by using the Tikhonov-regularized least-squares at which the regularization parameter for each inversion was found using the L-curve criterion (see Reference [55] for details). The inversions through January 2003 to December 2015 (142 months) were performed using the respective zeroth-and third-order approximations of Equation (10), that is, OPM and IPM formulations, as the functional models. The results are summarized in Figure 4 in terms of RMS of the errors (RMSE) regarding the initial mass changes inputs.
In general, the OPM solution (zeroth-order approximation) presents larger RMSE values of the inverted masses regarding the IPM (third-order approximation), as summarized in Table 1. The IPM solution presented an area-weighted averaged RMSE of 2.90 mm against 23.82 mm for the original point-mass (reduction of approximately 87%) considering the recovered masses over the continents (i.e., the oceans were masked out in the forward and inverse modeling). For the OPM solution, large RMSEs concentrate at the coastlines, as can be seen, for example, at the Antarctica region ( Figure 4a). Such patterns are not apparent in the RMSE values of surface densities recovered by the IPM modeling ( Figure 4b). Furthermore, the evaluation only over the Antarctica region presents an RMSE value of 58.94 mm and 9.46 mm for OPM and IPM solutions, respectively (Table 1). Despite the high RMSE values over Antarctica, yet the OPM presents a good performance based on the ratio of the RMSE to the standard deviation of the observed data (RSR, cf. Reference [56]). (RSR values less than 0.60 are seen as indicators of good performance as discussed in Reference [56]).
However, edge effects are apparent in the inversion targeting only the continents as seen in Figure 4a (the same holds at the global scale, which has not been shown here). In order to reduce the edge effects, it was also performed an investigation by considering the parameters covering the whole Earth, which gives a total of 41,168 parameters in comparison with the initial 11,930 parameters over the continents. The edge effects shown in Figure 4a were substantially reduced at which the occurrences of the patterns are less and with lower amplitudes, as shown in Figure 4c in the case of the OPM solution. The RMSE values changed from 23.82 mm to 2.41 mm, considering the evaluation only over the landmasses (Table 1). Nevertheless, the IPM modeling still performs better than the OPM formulation in recovering the signals with almost no noise (0.16 mm), as shown in Figure 4d. Furthermore, the weighted RMSE values considering the whole Earth are 1.68 mm and 0.13 mm for OPM and IPM solutions, respectively. The averaged (weighted) RMSE values over Antarctica presented very good performance as depicted by the RSR in Table 1

Inversion of GRACE Level-2 Data Using the IPM Modeling
The feasibility of the IPM formulation based on a closed-loop simulation showed that it delivered better results than the OPM modeling (Section 3.2). Thereby, in this section, only IPM is considered to provide a practical application as a post-processing alternative for the GRACE Level-2 product (see Figure 1). The monthly ITSG-Grace2018 solutions (Section 2.2.1) were used to synthesize the respective monthly residual gravitational potential fields from Jan 2003 to Dec 2015. Monthly regional grids of gravitational potential were synthesized at 500 km altitude with spatial resolutions of 1 • -by-1 • given 64,800 bins as monthly observations. Each set of monthly fields formed the observations required in the inversion of Equation (10), in which the parameters (mass changes) were estimated. The mass changes were parametrized by the surface density values of 41,168 distributed, equal-area spherical panels, which were the unknown parameters (m). Figure 5a shows the resulting summary of the spatial distribution of the linear trends of e.w.h. series for each equal-area panel over whole AIS. The linear trends refer to the best fit of the series considering constant, linear trend, annual, semi-annual, and 161-days components regarding the period from January 2003 to December 2015. The spatial patterns of the estimated linear trends based on the IPM and those from mascon solutions generally agree within the AIS (    A multiple-comparison test [57] was considered to determining the similarities (if any) of the monthly series presented in Figure 6a. A nonparametric Kruskal-Wallis test [58] at a 95% confidence interval was applied to the monthly series of the four solutions to determine whether there were statistically significant differences among them. Figure 6b shows the box plots for the respective solutions at which their notches overlap, indicating that the actual medians of the mass changes series do not differ at 95% confidence. Subsequently, a follow-up test was conducted to identify which solution presented monthly values of mass changes resulting from a different distribution. The relative performance of the four solutions using the multiple-comparison test, based on the Tukey-Kramer procedure, indicated that there is no significant difference among the series. While the series of the four solutions are not significant different (Figure 6b), this may imply that they have the same level of accuracy over the AIS. Nevertheless, validation of the GRACE solutions in terms of e.w.h. requires independent datasets, which are not available. However, attempts in evaluating the different GRACE solutions have been proposed based on selecting a region with minimal mass variations, like Sahara desert [7], the analysis of the de-trended and de-seasoned residuals [59], the ensemble prediction and inter-comparison analysis [60], and the noise decoupling problem [61]. The assessment suggested by Groh et al. [59] was carried out to evaluate the four time series shown in Figure 6a. Specifically, the residuals of the de-trended and de-seasoned series (already discussed) were smoothed using a 13-month moving average, and the differences between both (original and filtered residuals) as shown in Figure 6c were used to assess the noise level of the solutions. The high-pass filtered residuals shown in Figure 6c present low correlation (autocorrelation). This indicates that more noise dominates the high-pass filtered residuals of the respective solutions. Figure 6d shows the box plots summarizing the high-pass filtered residuals at which the RMSEs are 5.02, 6.05, 5.81, and 5.77 mm, for CSR-M, GSFC-M, JPL-M, and IPM, respectively. That is, IPM ranks second best among the four solutions.

Discussion
Two of the main goals of this experiment were (i) investigating the point-mass modeling, and (ii) proposing an improved alternative solution. Thus, a refinement for the so-called point-mass solution was forwarded, which considers terms neglected in the integration of the surface integral (see Equation (4)) for the calculation of the gravitational potential. Because the surface integral cannot be solved in closed analytical form, a Taylor series expansion to its integrand was applied. By evaluating the series at a point equivalent to the geometric center (Q 0 ) of the spherical panels, the number of non-vanishing terms was reduced because only even terms remained after subsequent integration as shown by Equation (11). Furthermore, the zeroth-order term of the series corresponds to the OPM solution, and the additional higher-order terms represent the deviation of a panel from a point-mass. The results summarized in Table 1 show that the OPM modeling is within the range of the uncertainties expected for GRACE-derived e.w.h. fields. For example, Scanlon et al. [27] have reported uncertainties of GRACE-derived e.w.h. in the range of 5 to 40 mm for river basins larger than 100 × 10 3 km 2 (for smaller basins, these values are slightly higher). Groh et al. [59] have estimated uncertainties ranging from approx. 5.9 to 8.1 mm for AIS (considering AIS's surface area as approx. 12 × 10 6 km 2 ). Although the differences in the period and release of the GRACE Level-2 datasets between the studies, the uncertainties of GRACE solutions seem to be less than 10 mm at AIS. However, artifacts can be detected in the global map of the RMSEs (not shown here) of the inverted mass changes regarding OPM. For instance, Figure 4a shows a pattern of the RMSEs where large values concentrate along the AIS's coastlines (the same holds for all landmasses).
As the spatial patterns might be a consequence of edge effects in the inversion, computations considering parameters covering the whole Earth reduced the edge effects to a certain extent. For instance, the RMSE reduced in approximately 90% (Table 1) while solving the inverse problem using the OPM modeling considering the edge effects. Furthermore, considering higher-order terms in Equation (10), improve the point-mass (i.e., IPM) approach. For instance, considering terms up to third-order reduced the edge effects ( Figure 4) as well as the RMSE values (Table 1). In this study, it was considered only a third-order expansion of the surface integral (Equation (10)), and the contribution of higher-order terms (if any) was not investigated. Nevertheless, analyzing the remainder of the series expansion of the kernel given by Equation (5), performed at the Appendix B (see Equation (A13)), it is shown that the relative error due to the fourth-order approximation in terms of forwarded gravitational potential is less than 0.01%. (These figures are approximated 0.10% and 0.02% for the zeroth-and second-order approximations, respectively.) Therefore, objectives (i) and (ii) proposed in Section 1 were met. Moreover, with or without accounting for edge effects, the IPM solution went one better than the OPM solution ( Figure 4). The contribution of this exercise lies in the regional application of OPM modeling, for example, for estimating mass changes over Antarctica, which requires adding buffer zones in the inversion, which has not been discussed in the previous studies (e.g., Reference [38]).
Consequently, the IPM modeling provides an alternative procedure to recover AIS's mass changes using synthesized gravitational potential fields projected at GRACE's altitude as wished-for in objective (iii). Thus, in Section 3.2, it is presented the mass changes over AIS based only on IPM. At this stage, OPM was not considered in the inversion of Level-2 data since the comparison between both mass parametrizations would be necessary to apply a low-pass filter to the models in order to make them spectrally consistent with the synthesized residual gravitational potentials, which are band-limited up to d/o 120. This is important since OPM and IPM parametrizations could bias the short wavelengths of the estimated mass changes in different ways. It will be investigated in a follow-up study. In any case, the preliminary results depicted in Section 3.2 indicate that IPM ranks the second best with an RMSE of 5.77 mm. However, if CSR mascon solution (CSR-M) provides slightly better results over AIS (see Figure 6d), why should one consider the IPM modeling (the same holds for the OPM modeling) as a GRACE processing approach? Andrews et al. [26] and Scanlon et al. [27], among other authors, have shown that unconstrained mascon and unfiltered spherical harmonic solutions have the same variances at different frequencies. This might indicate that both solutions provide the same results for inverted mass changes, in which the differences are consequences of the constrain in the mascon solutions and filtering in the spherical harmonic solutions. Generally, when the whole time series of GRACE raw data are reprocessed in order to implement new corrections and/or reductions based on improved de-aliasing products, the first friendly solutions available to the public are Level-2 data. However, there is typically a large latency before new releases of corresponding mascon solutions are provided. Thus, as the IPM method is based directly on Level-2 data, users immediately benefit from new releases to obtain improvements in the estimated mass changes.
Furthermore, while there are three centers that provide mascon solutions to the public [15,20,25], there are numerous research groups (see, e.g., References [3][4][5][6][7][8][9][10][11][12][13][14]) from three continents that generate Level-2 data based on different approaches and methods in order to make the most out of the GRACE measurements. Additionally, the Level-2 data from the measurements of the state-of-art satellite gravimetry GRACE Follow-On (GRACE-FO) is already available, extending the legacy of the former GRACE mission. Mascon solutions based on GRACE-FO might show up anytime, and, meanwhile, the users can use GRACE-FO Level-2 data in the framework of IPM modeling (or the OPM if considered the edge effects) as a viable processing methodology akin to mascon solutions as shown for AIS (Figure 6d) The mass change estimations over AIS considered only the Level-2 dataset represented by the ITSG-Grace2018 solution. The reason for this was that the assessment of ITSG-Grace2018 conducted by Kvas et al. [8] has shown up 46% lower noise level than the Release 06 solutions produced by CSR, GFZ, and JPL over the quite ocean areas. Furthermore, these authors have investigated the performance of ITSG-Grace2018 over 33 river basins (basins at which the noise level does not affect the signal), and it has been reported that ITSG-Grace2018 exhibits the same signal content like CSR, GFZ, and JPL with variability at the range of ±2 mm. This indicates that the choice of one or another solution would not impact the finds in Section 3.2 up to a certain point. Nevertheless, the inversion using all (or most of them) Level-2 solutions produced by other processing centers and institutions could lead to a higher generalization of the results. Although this is a small study focusing on the entire AIS, the results in Section 3.2 cannot be generalized to the sub-domains of AIS, for example, West AIS (or other regions). Hence, a global assessment of the inversion of mass changes parameterized by IPM still necessary. To this end, the inversion at a global scale using Tikhonov-regularized least-squares may be computationally expensive for the most of users since it depends on the computation of the decomposition matrices. Solutions to the least-squares problem may be calculated very easily using conjugated gradient-like methods on the normal equations (e.g., Reference [62]). Alternatively, one could consider a regional application at which the parameters conveying the target region (or basin) plus a buffer zone could be estimated. Although the results presented here are based on Level-2 data, the IPM can be used in the same fashion way to invert Level-1b data as, for example, the approach used by JPL-M (cf. Reference [15]).

Conclusions
In this paper, it was investigated the overall performance of the original point-mass (OPM) formalism and proposed an improvement (IPM) as well. For instance, a closed-loop simulation using inverse modeling showed that the IPM solution outperformed the results based on OPM modeling at global and regional scales. The OPM is more susceptible to the edge effects in the inversion than IPM, which are concentrated along the coastlines. For example, the RMSE values along the coastal regions are reduced to a large extent while parametrizing the surface densities with IPM modeling. That is, considering the inversion over the entire Earth, it upgrades the OPM's performance to the level of that of the IPM solution. This means that the application of OPM requires an extended area beyond the target region's domain (for example, five arc-degrees). Regarding the AIS's mass changes, the overall patterns of the linear trends of GRACE-era mass changes considering the IPM modeling were estimated at which a linear trend of −7.58 mm/year are in agreement with those of the mascon solutions (CSR-M, GSFC-M, and JPL-M). Noise assessment of the IPM solution has shown RMSE values at the same level as the mascon solutions. Therefore, it is concluded that all the variabilities of the different GRACE solutions analyzed here indicate the level of noise in GRACE data over AIS. Albeit mascon solutions are available for the public and being a friendly product for the most of users, the IPM associated with synthesized potential from Level-2 data can be seen as a fast and quick means to assess mass changes at regional or global scales. of China (2018YFA0605402). V.G. Ferreira thanks for the support from National Natural Science Foundation of China (Grant No. 41574001). The ITSG-Grace2018 GRACE Level-2 data used in this study can be find online as obtained from the Institute of Geodesy of the Graz University of Technology (https://www.tugraz.at/ institute/ifg/downloads/gravity-field-models/itsg-grace2018/), the relevant references have been provided in Section 2.2. The degree-1 and degree-2 coefficients, as well as the GIA corrections and the Tellus land grids (GFZ) were retrieved from GRACE Tellus website (https://grace.jpl.nasa.gov/data/get-data/) and further details can be found in Section 2.2 regarding the relevant references. Furthermore, mascon solutions data described in Section 2.2 were obtained from the CSR GRACE RL06 (http://www2.csr.utexas.edu/grace/RL06_mascons.html), JPL GRACE/GRACE-FO RL06v2 (http://grace.jpl.nasa.gov), and GSFC GSFC.glb.200301_201607_v02.4-GeruoA (https://earth.gsfc.nasa.gov/geo/data/grace-mascons) and the respective references are available therein. V.G. Ferreira also thanks Tengke Sun for the discussions during the early version of this study. The authors thank the three anonymous reviewers for their valuable comments that helped the quality of the final manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. The Surface Integral via Analytical Integration over the Radial Direction
The Newton integral can be subdivided into a vertical 1-D integral (i.e., a radial integral) and a surface integral as: where dσ = cos ϕ dϕ dλ is the surface element, and σ is the parameter domain of the unit sphere. The Euclidean distance is given by Equation (2). The integral in brackets is an integral of the type R x, ax 2 + bx + c dx, which is expressible in terms of elementary functions. The solution of such integrals is of type [63] (p. 1035): with X = ax 2 + bx + c and a > 0 [63] (p. 1034). Hence, considering Equations (A3) and (A4), the radial integral within Equation (A1) gives: R+∆r R r 2 dr = 1 2 (r + 3r cos ψ) + r 2 (3cos 2 ψ − 1) × ln + r − r cos ψ R+∆r R , (A5) and the potential in Equation (A1) results in where N is the integral kernel according to Equation (A5), that is, N(ϕ , λ ) =R( − ) + ∆r + 3r cos ψ( − ) + r 2 (3cos 2 ψ − 1) ln + R + ∆r − r cos ψ + R − r cos ψ , where ∆r is the e.w.h. responsible for the gravitational attraction V, and and are the Euclidean distances given as: = r 2 + (R + ∆r) 2 − 2r(R + ∆r) cos ψ (A8) respectively. The kernel (Equation (A7)) can be developed in Taylor series w.r.t. ∆r. Since ∆r/R is relatively small in hydrological applications, that is, in extreme cases where ∆r reaches 0.60 m (at annual time scales) and considering the mean radius of the Earth as 6378 km, it gives |∆r/R | < 10 −7 , and one may neglect the nonlinear terms in this series expansion. Since N is the integral of the function r 2 / w.r.t. r (see Equation (A5)), the first-order derivative of N is simply r 2 / . Therefore, Equation (A7) becomes where the Taylor point r = R has been chosen; hence, becomes , as given by Equation (A9).
Equation (A13) gives what is left over after evaluating the n-th order for some (ϕ , λ ). Here (ξ, η) is an intermediate point line segment joining the two points (ϕ 0 , λ 0 ) and (ϕ , λ ). That is, the difference between the values of the function K at the points (ϕ , λ ) and (ϕ 0 , λ 0 ) is just the result for the differential at an intermediate point (ξ, η) on the line segment joining the two points.