Global Glacial Isostatic Adjustment Constrained by GPS Measurements: Spherical Harmonic Analyses of Uplifts and Geopotential Variations

: In addition to studies of sea level change and mantle rheology, reliable Glacial Isostatic Adjustment (GIA) models are necessary as a background model to correct the widely used Gravity Recovery and Climate Experiment (GRACE) monthly gravity solutions to determine subsecular, nonviscous variations. Based on spherical harmonic analyses, we developed a method using degree-dependent weighting to assimilate the Global Positioning System (GPS) derived crustal uplift rates into GIA model predictions, in which the good global pattern of GIA model predictions and better local resolution of GPS solutions are both retained. Some systematic errors in global GPS uplift rates were also corrected during the spherical harmonic analyses. Further, we used the reﬁned GIA uplift rates to infer the GIA-induced rates of Stokes coe ﬃ cients (complete to degree / order 120) relying on the accurate relationship between GIA vertical surface deformation and gravitational potential changes. The results show notable improvements relative to GIA model outputs, and may serve as a GIA-correction model for GRACE time-variable gravity data.


Introduction
During the last glacial maximal (~18,000 years ago), large portions of the Northern and Southern Hemispheres were covered with great sheets of ice, which caused isostatic depressions of the land below and around the ice but bulges surrounding these indentations. Though the ice retreated long ago, much of North America, Scandinavia, Greenland, and Antarctica are still rising where the massive layers of ice pushed it down, while the bulges are still slowly sinking. This ongoing movement of land is a viscous response to its ice-age burden and is called the Glacial Isostatic Adjustment (GIA). As an adjustment process of the Earth towards an equilibrium state with retreated ice, GIA is evident in various phenomena, such as geopotential variations, Earth surface displacements, global sea-level change, geocenter variations, and Earth's rotational changes, etc., which have been studied to infer the extent and amount of the former ice masses, to reconstruct the sea level during a glacial cycle, and to constrain rheological properties of the Earth's interior [1][2][3][4][5][6][7][8].
The GIA-induced vertical crustal deformations, or crust uplifts, cause mass redistributions within the Earth and thus lead to changes in geopotential. On decadal and shorter time scales, GIA-induced changes in geopotential and crustal uplift can both be regarded as linear drifts. Therefore, we are basis of 14 C dating [13,14]. The ICE-6G_D (VM5a) was a revised version of ICE-6G_C (VM5a), by replacing the NOCS ETOPO2 bathymetry model with the BEDMAP2 and rerunning their GIA codes [15]. Hereafter, the two models, respectively, are referred to as ICE6GC and ICE6GD for convenience. One can refer to Figure 1a,b for the uplift rates predicted, respectively, by the two models.
Using the Australian National University (ANU) group's CALSEA software package, as well as the ICE-6G ice thickness history and VM5 Earth rheology as the inputs, Purcell et al. [28,29] recalculated the present-day GIA-induced uplifts and geopotential changes. Although sharing the same models of the deglaciation history and mantle viscosity with ICE-6G_C/D, Purcell et al. adopted different topographic data which are a combination of the GEBCO_08 topographic database (version 20100927, http://www.gebco.net) and the BEDMAP2 data sets [30]. However, the model data released by Purcell et al. [28,29] are actually not the same (see Figure 1c,d). Hereafter, we termed these two ICE-6G variants as PUR16 and PUR18, respectively. The PUR16 and PUR18 data sets are available as supporting information of [28] and [29], respectively.
Unfortunately, all the four model outputs provide no information on their errors. However, one can see some notable differences mainly around Canada, Greenland, and Antarctic from Figure 1.

GPS Data and 2D Fitting
The data for uplift rates (together with their uncertainties) released by Schumacher et al. [16] are derived from 4072 GPS sites (the sites were selected based on prior information from the GIA forward models to exclude tectonic signals; see Figure 2 for the site distributions over the globe), and tested with 13 global GIA forward model solutions. This data set is derived from the global GPS data set from the Nevada Geodetic Laboratory (NGL) at the University of Nevada, Reno (UNR), as well as

GPS Data and 2D Fitting
The data for uplift rates (together with their uncertainties) released by Schumacher et al. [16] are derived from 4072 GPS sites (the sites were selected based on prior information from the GIA forward models to exclude tectonic signals; see Figure 2 for the site distributions over the globe), and tested with 13 global GIA forward model solutions. This data set is derived from the global GPS data set from the Nevada Geodetic Laboratory (NGL) at the University of Nevada, Reno (UNR), as well as the A-NET and G-NET regional GPS data sets, respectively, for Antarctica and Greenland, with outliers, jumps, pole tides, and elastic response of the solid Earth to global changes in ice sheets, glaciers, and atmospheric-ocean loading all corrected; further, these GPS data were converted from the center of mass of the total Earth's system (CM) reference frame to the centre of mass of the solid Earth (CE) frame (more details can be found in [16,31] glaciers, and atmospheric-ocean loading all corrected; further, these GPS data were converted from the center of mass of the total Earth's system (CM) reference frame to the centre of mass of the solid Earth (CE) frame (more details can be found in [16,31], while the corresponding GPS data set for GIA uplift rates is available at https://doi.pangaea.de/10.1594/PANGAEA.889923). Due to uncertainties in the GIA models and possible errors in GPS measurements and data processing, the agreement between the model predictions and the uplifts derived from GPS data (see Figure 3 for an example) is limited. In order to make the best of the GPS-derived uplift rates, we applied a 2D cubic-spline weighted fitting to them using MATLAB and obtained the corresponding global gridded data (the weights are determined from the uncertainties of the GPS solutions; see Figure 4a for global grids and Figure 4b for land-only grids), which made the following spherical harmonic expansion more convenient, and removed some outliers as shown in Figure 3. Hereafter, we termed these gridded data as GPSinterp. Due to uncertainties in the GIA models and possible errors in GPS measurements and data processing, the agreement between the model predictions and the uplifts derived from GPS data (see Figure 3 for an example) is limited. In order to make the best of the GPS-derived uplift rates, we applied a 2D cubic-spline weighted fitting to them using MATLAB and obtained the corresponding global gridded data (the weights are determined from the uncertainties of the GPS solutions; see Figure 4a for global grids and Figure 4b for land-only grids), which made the following spherical harmonic expansion more convenient, and removed some outliers as shown in Figure 3. Hereafter, we termed these gridded data as GPSinterp.

Weighting the GIA Models
There are some differences among the GIA models as shown in Figure 1, due to the adoption of different bathymetry models and/or numerical codes. Further, the notable differences between PUR16 and PUR18 should be caused by some changes of the relevant parameters. In order to facilitate the refinement of GIA effects considering the fact that all these model outputs provide no information on their errors, a weighted average of these models will be used as the reference model, which is termed as ICE6Gavg hereafter. That is because the weighted average would be the best in most cases due to the theory of errors. Figure 3 for an example) is limited. In order to make the best of the GPS-derived uplift rates, we applied a 2D cubic-spline weighted fitting to them using MATLAB and obtained the corresponding global gridded data (the weights are determined from the uncertainties of the GPS solutions; see Figure 4a for global grids and Figure 4b for land-only grids), which made the following spherical harmonic expansion more convenient, and removed some outliers as shown in Figure 3. Hereafter, we termed these gridded data as GPSinterp.

Weighting the GIA Models
There are some differences among the GIA models as shown in Figure 1, due to the adoption of different bathymetry models and/or numerical codes. Further, the notable differences between PUR16 and PUR18 should be caused by some changes of the relevant parameters. In order to facilitate the refinement of GIA effects considering the fact that all these model outputs provide no information on their errors, a weighted average of these models will be used as the reference model, which is termed as ICE6Gavg hereafter. That is because the weighted average would be the best in most cases due to the theory of errors.
We noted from Figure 4b that after applying the landmask to GPSinterp, the GPS results are then comparable with the GIA model predictions as shown in Figure 1. As direct measurements of crust uplift rates, GPSinterp within land areas may be used to weight the GIA models through: where RMS[x] means the root mean square of x and the super script L indicates that only the land uplifts are considered while oceanic ones are ignored. The landmask is a generalized function which equals 1 over land areas and 0 for all other regions. For MOD = ICE6GD, PUR16, or PUR18, the corresponding RMS differences with respect to We noted from Figure 4b that after applying the landmask to GPSinterp, the GPS results are then comparable with the GIA model predictions as shown in Figure 1. As direct measurements of crust uplift rates, GPSinterp within land areas may be used to weight the GIA models through: where RMS[x] means the root mean square of x and the super script L indicates that only the land uplifts are considered while oceanic ones are ignored. The landmask is a generalized function which equals 1 over land areas and 0 for all other regions. For MOD = ICE6GD, PUR16, or PUR18, the corresponding RMS differences with respect to GPSinterp are 1.2914, 1.2519, and 1.2986 mm/year, respectively. Then, their relative weights are 0.3275, 0.3485, and 0.3239, respectively. The ICE6GC was not used as it is based on a different bathymetry model and contains a notable error near the Antarctic coast [28,29].
Remote Sens. 2020, 12, 1209 6 of 18 The uplift rates for the weighted average ICE6Gavg model are illustrated in Figure 5, while uplift rate differences of ICE6GD, PUR16, PUR18, and GPSinterp with respect to ICE6Gavg are illustrated in Figure 6. One can see that all of them agree with each other within a few mm/year even at the areas with maximal differences (note that strong GPSinterp signals over oceans are not realistic, and one should focus on its signals on land only). These comparisons confirm that ICE6Gavg should be reliable as the reference model, though the discrepancies between ICE6Gavg and GPSinterp are larger. That fact should be due to possible uncertainty of GPS data and errors in the models for ice history and mantle viscosity.   (b) difference between PUR16 and ICE6Gavg; (c) difference between PUR18 and ICE6Gavg; (d) difference between GPSinterp and ICE6Gavg, where the differences are reconstructed using spherical harmonic analyses. The zero-difference contours are indicated by white curves. Note that strong GPSinterp signals over oceans are not realistic, thus only its signals on land make sense and will be used.

Spherical Harmonic Expansion of the Crustal Uplift Rates
By comparing Figures 1-6, we can see notable differences between the GIA uplift rates derived from model predictions and GPS-based results, and GIA models are good at describing large-scale and especially global GIA patterns while GPS data are better at local details although the GPS rates can also be affected by unrelated local effects [17,18]. Bearing this in mind, we developed the following spherical harmonic analyses to assimilate GPS data into GIA model outputs (predictions).
The field of GIA-induced uplift rates can be expanded to a series of spherical harmonics, namely: .
where a is the mean equatorial radius, θ and λ, respectively, are the colatitude and east-longitude, P nm (cos θ) are the fully-normalized associated Legendre polynomials, and . c X nm and . s X nm are the coefficients for uplift rates and can be obtained by [32] through numerical integration. In Equations (2) and (3), the superscript X denotes MOD (for a certain GIA model, ICE6GC/D or PUR16/18 for the current study) or GPS (for GPS-derived GIA uplift rates).
In this study, all spherical harmonic coefficients are truncated to degree/order 120 in order to match the monthly GRACE gravity models.

Degree-Dependent Weighting
As GPS-based results tend to provide biased low-degree terms but are promising to provide better high-degree terms which describe local patterns, we may design a degree-dependent weight function w n so that the adjusted coefficients can be written as: where the superscript GPS = GPSinterp and MOD = ICE6Gavg. It is obvious that w n should increase gradually with degree n and satisfy the condition: Here, we proposed the following weight function which can meet the above conditions and produce acceptable results: where α is a parameter that controls the increasing rate of w n with respect to n. The curves of w n (α) for selected α are plotted in Figure 7. We noted from Figures 2-4 that the GPS results are adversely affected by the nonuniform distributions of GPS sites. In general, the GPS uplift rates for the United States, South America, and West Europe are much better than those for other areas (hereafter, A-areas and B-areas, respectively). In addition, GPSinterp over ocean regions can be very poor and thus, should be dropped. We applied the landmask to remove the GPSinterp signals over oceans. With landmask applied, we found it rather good to set α = 2 and α = 0.5 for the A-and B-areas, respectively. This choice not only helped make the best of the A-area GPS uplift rates, but also attenuated the anomalous B-area GPS results.
Here, we would like to stress some facts about the GIA signals in Greenland, which is in the Bareas though they are very important in GIA studies. First, as shown in Figure 2, almost all GPS stations on Greenland distribute only along the coast. Second, there are complex geodynamic phenomena (e.g., rising mantle plumes [1,7]) going on for Greenland which makes it difficult just relying on GPS data. Therefore, we could not obtain a notable improvement of Greenland GIA in this study.
Taking into account Equation (4), this procedure can be done by: We noted from Figures 2-4 that the GPS results are adversely affected by the nonuniform distributions of GPS sites. In general, the GPS uplift rates for the United States, South America, and West Europe are much better than those for other areas (hereafter, A-areas and B-areas, respectively). In addition, GPSinterp over ocean regions can be very poor and thus, should be dropped. We applied the landmask to remove the GPSinterp signals over oceans. With landmask applied, we found it rather good to set α = 2 and α = 0.5 for the A-and B-areas, respectively. This choice not only helped make the best of the A-area GPS uplift rates, but also attenuated the anomalous B-area GPS results.
Here, we would like to stress some facts about the GIA signals in Greenland, which is in the B-areas though they are very important in GIA studies. First, as shown in Figure 2, almost all GPS stations on Greenland distribute only along the coast. Second, there are complex geodynamic phenomena (e.g., rising mantle plumes [1,7]) going on for Greenland which makes it difficult just relying on GPS data. Therefore, we could not obtain a notable improvement of Greenland GIA in this study.
Taking into account Equation (4), this procedure can be done by: where L A [x] or L B [x] means applying the A-or B-landmask to x, and α 1 = 2, α 2 = 0.5 as discussed above, and is the degree-n component of . u X (X = GPS or ICE6Gavg). One can see that the systematic errors in GPS are automatically removed by our spherical harmonic analyses (Equations (6)- (8), where the (0,0) components of GPSinterp have zero weights) but may not be as effectively handled by other spatial filtering methods, such as empirical orthogonal function filtering or some Gaussian filtering [33][34][35][36][37].
The adjusted model is termed as ICE6Gavg (GPS).

Geopotential Variations Due to GIA
The rate of geopotential changes due to GIA can be expressed as [32,38]: where G is the gravitational constant, M and a are the mass and mean equatorial radius of the Earth, respectively, C nm (t) and S nm (t) are time-dependent fully-normalized Stokes coefficients.
A linear relation is generally assumed between . u GIA nm and δ . V GIA nm , that is: where R(n) can be obtained from a certain Earth model, such as PREM, with some theoretical approximations.
Wahr et al. [20,21] argued that for each degree n, the viscoelastic mode M0, primarily caused by buoyancy forces acting on the depressed lithosphere, dominates the total signals of crustal uplift rates. Then, they found that for a variety of possible ice geometries and time histories for both Antarctic and Greenland, as well as plausible mantle viscosity profiles, R(n) can be approximately obtained as: if only the effect of the M0 mode is considered. Wahr et al. [20] found that this approximation was able to recreate an uplift rate field for Antarctica to within ∼2 mm/yr compared to a field generated by their ice sheet modeling program. This empirical approximation has been used by van der Wal et al. [23], Tregoning et al. [24], and Wu et al. [25] for deriving GIA uplift rates from GRACE temporal spherical harmonic fields.
On the other hand, rather than relying on an approximate mathematical formulation restricted to the M0 mode, Purcell et al. [22] considered the rheology of the whole Earth and preferred to express the GIA-induced uplift as: which implies a more rigorous relation between the GIA-induced crustal uplift and the associated geopotential change: where h ve n and k ve n are degree-n viscoelastic load Love numbers. Purcell et al. [22] showed that Equation (13) is not sensitive to uncertainties in ice or Earth models and valid for a broad range of Earth and ice-load models, and uplift rates obtained by Equation (13) agree with the results of fully detailed forward modeling routines within 0.3 mm/yr, a significant improvement over the results of the Wahr et al. [20,21]. The R(n) functions for the Wahr and Purcell models are compared in Figure 8.  From Figure 8, one can see a notable difference between the two models and the difference increases with the degree n. Therefore, we chose the more accurate Purcell model but not the Wahr model as adopted by previous studies [23][24][25].

Improved GIA Model ICE6Gavg (GPS) for Uplift Rates by Assimilating GPS Data
Based on harmonic analyses as described by Equations (2) and (3), we obtained the coefficients X nm c  and X nm s  as described by Figure 9, as well as Tables 1 and 2. In Figures 9 and 12, the line number is relevant with the tabulating of degree n and order m (refer to Tables 3 and 4 for examples), and its definition is consistent with the rows in the data section of GRACE monthly gravity models. For a given degree n, the maximal of the line number is (n+1)(n+2)/2.  From Figure 8, one can see a notable difference between the two models and the difference increases with the degree n. Therefore, we chose the more accurate Purcell model but not the Wahr model as adopted by previous studies [23][24][25].

Improved GIA Model ICE6Gavg (GPS) for Uplift Rates by Assimilating GPS Data
Based on harmonic analyses as described by Equations (2) and (3), we obtained the coefficients . c X nm and . s X nm as described by Figure 9, as well as Tables 1 and 2. In Figure 9 and Figure 12, the line number is relevant with the tabulating of degree n and order m (refer to Tables 3 and 4 for examples), and its definition is consistent with the rows in the data section of GRACE monthly gravity models. For a given degree n, the maximal of the line number is (n+1)(n+2)/2.  Table S1 for the full table with coefficients with more digits up to degree and order 120.  From Tables 1 and 2 and Figure 9, one can see that the low-degree terms for GPSinterp are highly biased compared to ICE6GC/D and PUR16/18. Especially, the nonzero and positive (0,0) term, namely 0 0 0 0 X u ac =   , implies a systematic error in the global GPS solutions for uplift rates. That is, the Earth was expanding as a whole if 00 u  has the GPSinterp value as shown in Table 1, while it is more likely the Earth conserves its volume during the GIA process. Therefore, we need to apply to GPSinterp the degree-dependent weighting as described in Section 2.2.2, to eliminate or reduce the low-degree anomalies but retain some local detailed signals conveyed by higher-degree terms. By assimilating the GPS data (GPSinterp) into the averaged GIA model ICE6Gavg through Equations (6)-(8), we obtain the adjusted model for GIA uplift rates, namely ICE6Gavg (GPS), of which the adjusted coefficients are described in Figure 9 and Tables 1 and 2. Moreover, the corresponding adjusted global GIA uplift rates can be obtained from these adjusted coefficients by using Equation (2). From Figures 10 and 11, one can see that the adjusted uplifts in the areas with dense GPS sites got refined markedly (especially for the A-areas) with respect to the original model predictions. For other areas, GPS can only show some departures from GIA model predictions. From Tables 1 and 2 and Figure 9, one can see that the low-degree terms for GPSinterp are highly biased compared to ICE6GC/D and PUR16/18. Especially, the nonzero and positive (0,0) term, namely . u 00 = a . c X 00 , implies a systematic error in the global GPS solutions for uplift rates. That is, the Earth was expanding as a whole if . u 00 has the GPSinterp value as shown in Table 1, while it is more likely the Earth conserves its volume during the GIA process. Therefore, we need to apply to GPSinterp the degree-dependent weighting as described in Section 2.2.2, to eliminate or reduce the low-degree anomalies but retain some local detailed signals conveyed by higher-degree terms.
By assimilating the GPS data (GPSinterp) into the averaged GIA model ICE6Gavg through Equations (6)-(8), we obtain the adjusted model for GIA uplift rates, namely ICE6Gavg (GPS), of which the adjusted coefficients are described in Figure 9 and Tables 1 and 2. Moreover, the corresponding adjusted global GIA uplift rates can be obtained from these adjusted coefficients by using Equation (2). From Figures 10 and 11, one can see that the adjusted uplifts in the areas with dense GPS sites got refined markedly (especially for the A-areas) with respect to the original model predictions. For other areas, GPS can only show some departures from GIA model predictions.

Improved GIA-Induced Geopotential Variations
In this study, we chose R(n) as given by Equation (13) as it is more accurate than Equation (11).
Then, the geopotential changes δ s ADJ nm by Equation (10). The rates of Stokes coefficients for ICE6GC/D, PUR16/18, and GPSinterp can also be obtained in a similar way (see Figure 12 and Tables 3 and 4 for some comparisons).
From Figure 12 and Tables 3 and 4

Improved GIA-Induced Geopotential Variations
In this study, we chose R(n) as given by Equation (13) as it is more accurate than Equation (11). Then, the geopotential changes

Improved GIA-Induced Geopotential Variations
In this study, we chose R(n) as given by Equation (13) as it is more accurate than Equation (11). Then, the geopotential changes    Table S4 for the full table with coefficients with more digits up to degree and order 120.

Discussion
Recall that this study aims to provide an improved GIA-correction for GRACE monthly gravity models, whose spatial resolution is about 300 by 300 km (around 3 by 3 degrees at the equator). Therefore, the distances between every two GPS sites should be around 300 km (or 3 degrees) or less to match the resolution of GRACE. One can see that only the A-areas defined in Section 2.2 can meet this condition, which can justify the different treatments of the A-and B-areas adopted by this study. In fact, GPS may play a much more important role to study GIA-related issues provided GPS sites distributions get improved in the future.

Discussion
Recall that this study aims to provide an improved GIA-correction for GRACE monthly gravity models, whose spatial resolution is about 300 by 300 km (around 3 by 3 degrees at the equator). Therefore, the distances between every two GPS sites should be around 300 km (or 3 degrees) or less to match the resolution of GRACE. One can see that only the A-areas defined in Section 2.2 can meet this condition, which can justify the different treatments of the A-and B-areas adopted by this study. In fact, GPS may play a much more important role to study GIA-related issues provided GPS sites distributions get improved in the future.
On the other hand, as shown in Figures 2-4, though the global crustal uplift rates derived from GPS data (GPSinterp) may provide more local details, they are adversely affected by the nonuniform distributions of GPS sites, and can be very poor over ocean regions. Thus, the GPS uplift rates for the dense GPS site areas, namely the A-areas (including United States, South America, and West Europe), should have a higher weight (α = 2), and those for the less-dense GPS site areas, namely the B-areas (including all other lands), should have a lower weight (α = 0.5), while those for ocean regions should be at zero weight. This choice not only helped make the best of the A-area GPS uplift rates, but also attenuated the weaker B-area GPS results and removed the anomalous signals over oceans. All these were achieved by using Equations (6)- (8), which describe the full procedures of degree-dependent weighting. Admittedly, these α values are not well constrained since the data of another type (namely non-GPS) are needed to do so while only GPS can provide uplift rates at the global scale. However, for α >2, the resulting differences are not significant. Therefore, the results presented in this study should be reliable.
From Figures 10 and 11 and Tables 1 and 2, one can see that the . c ADJ nm and . s ADJ nm absorbed the merits of both GIA models and GPS data as designed. That is, the adjusted model ICE6Gavg (GPS) not only has a very similar global pattern ( Figure 10) as the three ICE6G variants, namely ICE6GD, PUR16, and PUR18, but also more details (Figure 11) than them. This proves the validity of the spherical harmonic analyses and degree-dependent weighting.
One may wonder why the low-degree coefficients for GPS uplift rates are incorrect. Physically speaking, . u 00 = a . c X 00 describes the uplift or sink rate of the Earth surface as a whole. Therefore, the GPS results as shown in the second row of Tables 1 and 2 and Figure 6d imply a net uplift rate at the global scale (Figure 6d shows that almost the entire Earth is uplifting, except for three small zones; if it were true, then the Earth is expanding as a whole), while it is more likely the Earth conserves its volume during the GIA processes. In addition, as shown below, coefficients . c X nm and . s X nm are related to the rates of Stokes coefficients, and thus these biased low degree coefficients . c X nm and . s X nm imply unrealistic changes in the Earth's total mass and inertial tensor, as well as geocenter, which are due to not only the nonuniform distributions of the GPS sites as shown in Figure 2, but also some systematic errors in the GPS solutions as indicated by Figure 6d.
Through Equation (13), the low-degree δ . C GPS nm and δ . S GPS nm are also problematic due to the significantly biased GPS-derived crustal uplift rates. Here, we stress the systematic errors of GPS uplift rates through the geophysical meanings of the degree-0, 1, and 2 Stokes coefficients: (1) The nonzero degree zero term implies a rather rapid change in the Earth's mass; (2) the anomalous degree one terms predict a wrong direction of geocenter motion; (3) the over-estimated . C 20 are more than twice of the observed values. The above reasons explained why the degree-dependent weight function should be applied. Therefore, the GPS-derived GIA uplifts are biased and should be corrected by using the method as described in Section 2.2.2.

Conclusions
In this study, we developed a method to combine and refine the crust uplift rates derived from GIA models and GPS observations, since GPS uplifts can be very good in some small regions but are very poor in global patterns, especially for the global systematic errors evidenced as the surface uplifting of the whole Earth, and the GIA model is on the contrary. This method is based on spherical harmonic analyses of both the modeled and observed crust uplift rates, and then a degree-dependent weighting of them. Therefore, our results absorbed the merits from both GIA models and GPS data, and removed the systematic errors in the global GPS crustal uplift rates.
Further, we derived the rates of Stokes coefficients based on the refined GIA uplift rates, as well as the accurate relationship between GIA uplift rates and gravitational potential changes proposed by Purcell et al. [22]. The results may be used as a GIA correction to various GRACE monthly geopotential models, especially the high-accuracy Mascon products as released by Watkins et al. [9], Wiese et al. [10] and Save et al. [11], and may help to generate an updated version of the multiple-data-based monthly geopotential model set LDCmgm90 as released by Chen et al. [12].
To sum, it may be useful to refine the resolution of the GIA uplift map by assimilating GPS data, through which the GIA-induced geopotential changes can also be improved. Doubtlessly, both GIA uplift rates and geopotential changes can be continuously improved provided there are more and more GPS sites distributed uniformly over the globe (especially for the B areas), with site distances around 3 degrees (about 300 km at the equator) or less. With the refined GIA-induced geopotential removed, GRACE data may be used to study other long-period variations in the Earth system.