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

Abstract

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 refined GIA uplift rates to infer the GIA-induced rates of Stokes coefficients (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.

1. 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 more interested in the rate of crust uplift^u˙GIA≡∂^uGIA(θ,λ,t)/∂tand rate of geopotential changeδ^V˙GIA≡∂[δ^VGIA(θ,λ,t)]/∂trather than^uGIAandδ^VGIAthemselves. Reliable^u˙GIAcan be used to infer mantle viscosity, andδ^V˙GIA are necessary as a background model to correct the widely used Gravity Recovery and Climate Experiment (GRACE) monthly gravity solutions for these secular variations [9,10,11,12].

Currently, the GIA-induced uplift rates can be obtained either from a certain GIA model [13,14,15] or from Global Positioning System (GPS) measurements [16]. However, the agreements between model-based and GPS-based results are usually limited due to a number of factors (some are discussed in the following sections). It is believed that the 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].

Arguing that the GIA-induced crust uplifts are dominated by the Earth’s viscosity mode M0 [19], Wahr et al. [20,21] proposed an approximate theory to infer changes in Stokes coefficients from the measured GIA uplift rates. Later, Purcell et al. [22] refined the theory relying on viscous load Love numbers (also see Section 2.3 for some details). While most studies [23,24,25] of the GIA-induced geopotential change are on the basis of the approximate theory of Wahr et al., Purcell et al. [22] and Jia et al. [26] pointed out that this method will lead to a ~15% relative error, which would weaken GRACE’s monitoring of global changes. However, both of the two methods ignore possible systematic errors in measured uplift rates as discussed in this study.

On the other hand, Argus et al. [13] and Peltier et al. [14,15] had assimilated GPS uplift rates into their GIA numerical simulations, leading to a notable improvement to their GIA models. However, only a few tens of GPS sites were used in their studies, which is obviously not enough for GIA as a global process. Guo et al. [27] compared 14 GIA models from different studies and concluded that the accuracy and consistency of GIA models need to be substantially improved to fully exploit space geodetic data, such as GRACE data, to enhance the constraints on an ice-sheet mass balance and the mass component of global sea-level change. Recently, Schumacher et al. [16] derived vertical crustal uplift rates from global GPS measurements and compared them with various GIA models, and also suggested notable discrepancies between different GIA models.

This study aims to refine the GIA uplift rates by assimilating global GPS results (but with some systematic errors corrected, as shown in Section 2), and then obtain refined rates of Stokes coefficients using 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 serve as a GIA-correction model for the GRACE time-variable gravity data (including Mascon solutions [9,10,11]), and can be improved in the future if the data of GPS-derived uplift rates are denser and more uniformly distributed over the globe.

2. Materials and Methods 2.1. Data and Models

In this study, we have adopted four variants of the ICE-6G GIA model, as well as a global GPS data set for GIA uplift rates recently released by Schumacher et al. [16] to obtain refined GIA-induced uplifts and geopotential changes. More details about these models and data can be found below.

2.1.1. GIA Models

ICE-6G_C (VM5a) and ICE-6G_D (VM5a) are the latest models of the glacial isostatic adjustment process in the ICE-NG (VMX) sequence from the University of Toronto [13,14,15]. New models of the deglaciation history and mantle viscosity were introduced by the two models (available at http://atmosp.physics.utoronto.ca/~peltier/data.php).

Based on the National Oceanography Centre, Southampton (NOCS), the ETOPO2 bathymetry model, the ICE-6G_C (VM5a) was constrained to fit the uplift rates at 42 sites estimated from GPS measurements, ice thickness change at 62 locations estimated from exposure-age dating, Holocene relative sea level histories from 12 locations estimated on the basis of radiocarbon dating, and age of the onset of marine sedimentation at nine locations along the Antarctic shelf also estimated on the basis of ¹⁴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.

2.1.2. 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], 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.

2.1.3. 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:

W_tMOD=1_σMOD2,σ=RMS[_vMODL−_vGPSinterpL],

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].

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.

2.2. Spherical Harmonic Analyses

2.2.1. Spherical Harmonic Expansion of the Crustal Uplift Rates

By comparing Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 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:

^u˙X(θ,λ)=a∑n=0∞∑m=0n[_c˙nmXcosmλ+_s˙nmXsinmλ]_Pnm(cosθ)

where a is the mean equatorial radius, θ and λ, respectively, are the colatitude and east-longitude,_Pnm(cosθ)are the fully-normalized associated Legendre polynomials, and_c˙nmXand_s˙nmX 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).

_c˙nmX=14π^a2∬S^u˙X(θ,λ)acosmλ_Pnm(cosθ)^a2dθdλ_s˙nmX=14π^a2∬S^u˙X(θ,λ)asinmλ_Pnm(cosθ)^a2dθdλ

In this study, all spherical harmonic coefficients are truncated to degree/order 120 in order to match the monthly GRACE gravity models.

2.2.2. 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_wnso that the adjusted coefficients can be written as:

(_{c˙nmADJs˙nmADJ})=_wn(_{c˙nmGPSs˙nmGPS})+(1−_wn)(_{c˙nmMODs˙nmMOD})=(_{c˙nmMODs˙nmMOD})+_wn(_c˙nmGPS−_{c˙nmMODs˙nmGPS}−_s˙nmMOD)

where the superscript GPS = GPSinterp and MOD = ICE6Gavg. It is obvious that_wnshould increase gradually with degree n and satisfy the condition:

{_wn=0,n=0_wn→1,n→+∞.

Here, we proposed the following weight function which can meet the above conditions and produce acceptable results:

_wn(α)=2π^tan−1αnπ180,

where α is a parameter that controls the increasing rate of_wnwith respect to n. The curves of_wn(α) for selected α are plotted in Figure 7.

We noted from Figure 2, Figure 3 and Figure 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:

^u˙ADJ=^u˙ICE6Gavg+_LA[∑n_wn(_α1)(_u˙nGPS−_u˙nICE6Gavg)]+_LB[∑n_wn(_α2)(_u˙nGPS−_u˙nICE6Gavg)]

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

_u˙nX(θ,λ)=a∑n=0∞[_c˙nmXcosmλ+_s˙nmXsinmλ]_Pnm(cosθ)

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).

2.3. Geopotential Variations Due to GIA

The rate of geopotential changes due to GIA can be expressed as [32,38]:

δ^V˙GIA(θ,λ,t)=GMr∑n=0∞^(ar)n∑m=0n[δ_C˙nmGIA(t)cosmλ+δ_S˙nmGIA(t)sinmλ]_Pnm(cosθ)

where G is the gravitational constant, M and a are the mass and mean equatorial radius of the Earth, respectively,_Cnm(t)and_Snm(t)are time-dependent fully-normalized Stokes coefficients.

A linear relation is generally assumed between_u˙nmGIAandδ_V˙nmGIA, that is:

[_{c˙nmGIAs˙nmGIA}]=R(n)[δ_C˙nmGIAδ_S˙nmGIA]

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:

^RWahr(n)=_u˙nmGIAδ_V˙nmGIA≈2n+12

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:

^u˙GIA(θ,λ)=a∑n=0∞ _{^h′nve ^k′nve}∑m=0n(δ_C˙nmGIAcosmλ+δ_S˙nmGIAsinmλ)_Pnm(cosθ)

which implies a more rigorous relation between the GIA-induced crustal uplift and the associated geopotential change:

^RPurcell(n)= _{^h′nve ^k′nve}

where _^h′nveand _^k′nve 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].

3. Results 3.1. 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˙nmXand_s˙nmX as described by Figure 9, as well as Table 1 and Table 2.

In Figure 9 and Figure 12, the line number is relevant with the tabulating of degree n and order m (refer to Table 3 and Table 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 Table 1 and Table 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˙00X, 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 Table 1 and Table 2. Moreover, the corresponding adjusted global GIA uplift rates can be obtained from these adjusted coefficients by using Equation (2). From Figure 10 and Figure 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.

3.2. 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δ_C˙nmADJandδ_S˙nmADJcan be obtained from the adjusted coefficients_c˙nmADJand_s˙nmADJ 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 Table 3 and Table 4 for some comparisons).

From Figure 12 and Table 3 and Table 4, one can see theδ_C˙nmADJandδ_S˙nmADJabsorbed the merits of both GIA models and GPS data as designed.

4. 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 Figure 2, Figure 3 and Figure 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 Figure 10 and Figure 11 and Table 1 and Table 2, one can see that the_c˙nmADJand_s˙nmADJ 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˙00X describes the uplift or sink rate of the Earth surface as a whole. Therefore, the GPS results as shown in the second row of Table 1 and Table 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˙nmXand_s˙nmXare related to the rates of Stokes coefficients, and thus these biased low degree coefficients_c˙nmXand_s˙nmX 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˙nmGPSandδ_S˙nmGPSare 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˙20are 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.

5. 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.

Figure 1. Glacial isostatic adjustment (GIA) uplift rates as predicted by four variants of the ICE-6G models. (a) ICE6GC; (b) ICE6GD; (c) PUR16; (d) PUR18. The zero-motion contours are indicated by white curves. Despite visible differences for some local details, all the four model predictions agree with each other at the global scale.

Figure 2. Distributions of the 4072 global positioning system (GPS) sites for the study of GIA crustal uplift rates.

Figure 3. Comparisons of GIA uplift rates from the ICE6GC model and GPS results. The ICE6GC uplift rates form a smooth surface while the discrete GPS results are similar to sharp 'spikes'.

Figure 4. Preliminary processing of the GPS-derived GIA uplift rates as provided by Schumacher et al. [13]. (a) GIA uplift rates from fitted GPS solutions; (b) the same as (a) but with the landmask applied. The zero-motion contours are indicated by white curves.

Figure 5. Uplift rates for the weighted average GIA model ICE6Gavg. The zero-motion contours are indicated by white curves.

Figure 6. Differences in uplift rates of ICE6GD, PUR16, PUR18, and GPSinterp with respect to the weighted average GIA model ICE6Gavg. (a) Difference between ICE6GD and ICE6Gavg; (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.

Figure 7. Degree-dependent weight wn for selected values of α.

Figure 8. R(n) functions for the Wahr and Purcell models, respectively.

Figure 9. Coefficientsc˙nmXands˙nmXfor the uplift rates as derived from ICE6GC/D, PUR16/18, GPSinterp, and ICE6Gavg (GPS). The coefficients for line number≤200 are displayed in small boxes to show more details.

Figure 10. Final results for GIA uplift rates: The ICE6Gavg (GPS) solution. The zero-motion contours are indicated by white curves.

Figure 11. Differences in uplift rates of ICE6Gavg (GPS) with respect to ICE6Gavg. The zero-difference contours are indicated by white curves.

Figure 12. Rates of Stokes coefficients due to GIA as derived from ICE6GC/D, PUR16/18, GPSinterp, and ICE6Gavg (GPS). The coefficients for line number≤200 are displayed in small boxes to show more details.

n, m	ICE6G-C	ICE6G-D	PUR16	PUR18	GPSinterp	ICE6Gavg (GPS)
0, 0	−0.0000	−0.0000	0.0000	0.0001	3.4973	0.0000
1, 0	0.9348	0.9700	1.3803	1.4527	2.3985	1.2694
1, 1	0.1554	0.1500	0.2337	0.2392	1.3595	0.2081
2, 0	4.0201	3.8636	3.7802	3.9165	6.7551	3.8517
2, 1	−0.5934	−0.6016	−0.2137	−0.2546	−3.3227	−0.3540
2, 2	−0.8500	−0.8470	−0.8839	−0.8875	−0.5817	−0.8730
3, 0	0.6715	0.8680	0.9535	1.0183	0.4444	0.9465
3, 1	0.3309	0.3130	0.3495	0.3135	3.0800	0.3259
3, 2	−1.5963	−1.5988	−1.6791	−1.7070	−2.9660	−1.6618
3, 3	0.0754	0.0756	0.0355	0.1017	2.2209	0.0700
4, 0	1.9807	1.7508	1.6660	1.7746	0.2638	1.7290
4, 1	0.3089	0.3510	0.3940	0.3507	5.1859	0.3659
4, 2	−2.6493	−2.6374	−2.7762	−2.9039	−0.8364	−2.7721
4, 3	0.1502	0.1423	0.1370	0.1643	−0.6214	0.1475
4, 4	0.1325	0.1352	0.1415	0.1552	3.6533	0.1439
5, 0	−2.5332	−2.2571	−2.2676	−2.4280	−4.6320	−2.3161
5, 1	0.6841	0.6520	0.7159	0.6802	1.2718	0.6834
5, 2	−3.1391	−3.1644	−3.3642	−3.5365	−2.5247	−3.3546
5, 3	−0.0278	−0.0229	−0.0560	0.0130	−0.4209	−0.0228
5, 4	0.9174	0.9055	0.9614	1.0851	0.8596	0.9832
5, 5	0.0034	0.0038	0.0082	0.0192	−1.8049	0.0103
6, 0	−0.2711	−0.5588	−0.7469	−0.8136	−1.0671	−0.7069
6, 1	0.5268	0.5785	0.6146	0.6067	0.4287	0.6002
6, 2	−3.8244	−3.8045	−3.9613	−4.0711	−3.3941	−3.9455
6, 3	−0.0672	−0.0759	−0.1010	0.0230	−0.4247	−0.0526
6, 4	1.2479	1.2476	1.3217	1.3811	−2.4670	1.3167
6, 5	0.0121	0.0165	0.0209	0.0123	1.4768	0.0167
6, 6	−0.0073	−0.0090	−0.0088	0.0025	−1.4879	−0.0052

1 See Table S1 for the full table with coefficients with more digits up to degree and order 120.

n, m	ICE6G-C	ICE6G-D	PUR16	PUR18	GPSinterp	ICE6Gavg (GPS)
0, 0	0	0	0	0	0	0
1, 0	0	0	0	0	0	0
1, 1	−0.3925	−0.3936	−0.5428	−0.5814	−0.2333	−0.5064
2, 0	0	0	0	0	0	0
2, 1	2.4902	2.5227	3.4276	3.6953	−2.0796	3.2179
2, 2	0.0732	0.0761	0.0752	0.0670	−1.3995	0.0728
3, 0	0	0	0	0	0	0
3, 1	−3.0718	−3.0483	−3.0044	−3.2498	−4.1242	−3.0983
3, 2	0.4888	0.4740	0.4718	0.5351	2.7229	0.4930
3, 3	0.6772	0.6705	0.6628	0.7411	−2.4413	0.6907
4, 0	0	0	0	0	0	0
4, 1	−2.7379	−2.8181	−2.9112	−3.1001	−3.2907	−2.9419
4, 2	0.3035	0.3135	0.3220	0.4008	0.0360	0.3447
4, 3	1.0494	1.0546	1.0758	1.1641	0.2258	1.0975
4, 4	0.3769	0.3635	0.4033	0.4435	3.4466	0.4033
5, 0	0	0	0	0	0	0
5, 1	−3.1112	−3.0053	−3.0913	−3.3399	−4.6518	−3.1437
5, 2	0.3057	0.2865	0.2900	0.4082	−1.1122	0.3271
5, 3	2.0740	2.0754	2.1677	2.3645	0.5202	2.2012
5, 4	0.0711	0.0750	0.0996	0.0729	1.7572	0.0829
5, 5	0.0839	0.0767	0.1060	0.0745	0.2649	0.0862
6, 0	0	0	0	0	0	0
6, 1	−0.5756	−0.7401	−0.8132	−0.7757	1.7817	−0.7771
6, 2	0.3991	0.4420	0.5284	0.6393	−1.4918	0.5360
6, 3	3.1258	3.1416	3.3148	3.4857	3.1034	3.3134
6, 4	0.2512	0.2517	0.3185	0.2648	−0.5003	0.2792
6, 5	−0.3714	−0.3669	−0.3886	−0.4056	0.5800	−0.3870
6, 6	0.0443	0.0427	0.0265	0.0340	−1.4207	0.0343

1 See Table S2 for the full table with coefficients with more digits up to degree and order 120.

n, m.	ICE6G-C	ICE6G-D	PUR16	PUR18	GPSinterp	ICE6Gavg (GPS)
0, 0	−0.0000	−0.0000	0.0001	0.0001	6.6450	0.0001
1, 0	0.5920	0.6143	0.8742	0.9201	1.5191	0.8039
1, 1	0.0984	0.0950	0.1480	0.1515	0.8610	0.1318
2, 0	1.4870	1.4291	1.3983	1.4487	2.4987	1.4247
2, 1	−0.2195	−0.2225	−0.0791	−0.0942	−1.2290	−0.1309
2, 2	−0.3144	−0.3133	−0.3270	−0.3283	−0.2152	−0.3229
3, 0	0.1732	0.2238	0.2459	0.2626	0.1146	0.2441
3, 1	0.0853	0.0807	0.0901	0.0808	0.7942	0.0840
3, 2	−0.4116	−0.4123	−0.4330	−0.4402	−0.7648	−0.4285
3, 3	0.0194	0.0195	0.0091	0.0262	0.5727	0.0181
4, 0	0.3890	0.3439	0.3272	0.3485	0.0518	0.3396
4, 1	0.0607	0.0689	0.0774	0.0689	1.0185	0.0719
4, 2	−0.5203	−0.5180	−0.5452	−0.5703	−0.1643	−0.5444
4, 3	0.0295	0.0279	0.0269	0.0323	−0.1220	0.0290
4, 4	0.0260	0.0266	0.0278	0.0305	0.7175	0.0283
5, 0	−0.4003	−0.3567	−0.3584	−0.3837	−0.7320	−0.3660
5, 1	0.1081	0.1030	0.1131	0.1075	0.2010	0.1080
5, 2	−0.4961	−0.5001	−0.5316	−0.5589	−0.3990	−0.5301
5, 3	−0.0044	−0.0036	−0.0088	0.0021	−0.0665	−0.0036
5, 4	0.1450	0.1431	0.1519	0.1715	0.1358	0.1554
5, 5	0.0005	0.0006	0.0013	0.0030	−0.2852	0.0016
6, 0	−0.0358	−0.0739	−0.0987	−0.1076	−0.1411	−0.0935
6, 1	0.0696	0.0765	0.0813	0.0802	0.0567	0.0794
6, 2	−0.5056	−0.5030	−0.5237	−0.5382	−0.4487	−0.5216
6, 3	−0.0089	−0.0100	−0.0134	0.0030	−0.0562	−0.0070
6, 4	0.1650	0.1649	0.1747	0.1826	−0.3262	0.1741
6, 5	0.0016	0.0022	0.0028	0.0016	0.1952	0.0022
6, 6	−0.0010	−0.0012	−0.0012	0.0003	−0.1967	−0.0007

1 See Table S3 for the full table with coefficients with more digits up to degree and order 120.

n, m	ICE6G-C	ICE6G-D	PUR16	PUR18	GPSinterp	ICE6Gavg (GPS)
0, 0	0	0	0	0	0	0
1, 0	0	0	0	0	0	0
1, 1	−0.2486	−0.2493	−0.3438	−0.3682	−0.1477	−0.3207
2, 0	0	0	0	0	0	0
2, 1	0.9211	0.9331	1.2679	1.3669	−0.7692	1.1903
2, 2	0.0271	0.0282	0.0278	0.0248	-0.5177	0.0269
3, 0	0	0	0	0	0	0
3, 1	−0.7921	−0.7861	−0.7747	−0.8380	−1.0635	−0.7990
3, 2	0.1261	0.1222	0.1217	0.1380	0.7022	0.1271
3, 3	0.1746	0.1729	0.1709	0.1911	−0.6295	0.1781
4, 0	0	0	0	0	0	0
4, 1	−0.5377	−0.5535	−0.5718	−0.6089	−0.6463	−0.5778
4, 2	0.0596	0.0616	0.0632	0.0787	0.0071	0.0677
4, 3	0.2061	0.2071	0.2113	0.2286	0.0443	0.2155
4, 4	0.0740	0.0714	0.0792	0.0871	0.6769	0.0792
5, 0	0	0	0	0	0	0
5, 1	−0.4917	−0.4749	−0.4885	−0.5278	−0.7351	−0.4968
5, 2	0.0483	0.0453	0.0458	0.0645	−0.1758	0.0517
5, 3	0.3278	0.3280	0.3426	0.3737	0.0822	0.3479
5, 4	0.0112	0.0119	0.0157	0.0115	0.2777	0.0131
5, 5	0.0133	0.0121	0.0168	0.0118	0.0419	0.0136
6, 0	0	0	0	0	0	0
6, 1	−0.0761	−0.0979	−0.1075	−0.1026	0.2356	−0.1027
6, 2	0.0528	0.0584	0.0699	0.0845	−0.1972	0.0709
6, 3	0.4133	0.4153	0.4382	0.4608	0.4103	0.4381
6, 4	0.0332	0.0333	0.0421	0.0350	−0.0661	0.0369
6, 5	−0.0491	−0.0485	−0.0514	−0.0536	0.0767	−0.0512
6, 6	0.0059	0.0056	0.0035	0.0045	−0.1878	0.0045

1 See Table S4 for the full table with coefficients with more digits up to degree and order 120.

Supplementary Materials

The following are available online at https://www.mdpi.com/2072-4292/12/7/1209/s1, Table S1. Full table for Table 1, Table S2. Full table for Table 2, Table S3. Full table for Table 3, Table S4. Full table for Table 4.

Author Contributions

Conceptualization, Y.Z., S.Y., and J.L. (Jiancheng Li); methodology, Y.Z. and J.R.; software, Y.Z. and J.L. (Jiesi Luo); validation, Y.Z. and J.L. (Jiesi Luo); formal analysis, Z.Y.; investigation, Y.Z., J.L. (Jiesi Luo), and Y.H.; resources, Z.Y.; data curation, Y.Z., J.L. (Jiesi Luo), and Y.H.; writing-original draft preparation, Y.Z.; writing-review and editing, J.R.; visualization, Y.Z. and J.L. (Jiesi Luo); supervision, J.L. (Jiancheng Li); project administration, S.Y.; funding acquisition, S.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China, grant number 41574017, 41874025, and 41474022, and the Hubei Provincial Natural Science Foundation, grant number 2019CFB555. The APC was funded by the National Natural Science Foundation of China, grant number 41574017.

Acknowledgments

Four anonymous reviewers have provided valuable comments and suggestions, which helped improve the presentations of this study. Insightful discussions with Wei Chen on harmonic analyses of the uplift rates are highly appreciated. We also thank W.R. Peltier, A. Percell, and M. Schumacher for making their GIA models or GPS solutions publicly available.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Spada, G.; Sabadini, R.; Yuen, D.A.; Ricard, Y. Effects on postglacial rebound from the hard rheology in the transition zone. Geophys. J. Int. 1992, 109, 683-700.

2. Anderson, J.B.; Shipp, S.S.; Lowe, A.L.; Smith Wellner, J.; Mosola, A.B. The Antarctic ice sheet during the last glacial maximum and its subsequent retreat history: A review. Quat. Sci. Rev. 2002, 21, 49-70.

3. Peltier, W.R. Global glacial isostasy and the surface of the ice-age Earth: The ICE-5G (VM2) model and GRACE. Ann. Rev. Earth Planet. Sci. 2004, 32, 111-149.

4. Kendall, R.A.; Mitrovica, J.X.; Milne, G.A. On post-glacial sea level-II. Numerical formulation and comparative results on spherically symmetric models. Geophys. J. Int. 2005, 161, 679-706.

5. Klemann, V.; Martinec, Z. Contribution of glacial-isostatic adjustment to the geocenter motion. Tectonophysics 2011, 511, 99-108.

6. Peltier, W.R. The history of Earth's rotation: Impacts of deep Earth physics and surface climate variability. In Treatise on Geophysics, 2nd ed.; Schubert, G., Ed.; Elsevier: Oxford, UK, 2015; Volume 9, pp. 221-279.

7. Stuhne, G.R.; Peltier, W.R. Reconciling the ICE-6G_C reconstruction of glacial chronology with ice sheet dynamics: The cases of Greenland and Antarctica. J. Geophys. Res. Earth Surf. 2015, 120, 1841-1865.

8. Caron, L.; Ivins, E.R.; Larour, E.; Adhikari, S.; Nilsson, J.; Blewitt, G. GIA model statis tics for GRACE hydrology, cryosphere, and ocean science. Geophys. Res. Lett. 2018, 45, 2203-2212.

9. Watkins, M.M.; Wiese, D.N.; Yuan, D.N.; Boening, C.; Landerer, F.W. Improved methods for observing Earth's time variable mass distribution with GRACE using spherical cap mascons. J. Geophys. Res. Solid Earth 2015, 120, 2648-2671.

10. Wiese, D.N.; Landerer, F.W.; Watkins, M.M. Quantifying and reducing leakage errors in the JPL RL05M GRACE mascon solution. Water Resour. Res. 2016, 52, 7490-7502.

11. Save, H.; Bettadpur, S.; Tapley, B.D. High resolution CSR GRACE RL05 mascons. J. Geophys. Res. Solid Earth 2016, 121, 7547-7569.

12. Chen, W.; Luo, J.; Ray, J.; Yu, N.; Li, J.C. Multiple-data-based monthly geopotential model set LDCmgm90. Sci. Data 2019, 6, 288.

13. Argus, D.F.; Peltier, W.R.; Drummond, R.; Moore, A.W. The Antarctica component of postglacial rebound model ICE-6G_C (VM5a) based upon GPS positioning, exposure age dating of ice thicknesses, and relative sea level histories. Geophys. J. Int. 2014, 198, 537-563.

14. Peltier, W.R.; Argus, D.F.; Drummond, R. Space geodesy constrains ice-age terminal deglaciation: The global ICE-6G_C (VM5a) model. J. Geophys. Res. Solid Earth 2015, 120, 450-487.

15. Peltier, W.R.; Argus, D.F.; Drummond, R. Comment on "An Assessment of the ICE-6G_C (VM5a) Glacial Isostatic Adjustment Model" by Purcell et al. J. Geophys. Res. Solid Earth 2018, 123, 2019-2028.

16. Schumacher, M.; King, M.A.; Rougier, J.; Sha, Z.; Khan, S.A.; Bamber, J.L. A new global GPS data set for testing and improving modelled GIA uplift rates. Geophys. J. Int. 2018, 214, 2164-2176.

17. Ray, J.; Griffiths, J.; Collilieux, X.; Rebischung, P. Subseasonal GPS positioning errors. Geophys. Res. Lett. 2013, 40, 5854-5860.

18. Ray, J.; Rebischung, P.; Schmid, R. Dependence of IGS products on the ITRF datum. In Reference Frames for Applications in Geosciences; Altamimi, Z., Collilieux, X., Eds.; Springer: Berlin, Germany, 2013; Volume 138, pp. 63-67.

19. Wahr, J.; Wingham, D.; Bentley, C. A method of combining ICESat and GRACE satellite data to constrain Antarctic mass balance. J. Geophys. Res. 2000, 105, 16279-16294.

20. Wahr, J.; Han, D.Z.; Trupin, A. Predictions of vertical uplift caused by changing polar ice volumes on a viscoelastic earth. Geophys. Res. Lett. 1995, 22, 977-980.

21. Peltier, W.R. The LAGEOS constraint on deep mantle viscosity: Results from a new normal mode method for the inversion of viscoelastic relaxation spectra. J. Geophys. Res. 1985, 90, 9411-9421.

22. Purcell, A.; Dehecq, A.; Tregoning, P.; Lambeck, K.; Potter, E.-K.; McClusky, S. Relationship between glacial isostatic adjustment and gravity perturbations observed by GRACE. Geophys. Res. Lett. 2011, 38, L18305.

23. van der Wal, W.; Wu, P.; Sideris, M.G.; Shum, C.K. Use of GRACE determined secular gravity rates for glacial isostatic adjustment studies in North-America. J. Geodyn. 2008, 46, 144-154.

24. Tregoning, P.; Ramillien, G.; McQueen, H.; Zwartz, D. Glacial isostatic adjustment and nonstationary signals observed by GRACE. J. Geophys. Res. 2009, 114, B06406.

25. Wu, X.; Heflin, M.B.; Schotman, H.; Vermeersen, B.L.; Dong, D.; Gross, R.S.; Owen, S.E. Simultaneous estimation of global present-day water transport and glacial isostatic adjustment. Nat. Geosci. 2010, 3, 642-646.

26. Jia, L.L.; Wang, H.S.; Xiang, L.W. Uncertainty of approximate relationship between GIA induced viscous gravity and radial displacement. Earth Sci. J. China Univ. Geosci. 2014, 39, 905-914.

27. Guo, J.Y.; Huang, Z.W.; Shum, C.K.; van der Wal, W. Comparisons among contemporary glacial isostatic adjustment models. J. Geodyn. 2012, 61, 129-137.

28. Purcell, A.; Tregoning, P.; Dehecq, A. An assessment of the ICE-6G_C (VM5a) glacial isostatic adjustment model. J. Geophys. Res. Solid Earth 2016, 121, 3939-3950.

29. Purcell, A.; Tregoning, P.; Dehecq, A. Reply to comment by W. R. Peltier, D.F. Argus, and R. Drummond on "An assessment of the ICE6G_C (VM5a) glacial isostatic adjustment model". J. Geophys. Res. Solid Earth 2018, 123, 2029-2032.

30. Fretwell, P.; Pritchard, H.D.; Vaughan, D.G.; Bamber, J.L. Bedmap2: Improved ice bed, surface and thickness datasets for Antarctica. Cryosphere 2013, 7, 375-393.

31. Lidberg, M.; Johansson, J.M.; Scherneck, H.-G.; Davis, J.L. An improved and extended GPS-derived 3D velocity field of the glacial isostatic adjustment (GIA) in Fennoscandia. J. Geod. 2007, 81, 213-230.

32. Heiskanen, W.A.; Moritz, H. Physical Geodesy; W. H. Freeman and Company: San Francisco, CA, USA, 1967; pp. 22-32.

33. Wahr, J.; Molenaar, M.; Bryan, F. Time variability of the Earth's gravity field: Hydrological and oceanic effects and their possible detection using GRACE. J. Geophys. Res. 1998, 103, 30205-30229.

34. Han, S.C.; Shum, C.K.; Jekeli, C.; Kuo, C.Y.; Wilson, C.; Seo, K.W. Non-isotropic filtering of GRACE temporal gravity for geophysical signal enhancement. Geophys. J. Int. 2005, 163, 18-25.

35. Wouters, B.; Schrama, E.J.O. Improved accuracy of GRACE gravity solutions through empirical orthogonal function filtering of spherical harmonics. Geophys. Res. Lett. 2007, 34, L23711.

36. Klees, R.; Revtova, E.A.; Gunter, B.C.; Ditmar, P.; Oudman, E.; Winsemius, H.C.; Savenije, H.H.G. The design of an optimal filter for monthly GRACE gravity models. Geophys. J. Int. 2008, 175, 417-432.

37. Davis, J.L.; Tamisie, M.E.; Elósegui, P.; Mitrovica, J.X.; Hill, E.M. A statistical filtering approach for Gravity Recovery and Climate Experiment (GRACE) gravity data. J. Geophys. Res. 2008, 113, B04410.

38. Chen, W.; Li, J.C.; Ray, J.; Shen, W.; Huang, C. Consistent estimates of the dynamical figure parameters of the Earth. J. Geod. 2015, 89, 179-188.