Present-day VLM and projected RSL change driven by deglacial mass change
Ice changes over the deglaciation (30ka-1ka) affect both present-day VLM and RSL projections. For predictions of both, we use output from a suite of 65 3D GIA model runs in which the Earth structure varies laterally and with depth (Supplementary Table 1). The majority of runs (60) are adopted from Pan et al.72. They vary in choice of two global ice histories, ICE-6G53 and ANU73, four seismic tomography models, which inform lateral variability in viscosity50,51,74,75, two models for variations in lithospheric thickness76,77, and two values for the scaling factor used to map between temperature perturbations in the mantle and viscosity (0.02/∘C and 0.04/∘C). The spherical average of each 3D viscosity profile is equivalent to the 1D Earth model generally paired with each ice history72. We note that for the ANU simulations, two different spherical average viscosity profiles were used to explore the non-uniqueness reported in Lambeck et al.73 (Supplementary Table 1).
We further adopt output from four 3D GIA model runs presented in Austermann et al.52 and one from Antwerpen et al.49. Austermann et al.52 inferred a model of Earth’s 3D viscosity structure from the Schaeffer & Lebedev50 seismic tomography model for the mantle above the transition zone (<400 km depth) and the French & Romanowicz51 model for the mantle below that. The Earth models explore differences in the average lithospheric thickness, average mantle viscosity, and scaling from temperature perturbation to viscosity perturbation to mimic dislocation creep. Earth structures are paired with the ice reconstruction from ICE-6G for 150 ka to present. Antwerpen et al.49 used the 3D viscosity perturbations from Austermann et al.52 paired with the Greenland-specific ice reconstruction Huy36. While Huy3 and the GrIS portion of ICE-6G were created using glaciological models that account for ice dynamics, ANU and the other global ice masses in ICE-6G and Huy3 (excluding the North American portion) were not. All 3D GIA calculations were performed with the GIA code SEAKON78, which is gravitationally self-consistent and solves for solid Earth deformation, geoid changes, and impacts to Earth’s rotational axis while accounting for the migration of shorelines3.
We use output from the suite of 65 GIA model runs described above to calculate the contribution of the deglaciation to present-day VLM and future RSL change. The rate of present-day VLM is taken as the rate of solid Earth deformation between 250 years before present (BP) and today. For future RSL change, we calculate the rate of RSL change between 250 years BP and today and assume that this is an appropriate estimate for the rate of change over this century. We extrapolate this rate to 2100 CE to obtain the deglacial contribution to future RSL. To account for the possibility that some model runs are more accurate representations of Earth’s internal structure and deglacial ice history than others, we weight the predictions from each model based on their fit to late-Holocene (7 ka to present) RSL observations around Greenland from the GAPSLIP database23. The dataset we use consists of both sea-level index points (n = 125), which give an estimate for the elevation of sea level at a point in the past, and limiting data points (n = 268), which can only tell us that sea level was above (marine-limiting) or below (terrestrial-limiting) a certain elevation. We use a weighted residual sum of squares (WRSS) calculation79,80 for each model run’s prediction to combine the uncertainties associated with the index and limiting data, according to:
$${{{{\rm{WRSS}}}}}_{{{{\rm{nm}}}}}=\left\{\begin{array}{ll}{(\frac{2{r}_{nm}}{{\varepsilon }_{n}})}^{2},\hfill & {q}_{n}=0 \hfill \\ – 2ln(0.5+0.5{{{\rm{erf}}}}({q}_{n}\frac{{r}_{nm}}{{\varepsilon }_{n}})),\hfill & {q}_{n} \ne 0\hfill\end{array}\right.$$
(1)
where erf represents the error function, rnm is the residual between the modeled and observed RSL at the age and location of data point n and model m, εn is the uncertainty in the observed RSL, and qn refers to the type of datum (1 for terrestrial limiting, −1 for marine limiting, and 0 for index). For a more detailed explanation of the equation for WRSSnm, please refer to Creel et al.79 and references within. The weight for each model is then calculated as the normalized inverse of WRSSnm summed over all data points. We use these weights to calculate the weighted average and weighted standard deviation of present-day VLM rates (Fig. 1). We also sample the deglacial contribution to future RSL change according to the model weights in order to combine them with the other contributions to projected RSL change.
Present-day VLM driven by ice changes since the beginning of the LIA
Ice changes since the beginning of the LIA (1ka–0ka) affect both present-day VLM and RSL projections. In this section, we describe the GIA model setup that captures the time period from the LIA to present-day, which we use to calculate present-day VLM. We calculate this contribution first since we can compare it to observations to better constrain uncertainties in the ice and Earth parameters of the GIA model (Table 1). These constraints are then propagated into our projections of RSL driven by ice-sheet changes from the LIA to 2100 CE.
Present-day VLM caused by the Earth’s response to LIA ice and ocean loading is predicted using a gravitationally self-consistent formalism that solves for solid Earth deformation, geoid changes, and impacts to Earth’s rotational axis while accounting for the migration of shorelines3. This 1D GIA model requires an input (1D) Earth structure and ice thickness history from the LIA to the beginning of the GPS record. Note that deformation associated with ice change during the GPS record is calculated separately and described in the next section. Calculations are run at spherical harmonic degree 512 for 12 timesteps from 1000 to 2022 CE. We vary lithospheric thickness, mantle viscosity, LIA mass anomaly, and LIA termination time following the posterior probability solutions for these parameters from Adhikari et al.16.
Adhikari et al.16 combined estimates of mass balance for the GrIS17 and its peripheral glaciers18 to create a history of post-LIA ice thickness anomaly relative to the year 2017 CE. We sample the mass anomaly of the ice sheet during the LIA from a range centered around the posterior mean (14,860 ± 2670 Gt) calculated by Adhikari et al.16 (Table 1, n = 13). Its magnitude is assumed to be constant throughout the duration of the LIA. The Greenland-wide termination time of the LIA is varied within the range from 1800–1910 CE (Table 1, n = 9). Adhikari et al.16 find that changing the preceding Medieval Warm Period (~ 1000 CE) mass anomaly and the LIA inception time do not affect their results significantly, so these parameters are kept constant in our models at 3673 Gt and 1450 CE, respectively. Ice change that approximately spans the time frame of GNET observations is set to zero to not double-count for coeval elastic deformation described in the next section.
We use a 1D homogeneous Maxwell model overlain by an elastic lithosphere, for which we vary the thickness. Our mantle viscosity ranges from 4 × 1019 Pa s to 25 × 1019 Pa s (Table 1, n = 12) and our lithospheric thickness varies from 60 to 240 km (Table 1, n = 10), with a higher sampling density around the highest posterior likelihood of each parameter16. The elastic and density structure vary with depth according to the Preliminary Reference Earth Model[PREM; 81. We assume a laterally and radially homogeneous viscosity structure because the inferred mantle viscosity over centennial timescales varies minimally for the upper mantle below Greenland41, and it is unlikely that small-scale load changes, like those experienced in Greenland since the beginning of the LIA, invoke deformation within the lower mantle42. However, we also acknowledge that omitting 3D variations in Earth structure, especially those driven by lateral variations in lithospheric thickness, may have some impact on our results. Combining all parameters leads to 14,040 GIA simulations.
Observations of present-day VLM and their elastic corrections
To better constrain the parameters described in Table 1, we compare predictions of VLM to observations. Here we describe the observations that are being used and additional corrections that need to be applied.
Observations of VLM were obtained from a network of 57 permanent GNSS stations located around Greenland’s coast, termed GNET. We use rates of VLM determined by Berg et al.20 that were calculated over the period of available data for each station. The majority of GNET stations were established between 2007 and 2009 CE, though there are some that were installed in the mid-to-late 1990s, and their observational records run to 2022 CE. Positive uplift rates are observed at all 57 stations, with rates that range in magnitude from 1.3 to 17.5 mm yr−1. Following Adhikari et al.16, we remove two stations located in central East Greenland near the Kangerlussuaq glacier from our analysis. The mantle viscosity beneath these locations is likely lowered locally due to mantle flow associated with the Icelandic plume16,26,82,83, which we don’t allow for in our calculations here. Observed uplift rates at each station (VOBS) are assumed to be the sum of VLM contributions due to ice-mass changes during the deglaciation (VDEG), ice-mass changes since the beginning of the LIA (VLIA), ice-mass changes during the GPS record (VELA) and a reference frame correction (VRF; Supplementary Section 1):
$${{{{\bf{V}}}}}_{{{{\bf{OBS}}}}}\approx {{{{\bf{V}}}}}_{{{{\bf{DEG}}}}}+{{{{\bf{V}}}}}_{{{{\bf{LIA}}}}}+{{{{\bf{V}}}}}_{{{{\bf{ELA}}}}}+{{{{\bf{V}}}}}_{{{{\bf{RF}}}}}$$
(2)
For the contribution of elastic deformation to present-day uplift rates, VELA, we use elastic VLM rates modeled by Berg et al.20. They obtain rates of mass change that span the last two decades for the GrIS84,85, Greenland’s peripheral glaciers21, and Canada’s peripheral glaciers20. In response to these ice changes, the solid Earth beneath Greenland has been uplifting at rates that range from 1.5 mm yr−1 to 13.8 mm yr−1. The GrIS is generally the largest contributor to these signals, though the contribution from peripheral glaciers can be significant, especially in North and East Greenland, where their presence is greatest20. The elastic, deglacial, and reference-frame contributions to present-day VLM rates are incorporated into the framework of our Bayesian analysis (next section), so that we can directly compare predictions of LIA-driven VLM to observations.
Bayesian framework
We find posterior distributions for our input parameters of the LIA-GIA model, summarized in Table 1, by combining model predictions of VLM with observations using a Bayesian inversion. Our approach is similar to that used by Piecuch et al.86 and uses the paradigm from Berliner87, which constructs a model with three levels of equations. First, process-level equations represent how the quantity of interest—here, the present-day VLM rate arising from the LIA—evolves in space and time. Second, data-level equations describe how the available observations relate to the underlying process. Finally, parameter-level equations place prior constraints on the model’s parameters. For more details on hierarchical modeling, see Cressie & Wikle88. Our posterior distributions are updates to those found in Adhikari et al.16 since we used improved GNET data and estimated elastic contributions20, which cover the longest time frames available for each GNET station and account for the elastic rebound caused by peripheral glacier melt from Greenland and Canada, as well as through our use of a more comprehensive Bayesian analysis to compare this data to our modeled predictions.
Data level
We have data \({{{\bf{x}}}}={[{x}_{1},…,{x}_{N}]}^{T}\) at N = 55 locations, which are GNSS rates corrected for elastic and deglacial contributions. Thus, these data reflect present-day VLM due to the LIA plus a reference frame correction (hereinafter LIA’ VLM). We assume that the data x are noisy versions of the true LIA’ VLM rate \({{{\bf{v}}}}={[{v}_{1},…,{v}_{N}]}^{T}\) according to:
$${{{\bf{x}}}} \sim {{{\mathcal{N}}}}({{{\bf{v}}}},\Delta )$$
(3)
where ~ means “is distributed as”, \({{{\mathcal{N}}}}(a,b)\) is the multivariate normal distribution with mean vector a and covariance matrix b, and Δ = diag(\({\delta }_{1}^{2},…,{\delta }_{N}^{2}\)) is the diagonal matrix of the data error variances where we define \({\delta }_{i}^{2}\) as the sum of error variance of the ith GNSS data and the models of the elastic and deglacial contributions.
Process Level
We assume that the true LIA’ VLM rates v arise from GIA, as well as a residual large-scale (regional) process and a residual small-scale (local) process, both unrelated to GIA. Formally, we write:
$${{{\bf{v}}}} \sim {{{\mathcal{N}}}}({{{\bf{u}}}},{\epsilon }^{2}{\mathsf{I}})$$
(4)
where ϵ2 is the unknown spatial variance of the local VLM process unrelated to GIA, \({\mathsf{I}}\) is the identity matrix, and u is the large-scale LIA’ VLM process, which we write as:
$${{{\bf{u}}}} \sim {{{\mathcal{N}}}}(\alpha {{{\bf{1}}}}+{\mathsf{G}}{{{\bf{c}}}},\Omega )$$
(5)
where \({\mathsf{G}}\) is the [K × N] matrix of K = 14,040 GIA model solutions at the N locations, \({{{\bf{c}}}}={[{c}_{1},…,{c}_{k}]}^{T}\) is an unknown [K × 1] selection vector that picks out one of the GIA model solutions, and α1 and Ω are, respectively, the unknown mean vector and covariance matrix of the regional VLM process unrelated to GIA, such that:
$${\Omega }_{ij}={\omega }^{2}\exp \left(-\rho | {s}_{i}-{s}_{j}| \right)$$
(6)
where ω2 is an unknown variance, ρ is an unknown inverse length scale, and ∣si − sj∣ is the distance between the ith and jth locations; this structure assumes that the covariance between two locations decreases with increasing distance. We model c as a multinomial random variable:
$${{{\bf{c}}}} \sim {{{\mathcal{M}}}}({{{\boldsymbol{\pi }}}})={\prod }_{k=1}^{K}{\pi }_{k}^{{c}_{k}}$$
(7)
where \({{{\mathcal{M}}}}({{{\boldsymbol{\pi }}}})\) is the multinomial distribution with probabilities \({{{\boldsymbol{\pi }}}}={[{\pi }_{1},\ldots,{\pi }_{K}]}^{T}\). A helpful analogy is to imagine c as a (potentially weighted) K-sided die and π as the probabilities of rolling the respective sides of the die.
Parameter level
To complete the model, we place priors on the unknown parameters. We choose weak, uninformative priors so that posterior solutions are controlled more by the available data than prior beliefs encoded into the algorithm. We use a normal prior for α:
$$\alpha \sim {{{\mathcal{N}}}}(\widetilde{{\eta }_{\alpha }},\widetilde{{\zeta }_{\alpha }^{2}}),\,\widetilde{{\eta }_{\alpha }}=3.7\,{{\mathrm{mm}}}\,{{{\mathrm{yr}}}}^{-1},\,\widetilde{{\zeta }_{\alpha }^{2}}={\left(11.5\,{{\mathrm{mm}}}\,{{{\mathrm{yr}}}}^{-1}\right)}^{2}$$
(8)
where \(\widetilde{{\eta }_{\alpha }}\) is the prior mean and \(\widetilde{{\zeta }_{\alpha }^{2}}\) is the prior variance. Note that we use superscripted tildes to identify hyperparameters (that is, parameters on the prior distributions). We use inverse-gamma priors for ϵ2 and ω2:
$${\epsilon }^{2} \sim {{{{\mathcal{G}}}}}^{-1}(\widetilde{{\xi }_{\epsilon }},\widetilde{{\chi }_{\epsilon }}),\,\widetilde{{\xi }_{\epsilon }}=0.5,\,\widetilde{{\chi }_{\epsilon }}={\left(0.7\,{{\mathrm{mm}}}\,{{{\mathrm{yr}}}}^{-1}\right)}^{2}$$
(9)
$${\omega }^{2} \sim {{{{\mathcal{G}}}}}^{-1}(\widetilde{{\xi }_{\omega }},\widetilde{{\chi }_{\omega }}),\,\widetilde{{\xi }_{\omega }}=0.5,\,\widetilde{{\chi }_{\omega }}={\left(1.6\,{{\mathrm{mm}}}\,{{{\mathrm{yr}}}}^{-1}\right)}^{2}$$
(10)
where \(\widetilde{\xi }\) is the shape parameter and \(\widetilde{\chi }\) is the scale parameter of inverse-gamma distribution \({{{{\mathcal{G}}}}}^{-1}\). For ρ, we use a log-normal prior:
$$\log \rho \sim {{{\mathcal{N}}}}(\widetilde{{\eta }_{\rho }},\widetilde{{\zeta }_{\rho }^{2}}),\,\widetilde{{\eta }_{\rho }}=-6.9\,\log {{{\mathrm{km}}}}^{-1},\,\widetilde{{\zeta }_{\rho }^{2}}={\left(0.35\,\log {{{\mathrm{km}}}}^{-1}\right)}^{2}$$
(11)
where \(\widetilde{{\eta }_{\rho }}\) is the “mean” and \(\widetilde{{\zeta }_{\rho }^{2}}\) is the “variance” of the distribution. Finally, we use a symmetric Dirichlet prior for π:
$${{{\boldsymbol{\pi }}}} \sim {{{\mathcal{D}}}}(\widetilde{\mu })=\frac{\Gamma \left(\widetilde{\mu }K\right)}{\Gamma {\left(\widetilde{\mu }\right)}^{K}}{\prod }_{k=1}^{K}{\pi }_{K}^{\widetilde{\mu }-1}$$
(12)
where \({{{\mathcal{D}}}}(\widetilde{\mu })\) is the Dirichlet distribution with concentration parameter \(\widetilde{\mu }={K}^{-1}\) and Γ is the gamma function. Note that this prior on π gives equal weight to all GIA model predictions. In keeping with our earlier analogy, this is like starting with a fair (unweighted) die, which has uniform probability of rolling any side. However, through the data assimilation, our posterior inference on π corresponds to an unfair (weighted) die, which favors GIA solutions that correspond better to the observations.
Posterior distribution
Given Bayes’ rule, we assume the posterior distribution is:
$$ p({{{\bf{v}}}},{{{\bf{u}}}},{{{\bf{c}}}},{{{\boldsymbol{\pi }}}},\alpha,{\omega }^{2},\rho,{\epsilon }^{2}| {{{\bf{x}}}}) \propto \\ p(\alpha )p({\omega }^{2})p(\rho )p({\epsilon }^{2})p({{{\boldsymbol{\pi }}}})p({{{\bf{c}}}}| {{{\boldsymbol{\pi }}}})p({{{\bf{u}}}}| \alpha,{\omega }^{2},\rho,{{{\bf{c}}}})p({{{\bf{v}}}}| {{{\bf{u}}}},{\epsilon }^{2})p({{{\bf{x}}}}| {{{\bf{v}}}})$$
(13)
Here, we use the symbols p, ∣, and ∝ to indicate probability distribution, conditionality, and proportionality, respectively. We draw samples from the posterior solution using standard numerical methods89. For u, v, and most parameters, we draw samples using a Gibbs sampler, whereas for ρ, which has a nonstandard conditional posterior distribution, we use a Metropolis sampler. We perform 200,000 iterations with the sampler. We omit the first 100,000 iterations, which we regard as burn in, to eliminate startup transients and the effects of initial conditions. To reduce serial correlation, we thin the remaining iterations, only keeping one out of every 100 remaining samples. We analyze the remaining 1000 joint posterior samples. We also computed standard model diagnostics89 to validate the model and interrogate its solutions (Supplementary Section 3, Supplementary Figs. 6–9). The Bayesian model code was written and executed using the MATLAB numeric computing platform and takes ~20 min to run on a 2023 MacBook Pro (Apple M2 Pro Chip; 32 GB Memory). The final results of our Bayesian analysis are posterior distributions on the GIA model parameters shown in Table 1 in the form of the posterior vector c that is then used to inform our prediction of present-day VLM driven by the LIA (Fig. 3) and our projections of LIA-to-2100 CE RSL change (next section).
Projections of RSL change driven by LIA-to-2100 CE ice-sheet change
Projected mass changes from Greenland and Antarctica to 2100 CE are adopted from eight ice-sheet modeling groups working under the framework of the Ice Sheet Model Intercomparison Project for CMIP6[ISMIP6; 22,44. For each modeling group, we use projections based on future emissions scenarios RCP 2.6 and RCP 8.5. The projections are converted to ice thickness anomaly relative to 2017 CE and concatenated to the Greenland LIA ice history used to calculate present-day VLM caused by the LIA. We ran a total of 224,640 (14,040 LIA ice histories/Earth structures × 8 ISMIP6 modeling groups × 2 RCP scenarios) 1D Maxwell GIA simulations, where we vary the mantle viscosity, lithospheric thickness, LIA mass anomaly, and LIA termination time over the ranges listed in Table 1. Each model is run at spherical harmonic degree 512 for a total of 26 timesteps from 1000 CE to 2100 CE, with timesteps from 2020–2100 CE increasing to 5-year intervals. A control simulation, in which the GIA model is forced by future ice-mass anomalies projected for a fixed modern climate, is subtracted from the projections of RSL to correct for model drift that arises from the initialization strategies of the individual ice-sheet modeling groups22. This was done for each combination of GIA model parameters (i.e. required an additional 112,320 GIA simulations). To obtain the median and likely range of future RSL change across our suite of GIA simulations, we use the insight gained from our Bayesian inversion and multiply our modeled RSL values at each timestep with the posterior vector c (see previous section). We then calculate the median and likely range across the resulting posterior distribution.
Contributions to projected RSL change excluding ice-sheet mass changes and associated GIA
In addition to RSL change caused by variations in ice-sheet mass and related GIA, other processes will contribute to the RSL change experienced by coastal Greenland over this century. These processes include changes in global glaciers, thermal expansion and ocean dynamics, and changes in terrestrial water storage, excluding glaciers and ice sheets1. We refer to the combination of these as the “other processes” contribution. Note that the contribution from future ice change in Antarctica is already included in the previous section. Estimates for the RSL contribution of each of these processes to 2100 CE for SSP1-2.6 (Shared Socioeconomic Pathway 1 paired with an approximate radiative forcing of 2.6 W/m2 by 2100 CE) and SSP5-8.5 (Shared Socioeconomic Pathway 5 paired with an approximate radiative forcing of 8.5 W/m2 by 2100 CE) were obtained from the dataset of RSL projections associated with the IPCC AR6 that excludes the signal driven by long-term vertical land movement1,46,47,48. The projections of RSL provided by this dataset are given as percentiles of a probability distribution at 10-year timesteps from 2020 to 2100 CE relative to a 1995–2014 CE baseline. We linearly extrapolate these percentiles to 2017 CE and make predictions from 2020–2100 CE relative to 2017 CE (instead of 1995–2014 CE). For each process, we construct a probability distribution at each timestep based on the percentiles. We then draw 1000 samples for each process and each timestep and combine them to obtain a distribution for the other processes component. This allows us to calculate the median and likely range for this component (Fig. 4i–l). For our total estimate of RSL change around Greenland (Fig. 5), we combine the 1000 random draws from the other processes component with random draws from the posterior distribution of future RSL driven by LIA-2100 CE ice-sheet change and the contribution from deglacial ice change. For the latter we sample predictions according to the deglacial model weights. We then calculate the median and likely range of total projected RSL across these combined samples.