A new global GPS data set for testing and improving modelled GIA uplift rates

,


I N T RO D U C T I O N
Glacial isostatic adjustment (GIA) is the response of the solid Earth to past ice-ocean loading changes.It influences the Earth's shape, gravity field, geoid and axis of rotation via redistribution of mantle material.By consequence it also affects regional sea level patterns.GIA can be forward modelled (e.g.Peltier 2004), inverse modelled (e.g.Riva et al. 2009) and can be observed geodetically (e.g. with GPS; King et al. 2010).Models of GIA are, however, very different regionally (e.g.Martín-Español et al. 2016a) and often poorly constrained by observations.GPS measurements can be used to test various GIA models and the global GPS network now contains many thousands of stations.However, the time-series of vertical land motion (VLM) obtained from GPS stations contain other signals, such as jumps (e.g.due to earthquakes or hardware changes) and longerterm changes due to natural and anthropogenic processes including the elastic rebound of the Earth surface due to present-day ice mass loss, tectonics and local hydrology (e.g.groundwater pumping).In addition, the data are typically noisy including, in many cases, a significant subannual component.
The VLM from permanent GPS records have been used multiple times to help constrain global or regional models of GIA uplift (e.g.Whitehouse et al. 2012;Argus et al. 2014a;Peltier et al. 2015) or to solve directly for GIA using a data-driven approach (Wu et al. 2010).For such applications, special care is required for the interpretation of the VLM in ice-covered regions, such as Antarctica, Greenland and Alaska, as the elastic rebound due to present-day ice load changes can dominate the signal.To separate the elastic 2164 C The Author(s) 2018.Published by Oxford University Press on behalf of The Royal Astronomical Society.
Downloaded from https://academic.oup.com/gji/article-abstract/214/3/2164/5048688 by DTU Library -Technical Information Center of Denmark user on 13 July 2018 signal from the viscous GIA signal, requires high resolution spatial loading data (on the order of 20 km) as elastic deformation is highly localized (e.g.Spada et al. 2012).Several studies have noted significant differences between global and regional GIA forward model solutions and GPS VLM for Antarctica (e.g.Martín-Español et al. 2016a), for Greenland (e.g.Khan et al. 2016) or applied to tide gauge data (Wöppelmann et al. 2009;King et al. 2012).
The aim of this study is to produce a global vertical velocity field due to GIA that can be used either to test GIA models or, as intended here, to directly determine GIA from observations.Such a GPS data set can be used to perform a Bayesian update on GIA models, as performed in Martín-Español et al. (2016a) for Antarctica, but on a global scale.By combining data and prior information from models, their discrepancies can be reduced, and thus, geophysical processes of the Earth system can be better represented.
To achieve this, we must remove both noise and signal due to other geophysical processes as outlined above.We explain the steps and methodology we use to achieve this and present a comparison with an ensemble of 13 global GIA estimates.Here, we use the GPS data set of the Nevada Geodetic Laboratory (NGL) as the starting point for providing an observational estimate of global GIA VLM.A novel fully-automatic post-processing strategy is developed to deal with the challenges of GPS time-series analysis in general, and for GIA purposes in particular, including outlier and jump detection, atmospheric mass loading correction, elastic signal correction and filtering for stations where other sources of VLM are likely to dominate over GIA.In order to accurately account for the elastic response of the Earth's crust over Antarctica and Greenland, separate data sets are used that have already been corrected for the contemporary ice mass loading impact on elastic deformation using high-resolution ice mass balance time-series (Khan et al. 2016;Martín-Español et al. 2016b).We compare our novel global GPS data set, denoted as GlobalMass (GM, after the project title) GPS data set in this paper, with 13 global GIA solutions that have been previously compared (Guo et al. 2012).

Global GPS data set from NGL
The GPS data provided by the Nevada Geodetic Laboratory (NGL) at the University of Nevada, Reno (UNR) were used in this study.The selected data set is provided as north, east and up components for more than 15 700 GPS sites in the IGS08 reference framework, with its origin in the centre of mass of the total Earth system (CM).The locations of the GPS sites are available, as well as a database of jumps occuring in the GPS time-series, though outliers and jumps are not removed from the time-series.In addition, no correction due to atmospheric mass loading is applied.Details on the data set and all applied conventions are documented on the NGL webpage (http://geodesy.unr.edu/).-Español et al. (2016b) use the A-NET GPS data in conjunction with additional data sources, such as GRACE and satellite altimetry, to solve for the mass balance of the Antarctic ice sheet and estimate the regional GIA.The GPS data are provided in the ITRF2008 reference frame, which is effectively identical to IGS08 (Rebischung et al. 2012).Corrections due to atmospheric mass loading and the solid Earth elastic response to present-day surface mass changes are already applied.For the latter, ice mass trends are derived from a rigorous statistical combination of remotely sensed gravity, altimetry and GPS observations within a Bayesian Hierarchical Model (BHM, Zammit-Mangion et al. 2015;Martín-Español et al. 2016b).In this study, corrected uplift rates at 65 GPS sites between 2009 and 2013 are used.

Regional GPS data set for Greenland
For Greenland, we also use the G-NET GPS data set of 54 sites corrected for elastic VLM.The starting points of the GPS timeseries vary between 1995 and 2010 and data are considered until 2015.The data are processed in the IGS08 reference frame.Details of the post-processing of GPS data and the estimation of the elastic correction for Greenland are provided in Khan et al. (2016).

Global GIA forward model solutions
A variety of global GIA forward model solutions are used in geodetic and geophysical studies to address changes in the surface mass balance, sea level, and solid Earth.These models differ in their two main assumptions about the deglaciation history and the viscoelastic solid Earth structure and rheology.Different combinations and model parameter assumptions lead to different spatial fields of GIA.In Guo et al. (2012), 13 GIA forward model solutions and one data-driven solution have been compared.The authors provided eleven of the GIA solutions (Pel-4-VM2 (with ice loading history Guo et al. 2012, Table 1 for details), and two additional solutions (Pel-6-VM5 (IH=ICE-6G) and Pur-6-VM5 (IH=ICE-6G) that both use the ICE-6G C model; details may be found in Peltier et al. 2015;Purcell et al. 2016).Only the ICE-6G model is tuned to fit GPS; the other ice histories are independent.It is important to note that the GPS data are not used in the models, but they are used to tune the models to fit the GPS data, the relative sea level curves, and other data.In total, we compare a set of 13 GIA forward model solutions with the GPS data set in this study.

Post-processing of the GPS time-series
In this work, we are interested in post-processing time-series at permanent GPS stations globally for GIA-related applications.Thus, we are only interested in the VLM and the post-processing steps described in this section are only applied to the vertical (or 'up') component of the provided GPS time-series.The horizontal components are not considered.The various data sets and post-processing steps are summarized in Fig. 1.In this section, data handling and corrections to the NGL time-series are described.In Section 3.2, the corrections applied to the calculated GPS uplift rates (i.e.velocities) for GIA purposes are presented.

Period considered for trend estimation
The estimation of linear trends from the elastic-corrected (see Section 3.2.2) vertical component of GPS data should be identical over arbitrary periods, since the GIA signal is assumed to generate a  constant velocity over the time periods considered here.To investigate this, we consider two periods, 2005-2015 (many GPS sites do not cover the full period) and the entire period for which GPS data are recorded (which varies considerably among the stations), and analyse the statistical significance of the difference between the trend estimates.The result shows that for the vast majority of the stations (98.9 per cent out of 3101 tested sites), considering a shorter or longer data period has no significant influence on the estimated VLM, that is the difference in trends lies within three times the uncertainties of the trend estimate with the larger uncertainty.Thus, we can conclude that for the majority of sites, the time period does not significantly influence the final GPS data set.Due to the effect that a changing climate over several decades might have on the trend estimation (for example, possibly from unresolved trends in atmospheric loading, i.e. signal that is not captured by the reanalysis product), we selected the period 2005-2015 for this study.

Excluding stations with short data record
Stations that have time-series with less than four years of daily data (based on the total number of data points, i.e. a threshold of 4 × 365 was used) are not considered in the post-processing of the global GPS data set.This value is based on previously published literature stating that at least four years of data are required to mitigate the influence of interannual deformations on the linear trend estimate (e.g.Santamaría-Gómez & Mémin 2015).Additional quality control tests were also applied.If a limited number of data points (less than half a year) exist with large gaps to the remaining data points in the time-series, these data points have been removed.In addition, if few data points (less than 2 weeks) exist between two jumps in the time-series, these data points have been removed, and only one jump is estimated.Excluding these data points helps to mitigate unwanted effects during the next post-processing steps, for example, unstable linear trend estimations.

Outlier detection
Outliers (i.e.data points that are far outside of the expected range of values) frequently occur in GPS time-series, for example, due to sudden icing and de-icing of antennas (spanning days to months) or antennas being buried by snow.Thus, a fully automatic outlier detection algorithm is applied to each 'jump period' separately.In this paper, the term jump period refers to the period between two consecutive reported jumps of a time-series.First, a linear fit considering only the intercept a 0 and linear trend a 1 is estimated for the vertical component of each GPS time-series y at given time t y = a 0 + a 1 t + r (1) and the residuals r of each data point with respect to the linear fit are calculated.A 99 per cent confidence interval is determined based on the root mean square (RMS) of these residuals.All data points that lie outside of the 99 per cent confidence interval are defined as outliers and are removed from the time-series.The selected confidence interval seems reasonable.A quality control (see Section3.2.1) showed that only at station name KSMV, a 95 per cent confidence interval would be more appropriate to remove all occurring outliers.The confidence interval for this station has been modified in the algorithm.
Additionally considering an annual signal represented by a sine and cosine term with parameters a 2 and a 3 for the outlier detection, that is y = a 0 + a 1 t + a 2 sin(2π t) + a 3 cos(2π t) + r, does not considerably influence the estimate except at one GPS site.At station PBRI, adding the annual signal results in a few outliers that are not detected and therefore bias the trend estimate.Thus, a simple linear fit seems to be better suited than including the addition of an annual Downloaded from https://academic.oup.com/gji/article-abstract/214/3/2164/5048688 by DTU Library -Technical Information Center of Denmark user on 13 July 2018 signal.This is likely because the stations selected do not contain a strong seasonal hydrological signal as they are not close to major aquifers or catchments.

Jump detection
The NGL group invested considerable effort in providing a database that reports jumps in the GPS time-series due to hardware issues (e.g.switching the GPS antenna) and due to geophysical signal (e.g.earthquakes).However, for the purposes of our work (GIA assessment) we found that, for a small number of stations, jumps were omitted from the NGL jump database.To detect these unreported jumps, we have designed and implemented an automatic procedure comprising the following steps (Fig. 2): (1) Calculating residual time-series for each station.
(2) Applying a moving average filter to the residual time-series.
(3) Computing the differences of the successive moving average values.
(4) Locating groups of jumps with large differences in the moving average values.
(5) Determining the accurate position of the detected jumps.
(1) To calculate the residual time-series r, a linear adjustment is applied considering intercept a 0 , linear trend a 1 , annual signal with parameters a 2 and a 3 , as well as the biases The matrix B contains the information on the location of the biases.For example, if jump j occurs at a time between t k and t k + 1 , the entries in row B j will contain zeros prior to the jump occurrence, that is in columns one to k, and ones after the jump occurrence, that is in columns k + 1 to the last column K (is known as the Heaviside function) The adjusted time-series is subtracted from the original time-series to calculate the residual time-series.(2) An unreported jump will significantly change the average of the residuals within a predefined window.The window length is set to 7 days, which implies that only one unreported jump occurs within 1 week.We consider this reasonable, as the vast majority of jumps are already reported in the NGL database.
(3) The differences between the 7-day moving average values provide information on the location of a jump.Thereby, one jump influences seven successively computed average values.
(4) To locate these groups of jumps, a threshold of 3σ y was selected for the chosen window length, with σ y representing the temporal variability of the complete available time-series y (often called the standard deviation of the time-series).Due to the high noise level within the GPS time-series, a minimum threshold of 3.5 mm was empirically defined.In case that the 3σ y value is smaller than the threshold for a time-series, that is that the time-series has a small temporal variability, it is replaced by 2σ y , which allowed for a better detection of unreported jumps.The points of time at which groups of jumps occur provide information on the number of unreported jumps and can be used to separate these groups.( 5) The maximum difference of the observed vertical component gives the exact location of each jump.Finally, the detected jumps are added to the NGL jump database file and this is used for all following computations.

Atmospheric mass loading
It has been shown in previous studies that deformations due to nontidal atmospheric, oceanic and terrestrial water mass loading have a significant impact on geodetic positioning time-series (e.g.Dach et al. 2011;Fritsche et al. 2012;van Dam et al. 2012;Santamaría-Gómez & Mémin 2015).The atmospheric mass loading has a spatial wavelength of about 1000 km and effectively adds noise and possibly small biases to the vertical linear trends derived from the NGL GPS time-series.Since atmospheric mass loading is well modelled, we use the data provided by the International Mass Loading Service to compute the loading at the GPS site locations considered in this study (http://massloading.net/).We selected the atmospheric re-analysis product MERRA2 (available from 1980 onwards) and download the product as global 2 × 2 maps with a temporal resolution of 6 hours.The time-series need to be in a reference frame with the origin in the centre of the solid Earth system mass to be consistent with the GPS data set at non-secular periods (Dong et al. 2003), which is realized by also downloading gridded global maps representing the degree-1 terms.

Trend and bias estimation
An automatic trend and bias estimation has been implemented, in which the information from the extended jump database (see Section 3.1.4)is used for bias estimation at each jump location.A linear adjustment is performed for each time-series, with intercept a 0 , linear trend a 1 , cyclic patterns (represented by a sine and cosine term) with parameters a 2 and a 3 for the annual and parameters a 4 and a 5 for the semi-annual signal, as well as biases b from the extended jump database (see Section 3.1.4)considered as parameters (similar to, e.g.Roggero 2012) y = a 0 + a 1 t + a 2 sin(2π t) + a 3 cos(2π t) +a 4 sin(4π t) + a 5 cos(4π t) Matrix B is defined in eq. ( 3).The GPS time-series exhibit temporally correlated (non-Gaussian) noise, which persists in the residual term after ordinary least squares fitting, as identified using the empirical autocorrelation function.The ordinary least squares fitting uncertainties for coefficient estimators are derived on the basis that the residual process is independent and identically distributed, and are unreliable if the residual is autocorrelated.Often, power law noise has been estimated to represent the noise level more realistically (e.g.Bos et al. 2013).Alternatively, Khan et al. (2016) used 30-day averages of the daily vertical solutions to consider temporally correlated noise.The RMS of the 30 daily values to the 30-day average were calculated to represent uncertainties of monthly values and these were used to propagate uncertainties.This approach requires jumps to be removed prior to the estimation of the 30-day average to avoid introducting biases.Therefore, the uncertainty of the jump detection is not included in the final estimation of trend uncertainties.
In this paper, we follow Khan et al. (2016), but rather than fitting a specific parametric form to the autocorrelation, we adopted the simple and robust alternative of thinning the time-series, to the point where the correlation between consecutive values was small enough to approximate the independent and identically distributed condition.Examination of some of the residual autocorrelation functions suggested that thinning to a time-step of 15 days would be adequate, which we confirmed by refitting the regressions after thinning, and recomputing the residual autocorrelation functions.We adopted thinning to a time-step of 15 days as a simple rule for all  GPS time-series.The final linear trend uncertainties are therefore somewhat larger (but likely more realistic) than those following the approach in Khan et al. (2016).

Post-processing of GPS uplift rates for GIA purposes
The generic post-processing steps explained in Section 3.1 have to be performed for any GPS time-series analysis and are independent of the application.In this section, we focus on the specific postprocessing strategy we have developed to generate a global GPS data set suitable for GIA applications.The strategy accounts for different physical processes, for example hydrological mass loading, tectonics and earthquakes (Sections 3.2.1 and 3.2.3),as well as elastic deformation of the Earth crust (Section 3.2.2) and the change of the rotational pole (Section 3.2.4).We also explain how corrections are applied to separate the GIA signal from other signals, and how GPS sites where VLM does not primarily reflect the GIA signal are excluded from our data set.Some of our corrections were not available as time-series, nor were some of the GPS data.Consequently, when applying corrections we used rates throughout and in some cases this could introduce biases due to time-variability of the underlying signal being corrected.

Station selection based on prior information from GIA forward models
A first selection of NGL GPS sites is performed based on prior information from a set of 13 global GIA forward model solutions to ensure that only stations that primarily represent the GIA signal are considered in the global data set.A mean value of the 13 GIA models and the model spread are calculated for each 1 • × 1 • grid cell globally.Then, the globe is divided into seven regions based on the signal strength of GIA and other geophysical signals, for example elastic rebound of the Earth, tectonics or hydrology (Fig. 3a): (1) Northern Europe and Asia, (2) North America, (3) mid latitudes, (4) Antarctica, (5) Greenland, (6) Gulf of Alaska, (7) Iceland and (8) Hawaii.The observed uplift rates in some parts of regions (1) and (2) are dominated by the GIA signal (e.g.Scandinavia and Canada), while region (3) is marginally influenced by GIA.In regions (4) and ( 5) a mixture of GIA signal and elastic response to present-day mass loss is observed, which have to be appropriately accounted for.In addition, in regions ( 6) and ( 7), a mixture of viscous and elastic response of the Earth's crust, as well as tectonic signals are measured, which are difficult to separate.Regions (4) to ( 7) are defined based on the Randolph Glacier Inventory (version 6.0) order 1 polygons (RGI Consortium 2017).
For each of the regions ( 1)-( 3), the minimum and maximum values of the mean GIA uplift rates, that is the mean of the 13 GIA forward models, are calculated and three times the model spread is either subtracted or added to the values to produce conservative estimates of the range of realistic GIA-related uplift values.We selected three times the empirical standard deviation (i.e.model spread), which corresponds to a 99.7 per cent confidence level, and thus, we ensure a pessimistic estimate of the range limits and do not exclude stations that actually contain a GIA signal.The ranges are, for region (1), between −6 and 16 mm yr −1 , region (2) between −10 and 24 mm yr −1 and region (3) between −4 and 4 mm yr −1 .These limits are used to define which of the GPS sites are selected for the global data set for GIA purposes.GPS sites that exhibit much larger trends are excluded, that is 12.9 per cent of the GPS sites (Fig. 3b).It is evident that for some regions a cluster of GPS sites is excluded, for example in Japan (Japan, Izu Ogasawara and Ryuku trench), Malaysia (Java (Sunda) trench), New Zealand (Kermadic trench), Italy (fault line between European and African plates), Chile (Peru-Chile trench) and California (San Andreas fault) which are assumed to be associated with tectonic VLM.In addition, the strong VLM trends in California might be associated with tectonics (Kreemer et al. 2010) and the extreme hydrological drought between 2012 and 2016 (Argus et al. 2014b;Borsa et al. 2014;He & Gautam 2016).In this study, the data for regions (4) and ( 5), as explained earlier, are replaced by external data sets, and data for regions (6), ( 7) and ( 8) are excluded from the data set (see Section 3.2.2 for details).

Elastic correction for ice-covered regions and far field
Special care is required for the interpretation of the uplift rate in ice covered regions as the elastic rebound due to present-day ice load changes can dominate the signal.High-resolution ice mass balance time-series were used to estimate the contemporary ice mass loading impact on elastic deformation.A correction in terms of uplift rates was applied to the GPS uplift rates in Antarctica and  Greenland, in order to accurately account for the elastic response of the Earth's crust in these regions.For Antarctica, we use the elastic corrected data set described in Section 2.2 with 65 GPS sites, and for Greenland, the data set described in Section 2.3 with 54 GPS sites.It is clearly visible in Fig. 4 that neglecting the elastic response would significantly contaminate, and in some areas dominate, the estimated GIA uplift rates.In addition to Antarctica and Greenland, it might be necessary to correct for the elastic rebound of the Earth's crust also in other glaciated regions, for example the Gulf of Alaska, Iceland Svalbard and Canadian Arctic.Areas underlain with low viscosity mantle, such as the Gulf of Alaska, Iceland and Patagonia, no longer experience substantial GIA signal in response to the last glacial maximum mass loss (Dietrich et al. 2010;Sato et al. 2011;Auriac et al. 2013).However, the viscous response of the past century since the Little Ice Age is visible in the data (Sato et al. 2011).The GPS data include a mixed signal of this viscous response, of elastic response to present-day mass losses and of tectonics.In Alaska, the latter can constitute a significant part of the signal (Fu & Freymueller 2013).Thus, GPS sites in this region are not included in our final global data set (138 GPS sites excluded).Similar issues can be observed for Iceland, and the corresponding 7 GPS stations are excluded.It is important to note that GRACE detects the viscous response of the Earth's crust for these regions, which will not be reflected in our global GPS data set, since it is currently not possible to accurately separate the different physical sources observed at these GPS sites in the Gulf of Alaska and Iceland.Riva et al. (2017) reported that the elastic response in Antarctica and Greenland has a long wavelength influence, that is influencing VLM thousands of kilometres from the poles.They used yearly mass losses from glaciers, ice caps, Greenland and Antarctic ice sheets between 1902 and 2014 to determine the solid Earth elastic response to surface load changes locally and for the far field.The sea level equation was solved in a centre of mass of the Earth system (CM) frame for each load and each year, meaning that the elastic response to ocean loading is also accounted for in this data set.The sum of all loads defines the total elastic response and was used to estimate linear trends over different time windows.The authors found that maximum uplift rates in the far field increased to about 1.0 mm yr −1 in the last decade affecting especially the Northern Hemisphere.The data have been published (https://data.4tu.nl/repository/uuid:fb667e7a-52f 3-4876-8cab-ae7a2ddaf0db; last access: 2017 October 02) and used in this study to correct for the elastic far field effect of present-day mass losses in regions (1), ( 2) and (3).It is important to note that this correction includes the effect of changes in the rotational pole (see Section 3.2.4 for details).Based on the grids of Riva et al. (2017) spanning 2005-2014, we computed uplift rates at the around 4000 considered NGL GPS sites and show them in Fig. 5 (no data are available for 2015 and it is assumed that the elastic uplift rate does not change during this year).No uncertainties are provided for this correction and are therefore not reflected in the uncertainty estimate of our GPS data set.

Median filter for removal of local and non GIA signals
The non-tidal oceanic and hydrological loading have a similar effect to the atmospheric loading on the GPS time-series but both are less well modelled in general (Santamaría-Gómez & Mémin 2015), which means that the loading computations are not as accurate.For example, it would be possible to subtract the hydrological loading by removing simulations from global hydrological models with daily time steps as was done by, for example in Simon et al. (2017) for North America.However, numerous studies showed that these models are not in close agreement for many regions worldwide.The representation of long-term trends, in particular, is a common problem among hydrological model simulations (Döll et al. 2014), which is the component in which we are primarily interested.Therefore, a hydrological mass loading correction would possibly introduce even more uncertainties to our analysis.In addition, a spatial resolution mismatch occurs between GPS observations (point-wise) and models (grid-based, usually provided on 0.5 • or 1 • global grids) which would further increase uncertainties.Thus, no explicit correction is applied for hydrological mass loading in this study but instead we perform a spatial filtering strategy to select stations that are predominantly influenced by the long wavelength GIA signal and to exclude stations that are affected by local to regional hydrology (such as groundwater pumping).This is consistent with Santamaría-Gómez & Mémin (2015), who reported that 'unless hydrological loading models are improved, there is presently no robust solution to mitigate these uncertainties other than to increase the time-series lengths'.
To ensure that stations influenced by local effects, such as land hydrology and earthquakes (not already filtered out) are removed, we used a low-pass spatial filtering approach.For this, a median filter with a radius of 500 km was applied.We also tested various filters between 250 and 1000 km radius.We found that using a larger radius smoothed the signal too much while, in contrast, using a smaller radius resulted in too few stations available to calculate the median value.Therefore, 500 km was chosen as a compromise, which roughly represents the correlation lengths of the GIA signal in the forward models.The value at a particular NGL GPS site is assigned the median value within the predefined radius.The difference between the median value and the original value is used as a selection criterion to either keep or remove the station from the global GPS data set.First, stations with an opposite sign of the linear trend before and after the application of the median filter are removed.Subsequently, stations with linear trends that strongly differ before and after the application of the median filter are removed as well.For this, a threshold of three times the uncertainties of the original linear trend was used.This results in 40.1 per cent of the NGL GPS sites being excluded.

Changes in rotational pole
A rapid shift of the Earth's rotational pole towards the east is observed since around 2005, primarily driven by ice sheet melting in Antarctica and Greenland (Chen et al. 2013).King & Watson (2014) noted out that the International Earth Rotation Service's (IERS) elastic pole tide model does not correctly account for this effect.This model is usually used to correct short-term polar motion in geodetic analyses.Hence, the GPS uplift rates determined in this study will be influenced by the deformation resulting from the deviations in polar motion from its longer-term path with a maximum (minimum) of up to 0.25 mm yr −1 uplift (subsidence) over the U.S. Pacific Coast and South Africa (Europe and south Pacific islands; Fig. 6).We used the data set reported in King & Watson (2014) to estimate the effect of polar motion at the considered GPS sites worldwide.The yearly deformations from 1980 to 2015 due to polar motion are used to estimate linear trends over the data period considered at each site, that is 2005-2015 for the NGL sites, 2009-2013 for the sites in Antarctica, and min./max.1995/2010 to 2015 depending on the sites in Greenland.These values are subtracted from the GPS uplift rates of the external data sets for Antarctica and Greenland (but not for the NGL sites since the effect is already included in the correction for the far field elastic response described in Section 3.2.2).

Assessment of GIA forward model solutions
The GPS data are provided in a reference frame with the origin in an approximation to the centre of mass of the total Earth's system (CM; comprising solid Earth and fluid components, such as atmosphere, oceans and land hydrology) on long-term scales (Dong et al. 2003).In contrast, the GIA forward models provide information in a reference frame with the origin in their own realization of the   centre of mass of the solid Earth (CE).The two may be transformed by a translation of the geocentre (e.g.Argus et al. 2014a).A consistent assessment of the GIA forward model solutions using the global GPS data set is therefore only possible if the origin shift is estimated and removed or a model is applied (Argus & Peltier 2010;King et al. 2012).Thus, we estimated the origin translation between our developed global GPS data set and each of the 13 GIA forward model solutions.The spatial distribution of the GPS sites is far from equal, and this might bias the estimation of the geocentre motion rate.Therefore, we follow the approach described in (King et al. 2012, Text S1 of the auxiliary material) using 20 • × 10 • windows (the authors proposed 30 • × 10 • windows in their paper), i.e. 18 windows in longitudinal and latitudinal direction, and calculated the velocity block median values from the GPS sites contained in each window.These are then converted from the local topocentric systems (setting the velocities of the north and east component to zero) to geocentric cartesian coordinates.The same is repeated for the uplift rates of all provided GIA models but only at the locations of the GPS sites.Linear trends of the CM and CE geocentre in the geocentric coordinate system are derived by averaging over all velocity block median values of the GPS data and the GIA model values, respectively.The difference of the two geocentres is calculated each time, that is CE-CM, and defines the geocentre motion (Table 1).It is obvious that the geocentre motion estimate is sensitive to the slection of the GIA forward model.The estimate is converted back to north, east and up components in local topocentric systems at all GPS sites.Subsequently, the derived corrections of the up component are removed from the GPS uplift rates, that is up C E G P S = up C M G P S − up C E−C M , to transform the data set into a reference frame with the origin in the CE to be consistent with the GIA forward model solutions (see Fig. 7 for geocentre motion with respect to the ICE-6G model).We did not estimate uncertainties for this correction, since no uncertainty information is provided for GIA forward model solutions.Considering an unceratinty estimate would nonetheless have a very small influence on the GPS uplift rate uncertainties.In summary, the effect of the frame origin transformation on the GPS uplift rates is very small (less than ±0.2 mm yr −1 for most GIA models as also found in King et al. 2012).

The novel GM GPS data set
The final post-processed global GM GPS data set is provided in a reference frame with its origin in the centre of mass of the solid Earth (CE).The estimated GIA uplift rates are shown in Fig. 8a).It evident that the GIA uplift rates are large in regions that are well known to have been covered by ice during the Last Glacial Maximum approximately 22 000 years ago.In Antarctica, the highest GIA uplift rates occur in the west with rates of up to around 20 mm yr −1 , and in Greenland the largest rates of 12 mm yr −1 occur close to the Iceland Hot Spot track near Kulusuk at the southeastern coast (Khan et al. 2016).Similarly, in Scandinavia and Canada, maximum rates of up to around 14 mm yr −1 are observed.In contrast, at the GPS sites located in lower latitude bands only small trends are present, most of which are statistically insignificant based on a 99 per cent confidence level (shown by hollow circles in Fig. 8a).The associated uncertainties of the GIA uplift rates are shown in Fig. 8b).The smallest uncertainties of mostly less than 0.5 mm yr −1 are estimated for the majority of GPS sites in North America, Europe, many parts of Australia and South Africa.This might be due to the fact that the data records in these regions are usually of very good quality and due to the relatively small trends.Higher uncertainties are obtained for stations along the western coast of South America, the majority of stations in Malaysia and Japan, and for some stations in New Zealand.In these regions, frequently reoccurring earthquakes might influence the quality of the GPS data records, for example many small jumps might be recorded that are not entirely detected and removed from the time-series in the post-processing.In addition, a number of stations in Antarctica, Greenland and Scandinavia show larger uncertainties resulting from the strong GIA uplift rates.In summary, the data set shows a clean GIA signal at all post-processed stations, and is, therefore, suitable to investigate the behaviour of global GIA forward models.

Comparison to 13 GIA forward model solutions
We assess the ICE-6G forward model solution against the GM GPS data set, acknowledging that the latter are not necessarily a perfect representation of present-day VLM due to GIA.We believe, nonetheless, that they represent the best available data set to test the veracity of GIA forward models.The values in the five regions shown in Fig. 9a indicate the spatial RMS difference between the GPS data and the ICE-6G values at all GPS sites (GPS minus ICE-6G).The best agreement is identified for the mid latitudes with an RMS difference of 1.36 mm yr −1 , followed by Northern Europe and Asia, as well as North America with 1.71 and 1.94 mm yr −1 , respectively.The highest discrepancy occurs in Greenland and Antarctica with RMS differences of 3.70 and 3.84 mm yr −1 , respectively.This may be due to the use of a relatively low resolution ice mass timeseries for GPS elastic correction in the ICE-6G set up (Argus et al. 2014a).In addition, poorly constrained ice loading since the Last Glacial Maximum, as well as uncertainties in lower mantle rheology, contribute to the differences between GPS data and global GIA forward models.
The differences at the GPS sites show a clear pattern over Scandinavia that suggest that the ICE-6G model underestimates GIA uplift rates.In North America, the uplift in the north and the subsidence in the centre of the continent towards the south are both too small.At mid latitudes, several spatial patterns are visible.Along the east coast of Australia, New Zealand and Southern America, the model overestimates the GIA uplift rates, which was also reported in King et al. (2012) for Australia and New Zealand when considering the ICE-5G model but positive differences were found for South America.GPS rates in these regions are negative (Fig. 8), so it might be that the ICE-6G model underestimates the magnitude of the GIA signal or that the signal recorded by the GPS is not related to GIA.In contrast, in Southern Africa, Japan, and the west coast of North America, the model underestimates the uplift rates (also reported in King et al. 2012,for coast lines and shown here over inland regions also).In addition, it is apparent that in the centre of the South American continent positive differences occur in coastal regions, but negative differences appear inland.The ICE-6G model systematically overestimates the GIA uplift rates in many parts of Antarctica.Poor agreement between the observations and model is found in coastal Greenland.A systematic pattern can be identified with underestimated GIA uplift rates especially along the east and west coast of Greenland.The largest differences occur around the Iceland Hot Spot track near Kulusuk and near Kangerdlugssuaq (KUAQ; see Khan et al. 2016).In summary, GIA rates in Scandinavia, Canada and western part of North America are underestimated by the ICE-6G model, as well as along the east and west coast of Greenland, while especially in Antarctica rates are overestimated (see also Fig. 9b).
In addition to the detailed analysis of the ICE-6G model, we assess the agreement between the GPS rates and the mean of the 13 GIA forward model solutions (Figs 9c and d) to identify regions of systematic discrepancy between the observations and models.The spatial pattern of the differences is similar to the one shown in Fig. 9a.However, there are less negative values in Greenland, which means that the model mean tends to underestimate the observed GPS rates while ICE-6G tends to overestimate them.It is evident that the model mean systematically underestimates the GPS rates in Antarctica.The scatter plots in Figs9b and 9d show that the GM GPS data is systematically smaller than the 13 GIA models, except for Northern Europe/Asia.One explanation for this difference is that something is systematically missing from the 13 GIA models    such as, for example, that none of the GIA models take into account ice loss since the end of the little ice age (LIA; 1900).Kjeldsen et al. (2015) showed ice loss since LIA is equivalent to 25 mm of global mean sea level rise, with largest loss in southeast and northwest Greenland, coincident with the largest difference in Figs 9a and 9c.
We subsequently assessed 13 global GIA forward models against the GM GPS data set, and summarize the results in Table 2.It is evident that no model performs better than any other model in all regions.Globally, the Pur-6-VM5 model shows the best overall agreement (RMS difference is 1.66 mm yr −1 ), followed by the ICE-6G model (Pel-6-VM5 in Table 2) with a RMS difference of 1.67 mm yr −1 .The largest RMS differences are found for Antarctica and Greenland for all analysed models where they are on the order of two to three times greater than differences for other continental areas.Since the RMS difference does not provide insights into whether the GIA models over-or underestimate the magnitude of the observed rates, we added the difference of the GPS rates and the GIA model rates averaged over the entire globe and the five regions indicated in Fig. 3 (Table 2).This represents a bias between the observations and the model, and shows whether the GIA model over-or underestimates the (spatial) mean GIA rates for the different regions.The bias should be close to zero if model and observations agree well, at least at the GPS locations.For all regions except Antarctica, the models seem to generally overestimate the mean GIA signal, while in Antarctica 12 out of 13 models underestimate the mean GIA signal.We also analysed the mean of the 13 GIA models, which statistically should be always closer to the observations (if the number of models is large enough and they are independent, which is only partially the case here due to common ice history models).In terms of RMS differences, this is true for all regions.However, in terms of the bias, this does not hold for Antarctica, which highlights the large discrepancies between the 13 forward model solutions and the need for additional, for example, observational, information.

C O N C L U S I O N S A N D O U T L O O K
https://doi.pangaea.de/10.1594/PANGAEA.889923In this paper, we developed a novel global GPS data set for GIA applications: the GlobalMass (GM) GPS data set.We developed and demonstrated a fully-automatic post-processing method for GPS time-series, using the NGL GPS database of more than 15 000 sites as a starting point.We applied various corrections for both geophysical and instrumental artefacts, e.g.due to the elastic rebound of the solid Earth and changes of the rotational pole, as well as GPS site selection strategies, to obtain a GPS data set that can be used to characterize the GIA signal globally.The final GM GPS data set provides a clean GIA signal for more than 4000 sites and is, therefore, applicable to both assess global GIA forward models and to improve them.Next, we compared the GM GPS data set with 13 GIA forward models.Significant discrepancies were found, especially for Antarctica and Greenland resulting from uncertain mantle rheology and (recent) ice loading history, as well as the use of spatially coarse resolution present-day ice mass balance estimates to correct for elastic VLM.It should also be mentioned that mantle rheology is thought to be strongly heterogeneous across Antarctica (e.g.van der Wal et al. 2015;Hay et al. 2017;Nield et al. 2018).Although the Pur-6-VM5 and ICE-6G models showed the best agreement with the GM GPS data set globally, different models perform better than others depending on the region considered.The GM GPS data set will Downloaded from https://academic.oup.com/gji/article-abstract/214/3/2164/5048688 by DTU Library -Technical Information Center of Denmark user on 13 July 2018 Table 2. Assessment of the 13 global GIA forward models and their mean in terms of spatial RMS differences in mm yr −1 and in terms of bias (i.e. the difference of GPS rates and the GIA model rates) averaged over the entire globe and regions 1-5 (Fig. 3), that is GPS minus forward model.The best agreement with the novel GM GPS data set is shown by underline, while the worst performance is shown in bold.The model Pel-6-VM5 is identical to the ICE-6G model.

Figure 1 .
Figure 1.Overview of considered data sets and implemented post-processing steps to generate a global GPS data set for GIA purposes.

Figure 2 .
Figure 2. Illustration of the automatic algorithm for jump detection as described in Section 3.1.4.The numbers indicate specific steps described in the text.The vertical lines indicate the position of jumps: red lines show jumps reported in the NGL jump database; and the magenta line shows an unreported jump.

Figure 3 .
Figure 3. (a) Mean field of GIA uplift rates from the 13 GIA forward model solutions, and definition of eight regions based on their response to GIA and other geophysical signals: (1) Northern Europe and Asia, (2) North America, (3) mid latitudes, (4) Antarctica, (5) Greenland, (6) Gulf of Alaska, (7) Iceland and (8) Hawaii.(b) Stations that are excluded from the data set for GIA purposes after selecting GPS sites based on prior information from the GIA forward model ensemble.GPS sites in region (1) are marked by squares, in region (2) by crosses and in region (3) by circles.

Figure 4 .
Figure 4. Data set for Greenland: (a) The original GPS uplift rates were corrected due to (b) the elastic response of the solid Earth due to present-day mass losses resulting in the (c) GIA uplift rates.

Figure 6 .
Figure 6.Correction for the change in rotational pole in mm yr −1 at all considered GPS stations worldwide.The uplift rates due to change in rotational plate are estimated between 2005-2015 for the considered NGL sites, between 2009-2013 for sites in Antarctica, and between 1995/2010-2015 for sites in Greenland (depending on considered data record at each site).The correction is only applied to GPS sites in Antarctica and Greenland, since the effect is already included in the correction for the far field elastic response at the NGL sites.

Figure 7 .
Figure7.Effect of the difference between the centre of mass of the solid Earth (CE) and the centre of mass of the total Earth's system (CM), that is CE-CM, on the GPS uplift rates.This correction is subtracted from the estimated GPS uplift rates in the CM frame to obtain the values in the CE frame.

Figure 8 .
Figure 8.The global GlobalMass (GM) GPS data set after removal of elastic effects and median filtering at 4072 GPS sites: (a) estimated GIA uplift rates in mm yr −1 and (b) their uncertainties.The statistically significant rates in (a) are shown by solid (filled) circles, while not statistically significant rates are shown with hollow circles based on a 99 per cent confidence level.

Figure 9 .
Figure 9. Assessment of the ICE-6G model and the mean of the 13 GIA forward models using the novel GM GPS data set: (a) Differences between GPS data and the simulated values at 4072 considered GPS sites, that is GPS minus ICE-6G.Numbers indicate the spatial RMS differences averaged over the five regions in mm yr −1 .Data in the Gulf of Alaska, Iceland and Hawaii have been excluded.(b) Scatter plot of the GPS data versus the simulated values.Colours in (b) indicate the five regions in (a).The coloured lines indicate the best linear fit for the five regions.(c) Same as (a) but using the mean of the 13 GIA forward models.(d) Same as (b) but using the mean of the GIA models.

Table 1 .
Downloaded from https://academic.oup.com/gji/article-abstract/214/3/2164/5048688 by DTU Library -Technical Information Center of Denmark user on 13 July 2018 Geocentre motion in terms of cartesian coordinates (in mm yr −1 ) estimated between the CE origin realized by each individual GIA forward model solution and the CM origin realized by the GPS data set, that is CE-CM.
Downloaded from https://academic.oup.com/gji/article-abstract/214/3/2164/5048688 by DTU Library -Technical Information Center of Denmark user on 13 July 2018 Martín-Español et al. (2016a)k to develop and extend the approach ofMartín-Español et al. (2016a).The framework will combine global geodetic data on GIA, ice mass and hydrological mass changes, changes in sea level (e.g. from GRACE and altimetry) and prior information from geophysical models to allow new insights about the different contributors to sea level rise on a regional and global scale.The GM GPS data set is available at .M. Schumacher, J. Rougier, Z. Sha, and J.L. Bamber are grateful for the financial support provided by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement no 694188 (GlobalMass project).J.L. Bamber was also supported by a Leverhulme Trust Fellowship (RF-2016-718) and a Royal Society Wolfson Research Merit Award.M.A. King was funded by an Australian Research Council Future Fellowship (project number FT110100207) and supported by the Australian Research Council Special Research Initiative for Antarctic Gateway Partnership (Project ID SR140300001).We would like to thank C. Zhang, J.Y. Guo, W. Huang, C.K. Shum and W. van der Wal for providing the ensemble of 13 GIA forward model solutions.
A C K N O W L E D G E M E N T S