Decoding the origins of vertical land motions observed today at coasts

Abstract

In recent decades, geodetic techniques have allowed detecting vertical land motions and sea-level changes of a few millimetres per year, based on measurements taken at the coast (tide gauges), on board of satellite platforms (satellite altimetry) or both (Global Navigation Satellite System). Here, contemporary vertical land motions are analysed from January 1993 to July 2013 at 849 globally distributed coastal sites. The vertical displacement of the coastal platform due to surface mass changes is modelled using elastic and viscoelastic Green’s functions. Special attention is paid to the effects of glacial isostatic adjustment induced by past and present-day ice melting. Various rheological and loading parameters are explored to provide a set of scenarios that could explain the coastal observations of vertical land motions globally. In well-instrumented regions, predicted vertical land motions explain more than 80 per cent of the variance observed at scales larger than a few hundred kilometres. Residual vertical land motions show a strong local variability, especially in the vicinity of plate boundaries due to the earthquake cycle. Significant residual signals are also observed at scales of a few hundred kilometres over nine well-instrumented regions forming observation windows on unmodelled geophysical processes. This study highlights the potential of our multitechnique database to detect geodynamical processes, driven by anthropogenic influence, surface mass changes (surface loading and glacial isostatic adjustment) and tectonic activity (including the earthquake cycle, sediment and volcanic loading, as well as regional tectonic constraints). Future improvements should be aimed at densifying the instrumental network and at investigating more thoroughly the uncertainties associated with glacial isostatic adjustment models.

Global change from geodesy, Satellite geodesy, Sea level change, Glaciology, Earthquake ground motions, Dynamics of lithosphere and mantle

1 INTRODUCTION

In response to changes from the solar system to its deep interior, the surface of the Earth is constantly moving up and down (Fig. 1). Under the influence of tides, the land surface oscillates with an amplitude of a few tens centimetres per day (e.g. Tomaschek 1957). Vertical land motions (VLMs) arise in response to changes of mass distribution within the atmosphere, oceans and continents forced by natural or anthropogenic mechanisms (e.g. Mangiarotti et al. 2001). Since the last glacial maximum (LGM ≈ 21 000 yr ago), the solid Earth responds to the melt of the continental ice sheets and is still in a state of isostatic disequilibrium (e.g. King et al. 2010). The ongoing glacial isostatic adjustment (GIA) causes changes in the Earth’s rotation rate, gravity field, geocentre and shape. The vertical displacement of the Earth’s crust can exceed 10 mm yr⁻¹ across previously glaciated areas (e.g. Sella et al. 2007; Milne et al. 2001). Driven by mantle convection, tectonic plates move at the Earth surface at the pace of a few centimetres per year. In the vicinity of plate boundaries, large VLMs are associated with seismic and volcanic events. During the earthquake cycle, tectonic stress accumulates along faults generating interseismic VLMs, which may range from ±1 to ±10 mm yr⁻¹ and last for decades (e.g. Ballu et al. 2011; Smith-Konter et al. 2014). The sudden stress release during the earthquake rupture may cause coseismic VLMs of tens of centimetres (e.g. Kanamori 1973). After major events, post-seismic VLMs may exceed ±10 mm yr⁻¹ to adjust stress changes in the crust and in the upper mantle (e.g. Freed et al. 2006; Satirapod et al. 2013). Vertically coherent deformation may be inferred down to depths of several hundred kilometres from geodetic, seismic and rock-mechanics observations (e.g. Bürgmann & Dresen 2008), suggesting that the mantle structure and flow play a dominant role in surface motions associated with earthquake cycles (e.g. Pollitz et al. 2001; Wang et al. 2012). Mantle convection also interacts with the Earth’s surface to generate a dynamic topography, responsible for metres per million years of elevation change over hundreds (Hoggard et al. 2016) to thousands of kilometres (e.g. Braun 2010; Flament et al. 2013).

Figure 1.

Causes of vertical land motions. The tides refer to solid Earth, ocean and pole tides mainly observed from hourly to annual time scales. The surface loading refers to atmospheric, oceanic and hydrological loading mainly observed at seasonal and yearly time scales. The anthropogenic influence refers to gas and groundwater pumping, dam building, mining and urbanization observed over the past 100 yr. The glacial isostatic adjustment (GIA) refers to the effects of present-day (PDIM) and past (PIM) ice mass changes, observed over the last 100 yr and 10 kyr, respectively. The tectonic activity refers to the earthquake cycle, sediment and volcanic loading and to the regional tectonic forces, observed over the past 10 yr, 10 kyr, 100 kyr and 10 myr, respectively. The mantle dynamics refers to the mantle control over the tectonic activity and to the dynamic topography, observed over the past 10 to 100 myr.

In coastal areas, sea-level variations result from a complex mix of climatic, oceanic and geophysical processes driven by natural and anthropogenic factors. Since tide gauges are bound to the coasts, their records are directly affected by VLM, which contributes to the considerable scatter of global mean sea-level rise so far assessed using the few long records available (e.g. Douglas 1991; Spada & Galassi 2012). In addition, VLMs contribute actively to the sea-level changes felt by coastal populations, as they can amplify, diminish or counteract the effects of climate-induced signals (e.g. Pfeffer & Allemand 2016). In many cases, for example, Torres Islands (Ballu et al. 2011), western Tropical Pacific (Becker et al. 2012), southern Europe (Wöppelmann & Marcos 2012) or Indian Ocean (Palanisamy et al. 2014), VLMs have been recognized as a dominant component of the total relative sea-level variations observed at coasts. An accurate determination of coastal VLMs is therefore necessary to understand the processes which contribute to sea-level changes from local to global scales, to appraise their impacts on coastal populations and make future predictions (e.g. Nicholls & Cazenave 2010).

In recent decades, several methods have been developed to evaluate coastal VLM based on direct GPS measurements (e.g. Wöppelmann et al. 2007), on combinations of geodetic observations (e.g. Cazenave et al. 1999), or on geodynamical models simulating tides (e.g. Scherneck 1991), surface loading (e.g. Mangiarotti et al. 2001) and GIA (e.g. Guo et al. 2012). While tidal and surface loading displacements due to atmospheric, oceanic and continental water mass circulation are predicted with an accuracy better than 1 mm yr⁻¹ (e.g. Penna et al. 2008; Santamaría-Gómez & Mémin 2015), much uncertainty remains on the response of the solid Earth to past ice mass (PIM) and present-day ice mass (PDIM) changes. Current models are still limited by an incomplete knowledge of the evolution of the cryosphere since the LGM and by inaccurate knowledge of the Earth’s rheological properties on time scales comparable with the Maxwell time of the mantle (e.g. King et al. 2010). Models imperfections result in a large range of predictions of surface displacements associated with PIM (e.g. Guo et al. 2012) and PDIM (e.g. Fleming et al. 2012) changes.

Geodetic measurements can bring new constraints on GIA models, both to explore the structural features of the Earth’s interior and to reconstruct past and present-day ice-melting history (e.g. Spada et al. 2012b; Mémin et al. 2014). This feedback is however limited by the sparsity of geodetic measurements and their limited accuracy (e.g. King et al. 2010). Besides, VLMs are influenced by various geodynamical processes (Fig. 1), making it pointless to search for a perfect agreement of observations with GIA models. Geodetic measurements contain valuable information about the vertical signature of plate tectonics associated with the earthquake cycle (e.g. Ballu et al. 2011; Satirapod et al. 2013; Smith-Konter et al. 2014) and regional geodynamical constraints (e.g. Kooi et al. 1998; Barlow et al. 2012; Serpelloni et al. 2013). The interactions between the lithosphere and the mantle explain much of the plate motions (e.g. Hager & O’Connell 1981), although the extent to which the mantle controls present-day surface deformation remains a major open question.

In this study, we take advantage of the global geodetic database ALTIGAPS, merging observations from multisatellite radar altimetry, tide gauges and GPS acquired over the past 20.6 yr (Pfeffer & Allemand 2016). VLMs are determined for 849 coastal sites with average error rates of about 0.6 mm yr⁻¹. Then, some of the causes of vertical coastal motions are investigated by modelling. The vertical displacements associated with the redistribution of sea, air, water and ice masses at the Earth’s surface are successively modelled using elastic and viscoelastic Green’s functions. To account for uncertainties on the GIA predictions, a wide range of rheological and ice loading (especially for PDIM changes) parameters are explored to predict VLMs. After that, models predictions are compared with observed VLMs, not to find the best model, but to investigate an ensemble of scenarios that could explain contemporary geodetic measurements. That way, observation windows are opened on active tectonic processes taking place within plates and at their margins. The results show that GIA corrections cannot alone account for vertical coastal motions affecting tide gauges. In turn, geodetic observations contain a wealth of information on global natural processes, involving the Earth’s external envelopes, the solid Earth and their interactions.

2 VLM OBSERVATIONS

VLM observations come from the ALTIGAPS database (Pfeffer & Allemand 2016), freely available online (http://julia-pfeffer.weebly.com). The ALTIGAPS database gives VLM trends and uncertainties from two different sources: (i) GPS measurements and (ii) a combination of satellite radar altimetry with tide gauges measurements. Different processing techniques are applied, which are in good agreement at global scale (fig. 6 from Pfeffer & Allemand 2016). Using several geodetic techniques allows covering a larger instrumental network at global scale, including 258 GPS stations and 591 tide gauges (Fig. 2). Some regions of the world remain poorly covered including South America, Africa, the Middle-East, South Asia and high-latitude regions.

Figure 2.

Locations of the measurements sites. The black crosses correspond to GPS station and the black dots to tide gauges. The coloured dots are the centres of the disks used to model PDIM changes. The colour scale corresponds to the ice-melting rates.

2.1 GPS data

At GPS stations, the VLM trends and uncertainties are estimated using the ULR5 solution provided by the University of La Rochelle Consortium (Santamaría-Gómez et al. 2012) available on the FTP server of SONEL (ftp.sonel.org). The data set covers the period from January 1995 to December 2010, and ensures a minimum length of 3 yr without discontinuities (considering both position and velocity offsets) and with data gaps not exceeding 30 per cent. The average length of the GPS time-series is about 10 yr. Vertical velocities are evaluated as the linear rate of change of height for 258 GPS stations. Formal uncertainties are assigned to take into account the noise content (parametrized by a power law plus a variable white noise model) in the residual position time-series with the maximum likelihood estimation technique. The uncertainties are given within one standard deviation (assuming a normal distribution, the 95 per cent confidence interval is given by 1.96 standard deviation). The same convention is kept for all VLM observations and predictions.

2.2 Combination of multisatellite radar altimetry with tide gauge data

Tide gauges are bound to the coast and measure the offset between the sea level and the land (i.e. relative sea level). Satellite radar altimetry observes sea level in a globally defined system referenced to the Earth’s centre of mass (i.e. absolute sea level). The combination of tide gauge with satellite altimetry measurements therefore enables (i) to isolate VLM, (ii) to place the results in the reference system used by satellite radar altimetry (ITRF 2008). VLM trends and uncertainties are extracted from the ALTIGAPS database using the same processing techniques as in Pfeffer & Allemand (2016).

2.2.1 Multisatellite radar altimetry data

Absolute sea-level changes are estimated using the weekly multimission grids produced by Ssalto/Duacs (DT-MSLA ‘All sat merged’) and distributed by AVISO¹ with support from the CNES (Centre National d’Études Spatiales). All classical corrections (which do not include GIA) are already applied to the grids. To ensure consistency with tide gauge data, the tidal and the atmospheric corrections applied by default to the DT-MSLA product are removed (Pfeffer & Allemand 2016). Absolute sea-level changes are then averaged within 1° of distance from tide gauges, and monthly averaged.

2.2.2 Tide gauge data

Relative sea-level changes are estimated using the Revised Local Reference (RLR) data from the Permanent Service for Mean Sea Level (Holgate et al. 2012; PSMSL, Permanent Service for Mean Sea Level 2014). Monthly data are selected at 591 tide gauges to fulfil a minimum completeness of 70 per cent. Data flagged for errors are not considered. To ensure consistency with satellite radar altimetry data, no correction has been applied for the tides, the atmosphere, nor for GIA (Pfeffer & Allemand 2016).

2.2.3 Linear trend estimation

Monthly VLMs are estimated at 591 tide gauges, based on the difference between absolute and relative sea levels over their common record period (January 1993–July 2013). Corrections for solid Earth tides are not applied since the monthly average eliminates the six largest components of the tidal signal (Chelton & Enfield 1986). Gaps in the time-series are filled using Nonequispaced Fast Fourier Transforms (NFTT) to remove high frequencies signals (periods shorter than 3 months) (Keiner et al. 2009). Seasonal signals are removed by least square adjustment of annual and semi-annual sinusoids. Long term trends are estimated simultaneously by least square adjustment of a linear trend. Individual error rates are estimated from the residuals of the regression and adjusted for serial correlation using an autoregressive model of order 1 (AR-1) (Pfeffer & Allemand 2016). The AR-1 noise model is shown to handle efficiently time correlation in annual sea-level data, but underestimated the errors in monthly data (Bos et al. 2014). Accurate error rates are retrieved with monthly sea-level data, when high frequencies signals are removed a priori (Pfeffer & Allemand 2016).

2.3 Combined VLM data set

VLMs are estimated at 258 GPS stations and 591 tide gauges, for a total of 849 coastal sites. Trend values range from −13 to +16 mm yr⁻¹ (Fig. 3a). A strong spatial variability in VLM is observed at all scales, including locally, which is evidenced for example by the dense instrumental network in Japan. High rates of crustal uplift are observed at high latitudes, while subsidence of coastal areas is more often observed at medium latitudes (e.g. the East American coast). Individual error rates range from 0.1 to 7.4 mm yr⁻¹, with a median value of 0.6 mm yr⁻¹ (Fig. 3b). The median error is lower at GPS stations (0.3 mm yr⁻¹) than tide gauges (0.8 mm yr⁻¹). Error rates only exceed 3 mm yr⁻¹ at 6 tide gauges.

Figure 3.

Vertical land motions observations: trend (a) and error rates (b).

Among the 591 tide gauges considered in this study, 113 are located at less than 10 km from a GPS station. The VLM trend values are obtained with independent techniques and are compared assuming that they should converge at neighbouring sites. The comparison indicate no bias (average difference = 0.05 mm yr⁻¹), high correlation coefficient (R = 0.77) and a root mean square error lower than 1.7 mm yr⁻¹ (Pfeffer & Allemand 2016). Part of the differences may be due to the use of different data periods for GPS and tide gauges, and to the distance (up to 10 km) separating ‘neighbouring’ stations. High spatial variability is indeed evidenced by the VLM observations. The root mean square error (1.5 mm yr⁻¹) slightly diminishes when considering shorter (1 km) distances between GPS stations and tide gauges (Pfeffer & Allemand 2016). An underestimation of individual error rates should also be envisioned (Pfeffer & Allemand 2016), though good agreement is found with recommended methods of error estimation (Bos et al. 2014).

3 VLM PREDICTIONS

Green’s functions are commonly used to calculate the vertical displacement of the Earth’s surface induced by loading (Farrell 1972). They characterize the Earth’s response to localized surface loads. The vertical displacement produced by a finite sized load is obtained by convolving the density of the load with the Green’s functions for a spherically symmetric, non-rotating and isotropic Earth:

\begin{equation} u_{\text{load}}(\theta , \lambda ,t)=\int _{-\infty }^t \text{d}t^{\prime } \int _{\Omega } G(\alpha ,t-t^{\prime }) {\cal L}(\theta ^{\prime }, \lambda ^{\prime },t^{\prime }) \text{d}A^{\prime } \end{equation}

(1)

with:

\begin{equation} \cos \alpha =\cos \theta \cos \theta ^{\prime }+ \sin \theta \sin \theta ^{\prime }\cos (\lambda ^{\prime }-\lambda ) \end{equation}

(2)

where u_load is the VLM predicted at colatitude θ, longitude λ and time t, G is the Green’s function, |${\cal L}$| is the load function, dA is the area element and Ω is the Earth’s surface. In the elastic case, the predicted VLM is instantaneous, so that G does not contain delayed viscoelastic effects (Spada & Stocchi 2007).

In order to model the vertical displacement of the Earth’s surface, we must have information on (1) the spatio-temporal variations of the load and (2) the rheological structure of the Earth, essential to compute Green’s functions. Model predictions generated with a range of rheological parameters and PDIM estimates result in a range of variations around the predicted average VLM (Table 1).

Table 1.

Parameters of the models used to predict VLM.

VLM source	Earth model	Load model
Surface loading	PREM	Atmosphere: ECMWF
	(SNREI type)	Ocean: ECCO
		Hydrology: ERA interim
PDIM	PREM	Uniform rates of ice-melt
	(SNREI type)	for eight glaciated regions
PIM	27 radial models	ICE-5G
	centred on VM2

3.1 Surface loading

The vertical displacements due to the redistribution of fluids at the Earth’s surface are estimated using atmospheric, oceanic and hydrological general circulation models. For the atmosphere, the operational European Center for Medium-range Weather Forecasts (ECMWF; Molteni et al. 1996) were used assuming an inverted barometer oceanic response. For the oceans, the Estimating the Circulation and Climate of the Ocean (ECCO) model was used (Wunsch et al. 2009), and ERA-interim was used for hydrology (Dee et al. 2011). The computations are based on the Preliminary Reference Earth Model (PREM) assuming a Spherically symmetric, Non Rotating, Elastic and Isotropic (SNREI) Earth (Dziewonski & Anderson 1981). Predicted 3-D surface displacements are provided in the centre of mass (CM) reference frame. Monthly averages are performed based on 6-hourly and 12-hourly data provided by the EOST loading service at the University of Strasbourg (http://loading.u-strasbg.fr/). Trends and associated error rates are then estimated using the same method described in Section 2.2.3 and the same period of time (January 1993–July 2013). The statistical distribution of the VLM predicted by surface loading is given in Fig. 4.

Figure 4.

Box plot of the predicted vertical land motions. Black circles are median values. Blue rectangles extend from the first to the third quartiles. Blue lines are whiskers, extending from the lowest value to the first quartile and from the third quartile to the largest value.

Predicted VLMs due to surface loading are lower than 0.3 mm yr⁻¹ in coastal areas (Fig. 4). The highest contribution comes from hydrological loading and reaches −0.21 mm yr⁻¹ in Antarctica. Ocean loading is slightly lower (−0.17 mm yr⁻¹ at Groote Eylandt, Australia), while atmospheric pressure loading only reaches 0.05 mm yr⁻¹. Averaged at global scale, the sum of these effects turns out to be close to zero (0.03 mm yr⁻¹). The effect of surface loading is small (a few 0.1 mm yr⁻¹), but not negligible compared to the error rates on VLM observations. Surface loading should therefore be accounted for to interpret correctly the VLM signal.

3.2 GIA associated with PDIM changes

The vertical displacement of the Earth’s surface induced by present-day ice melting is obtained by solving the sea-level equation (Farrell & Clark 1976) with an improved version (Spada et al. 2012a) of the open source program SELEN (Spada & Stocchi 2007). The Green’s functions are computed assuming a PREM structure (Table 1) and accounting for the effects of Earth’s rotation. The maximum harmonic degree is 128, corresponding to a horizontal resolution of ∼300 km. Predicted VLMs are provided in the CM reference frame.

VLMs are predicted for eight glaciated regions (Table 2) assuming a uniformly distributed ice loss. The melt rates are based on reconciled mass budgets from satellite gravimetry, satellite radar altimetry and from local glaciological records (Shepherd et al. 2012; Gardner et al. 2013) synthesized by Mémin et al. (2014). For Greenland and Antarctica, PDIM changes exhibit large error rates (Table 2), due to the nonlinearities of the melt over the past two decades and to the discrepancies between the various data sets available (table 1 and fig. 4 from Shepherd et al. 2012). Since the Earth model is elastic, large errors on PDIM changes are linearly mapped onto predicted VLMs. Final uncertainties on VLM predictions are therefore accounting for nonlinear melt rates over Greenland an Antarctica. For smaller glaciers (Alaska, Canada, Patagonia, Europe, Svalbard, Arctic), the PDIM changes were evaluated over a shorter time period (2003–2009) (Gardner et al. 2013), so that the uncertainty due to nonlinear melt rates may be underestimated (Table 2). PDIM changes were however shown to be relatively linear for the Alaskan, Patagonian, European and Arctic glaciers since the mid-1990s (fig. 4 from Dyurgerov & Meier 2005). The variation of the PDIM changes from 2003-2009 to 1993-2013 should therefore have little impact on our results. This assumption is not valid for the Andean and Canadian glaciers, where large variations in the ice melt rates were observed over the past two decades. However, few stations are located in the surroundings of these glaciers (Fig. 2), resulting in negligible effect on predicted VLMs (Supporting Information Appendix 1).

Table 2.

PDIM changes for the main glaciated regions based on Shepherd et al. (2012) and Gardner et al. (2013).

Region	Ice-melt	Error	Relative error
	(Gt yr−1)	(Gt yr−1)	(per cent)
Greenland	142	49	35
Antarctica	71	53	75
Alaska	64	20	31
Canada	60	8	13
Patagonia	29	10	34
Europe	15	2	13
Svalbard	12	2	10
Arctic	11	4	60

Predicted VLMs (Fig. 5) depend on the ice-melt rate (Table 2) and on the location of the stations (Fig. 2). The greater the melting is and the closer to glaciated regions the stations are, the higher predicted VLMs are. VLM predictions due to PDIM range from - 0.3 to 2.6 mm yr⁻¹, with error rates reaching 0.8 mm yr⁻¹ (Fig. 5). For most (90 per cent) of the stations (located relatively far from ice thinning glaciers), VLM predictions are ≤0.5 mm yr⁻¹ (Fig. 2). Predicted uplift rates exceed 1 mm yr⁻¹ for 34 stations located in Greenland, Alaska, West Canada, Svalbard and Norway (black stars in Fig. 5). For a detailed description of the surface displacement occurring in these regions, higher resolution models should be employed to accurately represent the PDIM changes in space and time (see Section 4.5). Elsewhere, VLMs due to PDIM changes are predicted with reasonable uncertainty compared to the observed signal (Fig. 3). For a detailed description of the contribution of each glaciated region, the reader is invited to refer to the Supporting Information (Appendix 1). The statistical distribution of the predicted VLM due to PDIM is given in Fig. 4.

Figure 5.

Predicted vertical land motions due to present-day ice mass (PDIM) changes: trend (a) and error rates (b). Black stars identify stations where predicted VLM exceed 1 mm yr⁻¹.

3.2.1 GIA associated with PIM changes

An ensemble of models is run using SELEN (Spada et al. 2012a) to predict at all sites the VLM due to the GIA following the melt of past ice sheets. In all runs, the ICE-5G (Peltier 2004) melting history has been employed. The earth model has an elastic lithosphere, a three layer mantle and a fluid homogeneous core. To account for the uncertainties in the Earth’s viscosity profile, each layer has a different viscosity value that represents the upper, average and lower bounds to the VM2 profile (Peltier 2004; Fig. 6). The lower mantle takes viscosity values of 1.6, 2.8 and 4.0 × 10²¹ Pa s. The upper mantle and transition zone takes viscosity values of 0.3, 0.4 and 0.5 × 10²¹ Pa s. The thickness of the lithosphere varies from 70, 90 to 110 km. The Earth model is incompressible and includes the rotational feedback on sea level (Spada 2017). Predicted VLMs are provided in the CM reference frame with a resolution corresponding to a maximum harmonic degree of 128.

Figure 6.

Viscosity profiles used in GIA models predicting the effects of PIM changes.

A total of 27 SELEN runs have been performed. At each site, the predicted VLM value corresponds to the average of the ensemble results. The standard deviation of the distribution is computed to estimate the error on predicted VLM values that reflects imperfect knowledge on the viscosity profile. The GIA associated with PIM changes account for the major part (coefficient of determination R² = 0.91) of the total VLM predictions (including surface loading, PDIM and PIM). According to our runs, the melt of former ice-caps generates VLMs ranging from −3.5 to 14.7 mm yr⁻¹ (Fig. 4), with error rates up to 4.8 mm yr⁻¹ (median value = 0.3 mm yr⁻¹). The statistical distribution of the predicted VLMs due to PIM is given in Fig. 4.

3.3 Total VLM predictions

Total VLM predictions (u^* in Fig. 7) are computed as the sum of VLM predictions associated with surface loading, PDIM and PIM changes:

\begin{equation} u^{*}=u_{\text{SL}}+u_{\text{PDIM}}+u_{\text{PIM}}, \end{equation}

(3)

where u^* ranges between −3.0 and 15.5 mm yr⁻¹ (median and average values are 0.0 and 0.3 mm yr⁻¹ respectively). The error rates are computed by adding the errors on surface loading, PDIM and PIM predictions in quadrature:

\begin{equation} \sigma _{u^{*}}=\sqrt{\sigma _{u_{\text{SL}}}^{2}+\sigma _{u_{\text{PDIM}}}^{2}+\sigma _{u_{\text{PIM}}}^{2}}. \end{equation}

(4)

Except for one station located along the Hudson Bay (⁠|$\sigma _{u^{*}}$| = 4.8 mm yr⁻¹), the error rates (⁠|$\sigma _{u^{*}}$|⁠) remain under 3 mm yr⁻¹ (median and average values are 0.3 and 0.4 mm yr⁻¹ respectively). The spatial distributions of |$\sigma _{u^{*}}$| are correlated with u^*; both display a long-wavelength spatial variability, with high values at high latitudes near past and present-day glaciated regions. The error rates consist in a range of VLM predictions around the average value, accounting for a range of viscosity and PDIM values.

Figure 7.

Total predicted vertical land motions: trend (a) and error rates (b).

4 COMPARISON OF VLM PREDICTIONS WITH OBSERVATIONS

4.1 Spatial filter

To be compared to the large scale features shown in VLM predictions (Fig. 7), observed VLMs are smoothed by averaging the values at the stations located at distances d_i ≤ L:

\begin{equation} u_{L} = \frac{\sum \limits _i(w_{i} u_{i})}{\sum \limits _i w_{i}}, \end{equation}

(5)

with weights

\begin{equation} w_{i} = \cos \left({\frac{\pi d_{i}}{2L}}\right), \end{equation}

(6)

where L is the averaging length. Assuming that the observed VLM trends u_i are independent, the error on the smoothed VLM trend is

\begin{equation} \sigma _{u_{L}} = \frac{\sqrt{\sum \limits _i(w_{i}^{2} \sigma _{u_{i}}^{2})}}{\sum \limits _i w_{i}} . \end{equation}

(7)

Different kind of spatial filters were designed to compare sparse geodetic observations to large-scale GIA predictions. For example, Calais et al. (2006) introduced the measurement accuracy in the weighting approach. Using eqs (5)–(7), Lebedev et al. (2003) have shown that the spatial filter is directly dependent on the network geometry. As L grows, VLM trend values are averaged over a larger number of stations to produce a smoothed value, which reduces the noise due to the local variability (Fig. 8). The efficiency of the weighted average filter depends however on the presence of nearby stations. Remote stations may be removed by imposing a minimum threshold w_min on the sum of weights, such as |$\sum \nolimits _i w_{i}$| > w_min. For w_min = 1, all the stations are included in the data set. As w_min grows, more stations are removed, leaving only the densely instrumented regions (US, Europe, Japan and Australia) in the data set (Fig. 8). If the averaging length is too short, high thresholds will purge the data set from a large part of its stations (Table 3). To limit the conflicts between two competing parameters (w_min and L), only the data sets containing at least 400 stations (Table 3) are considered in the following.

Figure 8.

VLM observations smoothed at 500 km: trend (a) and error rates (b). Imposing a threshold (w_min) on the spatial filter leads to exclude some stations from the data set. For w_min = 2, sole the stations circled in green, red and blue remain in the data set. For w_min = 5, sole the stations circled in blue and red remain. For w_min = 7, sole the stations circled in red remain.

Table 3.

Number of stations remaining in the smoothed VLM data set.

	wmin	wmin	wmin	wmin	wmin	wmin	wmin	wmin	wmin	wmin
L (km)	1	2	3	4	5	6	7	8	9	10
100	849	451	268	148	87	62	47	33	17	5
200	849	626	525	411	339	251	214	190	152	130
300	849	697	630	545	481	413	348	317	296	267
400	849	726	676	619	567	514	460	404	363	345
500	849	742	697	657	613	573	538	488	440	394
600	849	755	721	675	647	601	582	557	518	464
700	849	770	737	693	668	636	609	579	560	539
800	849	776	755	706	685	655	628	595	591	568
900	849	787	760	729	694	668	645	614	600	590
1000	849	799	761	739	704	687	659	626	616	600

4.2 Basic statistics

Taking into account the full data set (w_min = 1, Fig. 9), the maximum correlation coefficient (R = 0.7) is found for L = 500 km. At this averaging length, VLM predictions are able to explain half of the observed variance (R² = 0.49 in Fig. 9c) and root-mean-square error (RMSE) reaches its minimum value (1.6 mm yr⁻¹ in Fig. 9b). The average bias is lower than 0.3 mm yr⁻¹ regardless of the averaging length (Fig. 9d). The difference between predictions and observations reaches 10 mm yr⁻¹ at most at the Yakutat station (Fig. 10), due to the difficulty of reproducing the large crustal uplifts observed in Alaska (Fig. 3).

Figure 9.

Comparisons of vertical land motions predictions with observations: correlation coefficient (a), root mean square error (b), determination coefficient (c) and average bias (d).

Figure 10.

Difference between vertical land motions observations and predictions: residual trend (a) and cumulative error rates (b). Major plate boundaries from NUVEL-1 (DeMets et al. 1990) are shown with red lines.

Excluding the most isolated stations (w_min > 7), VLM predictions explain more than 80 per cent of the observed variance (Fig. 9c). The high crustal uplifts observed in Fennoscandia and Greenland are well reproduced by models. The moderate uplifts observed in Northern Europe and South America are also reasonably well predicted, as is the subsidence observed in North America below 50°N (Figs 3 and 7). It is expected for VLM predictions to be closer to observations in densely instrumented regions. Indeed, where more geodetic observations are available, better constraints are brought to the PIM, PDIM and rheological models, resulting in more accurate VLM predictions. Yet, even in well instrumented regions, some differences remain between VLM predictions and observations. Larger biases (up to 0.7 mm yr⁻¹) are observed with higher thresholds, which is due to (i) an over-smoothing of the VLM observations (high w_min are reached at large L values, see Table 3) and (ii) the reduction of the instrumental network to densely instrumented region (Fig. 8), where unmodelled geophysical processes (Section 4.4) may distort the comparison of VLM predictions with observations. In particular, greater importance is attributed to the VLM observations in Japan and West USA, which are known to be active tectonic regions.

4.3 Significance ratio

Residual VLM signal is observed worldwide (Fig. 10). VLM predictions do not explain the totality of the observed signal because of (i) imperfections in the models and in the instrumental network and (ii) the sources of displacement that have not been considered by modelling. In the first case, the consistency between VLM observations and predictions may be improved with the investigation of new rheological and ice-load models and with the densification of the instrumental network. In the second case, looking for a perfect agreement between predictions and observations would be pointless, as various geodynamical processes generate significant vertical displacement over a broad range of spatial and temporal scales (Fig. 1, Table 4). To identify the significant residual VLM, we compute an adimensional significance ratio (SR):

\begin{equation} \text{SR} = \sqrt{\frac{(u_{500}-u^{*})^{2}}{\sigma _{u_{500}}^{2}+\sigma _{u^{*}}^{2}}}, \end{equation}

(8)

where u₅₀₀, |$\sigma _{u_{500}}$|⁠, u^* and |$\sigma _u^*$| are given by eqs (5), (7), (3) and (4), respectively. Since the results are presented in a map (Fig. 11), all VLM data are kept (w_min = 1) as isolated stations can be easily identified.

Figure 11.

Non-dimensional significance ratio (SR). Black rectangles identify well-instrumented regions with significant residual VLM signal.

Table 4.

Order of magnitude of the VLM signal expected at decadal time scales. The references give some examples of the observed and predicted VLM and do not aim to provide a complete review of the associated geophysical processes.

Geophysical process	VLM signal	Scale of available models	Complementary evidence in the literature
Surface loading	0.1 mm yr−1	Global	e.g. Blewitt & Lavallée (2002)
			Santamaría-Gómez & Mémin (2015)
GIA (PDIM)	1 mm yr−1	Global	e.g. Fleming et al. (2012); Spada et al. (2012b)
			Mémin et al. (2014); Melini et al. (2015)
GIA (PIM)	1–10 mm yr−1	Global	e.g. Milne et al. (2001); Sella et al. (2007)
			King et al. (2010); Guo et al. (2012)
Eathquake cycle	1–10 mm yr−1	Local to regional	e.g. Fialko (2004); Smith-Konter et al. (2014)
			Ozawa et al. (2012); Satirapod et al. (2013)
Volcanic loading	1–10 mm yr−1	Local to regional	e.g. Moore (1970); Caccamise et al. (2005)
			Zhong & Watts (2013); Huppert et al. (2015)
Sediment loading	1–10 mm yr−1	Local to regional	e.g. Ivins et al. (2007); Törnqvist et al. (2008)
			Shennan et al. (2012); Kooi et al. (1998)
Tectonic forces	1–10 mm yr−1	Local to regional	e.g. Dokka (2006); Freymueller et al. (2008)
			Serpelloni et al. (2013); Anzidei et al. (2014)
Man influence	1–100 mm yr−1	Local to regional	e.g. Galloway & Burbey (2011)
			Syvitski et al. (2009); Cooper (2002)

High SR values are observed over nine well-instrumented regions, namely: Japan, Australia, West USA, Southeast USA, Great Britain, Hawaiian islands, Alaska, Southeast Asia and Northeast Mediterranean basin. In this case, the residual VLM signal can be explained neither by the uncertainty on observations, nor by the range of variability considered in VLM predictions. On the opposite, the residual VLM signal observed in the largest part of Europe and in Northeast America is comprised within the cumulative uncertainty on models and observations.

4.4 Observation windows on unmodelled geophysical processes

Significant differences between VLM observations and predictions are observed over nine well instrumented regions of the world (Fig. 11), defining observations windows on the geophysical processes unaccounted by global geodynamical models (Table 4). Zooms in these regions are provided in Supporting Information (Appendix 3). Such processes are not necessarily limited to these areas, but will elsewhere be masked by inaccurate VLM predictions and observations.

4.4.1 Anthropogenic influence

The exploitation of groundwater resources and the land modifications induced by human activities may generate significant VLM signal in coastal areas. In some places, fluid withdrawal such as groundwater, oil, or gas extraction can cause subsidence rates up to tens of cm yr⁻¹ (e.g. Galloway & Burbey 2011). Coastal areas, in particular deltas, are also sinking due to sediment shortage caused by dam building or sediment compaction caused by urbanization (e.g. Syvitski et al. 2009). The development of InSAR (Interferometric Synthetic Aperture Radar) techniques over the past two decades showed that man induced surface displacements are strongly nonlinear both in space and time (e.g. Bock et al. 2012; Chaussard et al. 2013; Tosi et al. 2016). Such processes are unlikely to explain the large scale patterns observed in residual VLM, but may contribute to the local variability observed in Japan, California and Great Britain (Fig. 10). Large subsidence rates are still observed in Wales due to the collapse of ancient chalk and coal mines (e.g. Cooper 2002; Bell et al. 2005). Besides, complex surface displacements associated with fault reactivation has been reported in British coal mines (Donnelly 2006). These phenomena, allied with gas, oil and groundwater extraction/storage (Bateson et al. 2015), may explain part of the residual VLM signal observed in Great Britain (Fig. 10). Various instrumental issues have also been faced at tide gauges in Great Britain, which may impact the observed VLM trends. Quality control flags and datum issues have been reported for the Cromer, Holyhead and Liverpool gauges. Besides, monthly values have been lost at the Barmouth, Llandudno and Newport gauges (total signal lost <30 per cent of the complete study period). Residual VLM signal should therefore be interpreted with caution for these stations and in short (<20 yr) tide gauge records in general (Wöppelmann & Marcos 2016).

4.4.2 Earthquake cycle

In the vicinity of plate boundaries, residual VLMs display large local variability, well evidenced by the dense instrumental networks available in Japan and the west coast of USA. Indeed, the release of strain by large earthquakes takes place in a variety of ways (Kanamori 1973), inducing VLMs before, during and after the seismic events (e.g. Ozawa et al. 2012). The surface displacements associated with the earthquake cycle have been discussed in detailed across Japan, since abundant seismic and geodetic data are available (e.g. Aubrey & Emery 1986; Sagiya et al. 2000; Suwa et al. 2006; Hashimoto et al. 2009; Ozawa et al. 2011). Along the San Andreas Fault System, interseismic strain accumulation was shown to cause VLMs in excess of 3 mm yr⁻¹ (e.g. Smith & Sandwell 2003, 2006), while vertical displacements of several tens of cm can be associated with the earthquakes events (e.g. Fialko 2004). Recent 3-D earthquake cycle deformation models were used to predict VLM along the San Andreas Fault System in good agreement with tide gauge observations (Smith-Konter et al. 2014). Significant residual VLMs are also observed in Southeast Asia (Figs 10 and 11), which are likely to be associated with the Sumatra–Andaman 2004 earthquake (Lay et al. 2005). Postseismic surface displacements are observed at several hundred kilometres from the trench (Trubienko et al. 2013), which is consistent with our results shown in Figs 10 and 11 and in theSupporting Information (Appendix 2).

4.4.3 Volcanic loads

Significant residual VLMs are observed over four islands (Hawaii, Maui, Oahu and Kaui) of the Hawaiian archipelago (Figs 10 and 11). Large residual subsidence rates are observed in Hawaii (from −1.3 ± 0.4 to −7.6 ± 0.9 mm yr⁻¹), which progressively decrease for older islands (from −0.1 ± 0.3 to −0.7 ± 0.3 mm yr⁻¹ in Maui). Slight residual uplift rates (from +0.2 ± 0.3 to +0.9 ± 0.2 mm yr⁻¹ in Oahu) are observed at large distances from the hotspot (from +0.2 ± 0.3 to +1.2 ± 0.4 mm yr⁻¹ in Kaui). This regional VLM pattern is related to the flexure of the lithosphere expected from volcanic loading at the Hawaiian hotspot (e.g. Moore 1970; Caccamise et al. 2005; Zhong & Watts 2013). Interactions with the mantle plume may also induce surface displacements (e.g. Zhong & Watts 2002), though the migration of the Hawaiian islands over the swell was shown to have minor effects on VLMs for Hawaii, Maui, Oahu and Kaui (e.g. Huppert et al. 2015).

4.4.4 Sediment loads

Sediment loading and unloading may cause significant residual VLM in coastal areas. Holocene sediment compaction have been shown to contribute for several mm yr⁻¹ of the subsidence observed in the Gulf of Mexico (e.g. Ivins et al. 2007; Törnqvist et al. 2008). Local tectonic processes (e.g. Dokka 2006) and fluid withdrawal (e.g. Kolker et al. 2011) were also invoked to explain the large subsidence rates observed in Texas and Louisiana. The moderate uplift rates observed in the North East of the Gulf were also confirmed by independent studies (Letetrel et al. 2015).

4.4.5 Tectonic forces

Significant residual VLMs are also observed at scales of a few hundred kilometres in active tectonic regions such as eastern Mediterranean or Alaska (Fig. 11).

The Mediterranean is undergoing complex surface displacements related to the convergence of Africa towards Eurasia, with possibly smaller intervening microplates (e.g. Carminati & Doglioni 2005; Anzidei et al. 2014). In the eastern Mediterranean, horizontal surface motions are dominated by a counterclockwise rotation of the Agean-Anatolian block with respect to Europe (e.g. Nocquet 2012). The residual subsidence observed in Greece may be related to diffuse extension from western Bulgaria up to the eastern Gulf of Corinthe (e.g. Pérouse et al. 2012). The retreat of the Hellenic slab under the Agean have also been suggested to explain the flow-like kinematics of the entire region from the South Balkan to continental Greece (e.g. Jolivet et al. 2015). At the North, the Adriatic region is often considered as an independent microplate bordered with active orogenic systems (e.g. Battaglia et al. 2004). At the northeast of the Adriatic sea, coastal areas are relatively stable or in moderate residual uplift, which may be related to the late orogenic strike-slip faulting in the South of the Alps (e.g. D’Agostino et al. 2008) and in the frontal zone of the Dinarides (e.g. Grenerczy et al. 2005). One should keep in mind that these zones are still moderately seismically active and that surface displacements are also expected from the earthquake cycle (e.g. Burton et al. 2004; Di Bucci & Angeloni 2013).

The situation is just as complex in Alaska, where abundant seismicity is produced by the subduction of the Pacific under the North American plates (e.g. Freymueller et al. 2008). Most of the surface displacements in Alaska can be described in terms of post-seismic deformation (e.g. Freed et al. 2006; Suito & Freymueller 2009), plate coupling and slip deficit (e.g. Fournier & Freymueller 2007), motions of large crustal blocks (e.g. Elliott et al. 2010), and slow-slip events (e.g. Fu et al. 2015). Evidence for deeper mantle influence on surface displacement has also been reported by Finzel et al. (2015). VLMs are also highly impacted by present-day ice melting and may be better explained with the use of high resolution models (e.g. Melini et al. 2015). Part of the uplift observed in southeast Alaska is also caused by the rapid viscoelastic relaxation due to ice unloading after the Little Ice Age (Larsen et al. 2005). This effect is not taken into account in our GIA predictions and may explain part of the residual signal observed. In spite of these uncertainties, significant residual uplift (SR > 10) is observed in Alaska (Fig. 11), which may be explained by one or a combination of the above mentioned processes.

4.4.6 Influence of the mantle

More surprisingly, significant residual VLMs are also observed at regional and continental scales along passive margins.

For example, significant residual uplift are observed along the Australian coasts. Residual uplift rates are greater in the South-East, but otherwise display remarkable spatial consistency at continental scale (Figs 10 and 11). Such observations are correlated with the surface topography showing high reliefs at the East, which supports the assumption of the tilting of Australia due to interactions with mantle convection processes (e.g. Heine et al. 2010). If the various geodynamical models reproducing inherited dynamic topography often diverge (e.g. Flament et al. 2013), they recognize the importance of mantle dynamics on the surface displacements observed along rivers profiles (e.g. Czarnota et al. 2014) or marine terraces (e.g. Sandiford 2007) in Australia.

The east coast of the United States is a typical passive continental margin, which is still slowly sinking from the forebulge collapse (e.g. Engelhart et al. 2009). Once corrected for GIA, significant residual VLMs remain in the South of the country (Figs 10 and 11). Thermally driven subsidence, sediment loading and compaction were already invoked to explain relative sea-level changes on secular (e.g. Simms et al. 2013) to 100 Myr time scales (e.g. Browning et al. 2008; Kominz et al. 2008). Recent mantle convection simulations also imply that dynamic topography may account for a large part of the geological structure of the southeastern US Atlantic coast (e.g. Spasojević et al. 2008; Rowley et al. 2013).

Such processes occur however on much longer time-scales (from a few million to a few hundred of million years), and are unlikely to explain the present-day VLM observed at coasts. The vertical displacement expected from dynamic topography is usually only a few metres per million years (e.g. Braun 2010; Flament et al. 2013), though faster rates (up to a few 0.1 mm yr⁻¹) have recently been inferred from stratigraphic observations (Hoggard et al. 2016).

4.5 Additional sources of uncertainties

Several sources of errors have to be considered in the VLM predictions induced by PIM changes, including the uncertainties related to (i) the Earth’s viscosity model, (ii) the spatial and temporal distribution of ice melting and (iii) the numerical implementation of the model features. In our analysis, we only explored the uncertainties on the Earth’s viscosity model. Here, we compared two set of VLM predictions taking into account two different ice models (ICE-5G and ICE-6G), two different viscosity profiles (VM2 and VM5a) and two different computing methods (Peltier 2004; Purcell et al. 2016). The resulting VLM predictions are similar to our analysis (Fig. 12). The largest differences in VLM predictions are observed in North America, especially in Greenland and northern Canada. Elsewhere, VLM predictions differ by less than 1 mm yr⁻¹. Regardless of the models considered, the same regions (Japan, Australia, West USA, Southeast USA, Great Britain, Hawaiian islands, Alaska, Southeast Asia and the Northeast of the Mediterranean basin) show significant residual VLM. In fact, when considering the models from Peltier (2004) and Purcell et al. (2016), the significance ratio is higher (see Supporting Information Appendix 2), which suggests that the residual VLM may be significant in other regions of the world.

Figure 12.

Comparison of VLM predicted with various GIA models.

Large errors can occur in the predictions of VLM associated with PDIM changes near thinning glaciers (Figs 2 and 5). The large residual VLM observed in Alaska are likely due for some part to inaccurate PDIM predictions (Figs 10 and 11), which should be performed at higher resolution (Bevis et al. 2016) and should also consider nonlinear PDIM variations. Large uplift rates (up to 7 mm yr⁻¹) were indeed predicted in Alaska using the REAR (Regional Elastic Rebound Calculator) software (Melini et al. 2015) and reconstructed PDIM changes taking into account the glaciers of the Randolph Glacier Inventory (Pfeffer et al. 2014) and the ice mass balance from Marzeion et al. (2012). Besides rapid (from secular to millennium time scales) viscous responses to ice mass changes should also be considered in future models (Larsen et al. 2005). Unaccounted PDIM changes cannot explain the residual VLM observed over Japan, Australia, West USA, Southeast USA, Great Britain, Hawaiian islands, Southeast Asia and Northeast Mediterranean basin (Fig. 11) since these regions are located far from ice thinning glaciers (Fig. 2).

Besides, numerous instrumental issues can affect the tide gauge and GPS records used to estimate VLM trends. While high accuracies can be reached for GPS records lasting at least 3 yr (Santamaría-Gómez et al. 2012), much longer records are needed for the combination of tide gauge and altimetry measurements. Wöppelmann & Marcos (2016) found that a common signal of 20 yr allowed reaching high accuracies (<1 mm yr⁻¹) on VLM trends. Our study period (20.6 year) is just above this limit, meaning that VLM trends of a few mm yr⁻¹ can be detected with our approach, but should be interpreted with care.

The uncertainties in the realization of the International Terrestrial Reference Frame (ITRF) may also have an impact on observed VLM rates (Wöppelmann & Marcos 2016). The accuracy of the ITRF2008, in which the VLM observations are referenced, was estimated at 0.5 mm yr⁻¹ (Collilieux et al. 2014). The stability and accuracy of the ITRF should therefore be considered for the interpretation of small (a few 0.1 mm yr⁻¹) geodetic signals (Blewitt 2015).

5 CONCLUSIONS AND PERSPECTIVES

Associated with a modelling approach, multitechnique geodetic database constitute an interesting tool to investigate the geophysical processes affecting the solid Earth, its external envelopes and their interactions. All VLM observations and predictions discussed here are provided in Supporting Information (Appendix 4).

The ALTIGAPS database evaluates VLM at coasts over the past two decades based on a combination of satellite radar altimetry, tide gauges and GPS measurements. The instrumental networks counts 849 stations unevenly distributed along the coastlines. The Southern hemisphere is sparsely covered, as well as high latitudinal regions. An ensemble of surface loading and GIA models was tested at global scale taking into account both PDIM and PIM changes. The VLM predictions explain about half of the total observed variance and 80 per cent of the large-scale variability (>300 km) observed in well instrumented regions of the world. VLM predictions follow closely the observations in Europe (including Fennoscandia) and North America. The history of past ice sheets is indeed reasonably well constrained in these regions, which leads to a better agreement of the GIA predictions with geodetic observations.

Significant residual VLMs are evidenced over nine densely instrumented regions. In good agreement with the literature, residual VLMs may be associated with the earthquake cycle (Japan, South East Asia, region of the San Andreas Fault), fluid withdrawal (Japan, California, Southeast Asia), coal mine exploitation (Great Britain), volcanic charges (Hawaii), sediment loading (Gulf of Mexico) and tectonic forces (Northeast Mediterranean basin, Alaska). The local to regional VLMs expected from these processes range from 1 to 10 mm yr⁻¹ (Table 4). Coastal subsidence may even reach 100 mm yr⁻¹ in densely populated areas, due to sediment shortage and compaction caused by fluid withdrawal, dam building and urbanization. In the absence of global models, the understanding of VLMs generated by man and by the tectonic activity (including the earthquake cycle, the sediment loading, the volcanic loading and regional tectonic constraints) relies on geodetic observations and on local to regional geophysical models (Fig. 1 and Table 4).

The understanding of the causes of VLMs at coasts is limited both by the availability of accurate geodetic data and geophysical models. A better spatial coverage would be obtained with the inclusion of additional geodetic measurements from satellite laser ranging, very-long-baseline interferometry and Doppler Orbitography and Radiopositioning Integrated by Satellite. Besides, various instrumental issues may corrupt the VLM signal extracted from the combination of tide gauge and altimetry measurements, implying that relatively short (<20 yr) records should be interpreted with caution. Longer period of observations could eventually be accessed with the use of differential tide gauges measurements available since the beginning of the 19th century (e.g. Kuo et al. 2008; Santamaría-Gómez et al. 2014; Letetrel et al. 2015). The interpretation of small (a few 0.1 mm yr⁻¹) geophysical signals is also limited by the current accuracy and stability of the ITRF (Collilieux et al. 2014; Blewitt 2015; Wöppelmann & Marcos 2016).

Concerning VLM predictions, higher resolutions models should be used near ice-thinning glaciers to better evaluate the VLM induced by PDIM changes in space and time. Rapid viscous responses to ice mass changes should also be taken into account to assess the impacts of the glaciers retreat since the Little Ice Age. Finally, a rigorous assessment of systematic errors in GIA modelling remains to be done, which would explore simultaneously the errors involved in the viscosity models, in the ice history models and in their numerical implementation. The deglaciation of Antarctica since the LGM is particularly poorly constrained (Whitehouse et al. 2012), which limits severely the relevance of GIA predictions at high latitudes in the Southern Hemisphere.

Acknowledgments

The authors thank two anonymous reviewers for their comments. This research benefited from financial support from the CNES (Centre National d’Etudes Spatiales, France) through the TOSCA committee fellowship and from the European Research Council within the framework of the SP2-Ideas Program ERC-2013-CoG, under ERC Grant agreement number 617588. GS is supported by a DiSPeA research grant (CUP H32I160000000005) and by Programma Nazionale di Ricerche in Antartide (PNRA 2013/B2.06, CUP D32I14000230005). AM was supported by an Australian Research Council Super Science Fellowship (FS110200045).

1

http://www.aviso.altimetry.fr/duacs/

SUPPORTING INFORMATION

Supplementary data are available at GJI online.

Appendix4.zip

gji_supplementary_1.pdf

gji_supplementary_2.pdf

gji_supplementary_3.pdf

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.