Projections of 21st-century sea-level fall along coastal Greenland

Domain: earth_science
Published: 2026-01-20
Editorial outcome: major_revision

Abstract
The Greenland Ice Sheet will significantly contribute to global mean sea-level rise this century. However, glacial isostatic adjustment is expected to cause regional sea-level fall around Greenland as the land rebounds and the gravitational pull of the shrinking ice sheet decreases. A fall in local sea level has implications for Greenlandic communities whose economy, near-shore infrastructure, and food security are vulnerable to coastal changes, and for the dynamics of marine-terminating glaciers. In this study, we use a glacial isostatic model paired with sea-level and vertical land motion observations to predict future sea-level change in this region. We find that local sea level will fall around the island by 2100 CE (relative to 2017 CE) reaching a median value and likely range (17th-83rd percentile) of −0.9 m [−1.2 m to −0.66 m] and −2.5 m [−3.8 m to −1.7 m] for Representative Concentration Pathways 2.6 and 8.5, respectively, with viscous effects contributing 20–40% of the predicted local sea-level signal associated with future ice change.

Introduction
As anthropogenic climate change continues, both the thermal expansion of ocean water and the melting of ice sheets and glaciers will contribute to a rise in global mean sea level (GMSL). The Greenland Ice Sheet (GrIS) is expected to be a primary contributor to rising sea levels, projected to cause a median value and likely range (17th-83rd percentile) of GMSL rise by 2100 CE relative to 1995–2014 CE of 60 mm [10 mm to 100 mm] for Shared Socioeconomic Pathway (SSP) 1–2.6 and 130 mm [90 mm to 180 mm] for SSP 5-8.5 1. However, little research has been done to explore how the melting ice sheet will change sea level around the coastlines of Greenland itself. Since Greenland has nearly 60,000 permanent inhabitants, all of whom live in towns and settlements along the coast, changing sea levels will impact the majority of Greenland’s population. Sea-level change around Greenland deviates substantially from the global mean mainly due to glacial isostatic adjustment (GIA), which is the deformational response of the solid Earth and associated modulation in Earth’s gravitational field and rotational axis to changes in ice and ocean mass through time 2,3. In the near field of a melting ice sheet, GIA causes a local sea-level fall, as the reduction in ice mass causes the Earth beneath it to rebound and a decrease in the gravitational attraction of the ice sheet causes the sea surface to lower 4. This contrast in the sign of expected local sea-level change from that of the majority of other vulnerable coastlines around the world motivates the need for region-specific projections of sea-level change for Greenland. In this study, we will use the term relative sea level (RSL) to denote local sea surface height relative to bedrock topography 5.

In order to predict future changes in sea level due to GIA around Greenland, we need to know past and future ice-mass change, as well as the rheology of Earth’s interior, which determines the magnitude and timescale of the solid Earth response. Regarding past ice-sheet changes, the GrIS peaked in extent around 16.5 kyr ago when it held an excess volume relative to present of approximately 5 m ice-equivalent sea level and extended over large swaths of the continental shelf 6. While the timing of its retreat from this position is believed to have varied spatially across the ice sheet 6, the majority of the GrIS receded to or in land of the present-day shoreline by 10 ka 6–9. The early-to-mid Holocene was marked by higher-than-present temperatures in the Northern Hemisphere 10,11 and much of the GrIS and its peripheral glaciers retreated past their contemporary extent 11–14. This was followed by the readvance of the GrIS past its present position during the Neoglacial cooling period, which began approximately 5000 yrs ago and culminated in the Little Ice Age (LIA) around 800–120 yrs ago 15. The mass of the GrIS and its peripheral glaciers is estimated to have been 14,862 ± 3758 Gt larger than present during the LIA 16–18. The onset of ice retreat from its LIA maximum extent varies spatially 17, but is generally centered around 1850–1900 CE 16. More recently, measurements from satellite altimetry reveal that the GrIS decreased in surface elevation along the majority of its margins, excluding those on the central east coast and a small area in the southeast, during the period 2012–2017 CE 19. Today, the rate of mass loss from the GrIS and its peripheral glaciers is 262 ± 39 Gt/yr 20 and 42 ± 6 Gt/yr 21, respectively. Though projections to 2100 CE show decreases in ice thickness along all margins of the GrIS, the southwestern margin is especially susceptible to change, partially reflecting the pattern of surface elevation change seen today 22. Accurately incorporating the evolution of the GrIS from the Last Glacial Maximum (LGM) to 2100 CE is fundamental to reliably modeling the viscoelastic response of the Earth and the reshaping of the ocean’s surface to these mass changes.

Knowledge of Earth’s internal structure, especially mantle rheology, is equally important for accurately capturing the GIA response to changing surface loads. In GIA studies, Earth’s mantle is often considered to be a radially-symmetric (i.e., 1D) viscoelastic Maxwell body 4. Solid Earth parameters are typically determined by tuning GIA models to fit paleo RSL data 6,23–25, modern GPS data 16,26,27, and Earth’s rotation parameters 28. Lecavalier et al. 6 find that GIA models forced by a GrIS deglaciation history (termed Huy3) derived from a dynamical ice-sheet model best fit RSL data in Greenland when using a 120 km thick elastic lithosphere and upper and lower mantle viscosities of 5 × 10^20 Pas and 2 × 10^21 Pas, respectively. However, these values are not consistent across all studies that tune their model parameters to fit Greenland RSL data. Tarasov & Peltier 24 find that their GrIS history performs well when coupled with a lithosphere of 90 km. Fleming & Lambeck 25 also show a preference for a thinner lithosphere (80 km), in addition to higher values for the viscosities of the upper mantle (4 × 10^21 Pas) and lower mantle (10 × 10^21 Pas). Each of these studies finds that their predictions of RSL are unable to match all of the observed RSL data, a consequence of uncertainties in both the ice reconstruction and prescribed Earth parameters. Regarding the latter, one assumption in all of these studies is that Earth’s properties vary with depth alone when evidence for a more complex, laterally-varying (i.e., 3D) Earth structure below Greenland is present [e.g. 29]. Milne et al. 30 predict Holocene RSL in Greenland using a 3D GIA model paired with Huy3 and find that this model can improve data-model residuals in some locations.

Contemporary observations of GIA include vertical land motion (VLM) measured with permanent Global Navigation Satellite Systems (GNSS). Present-day VLM recorded at the 57 permanent GNSS stations that span coastal Greenland, called GNET, is caused by deformation driven by ice changes that have occurred in the past (most importantly over the last deglaciation, i.e., last 26 kyr) and that are occurring today (i.e., during the recorded period) 31. These records are often corrected for the deformation caused by coeval ice change (assumed to be elastic), so that trends in the data can be directly compared to predictions of VLM driven by past ice changes from GIA models 16,20,26,32. However, these models, which are tuned to fit paleo RSL data, consistently underpredict VLM observed from GNET 16,20,26,30,32,33. One source of uncertainty in this comparison comes from the incomplete removal of the trend in elastic deformation. Berg et al. 20 consider both longer time series (~2007–2023) of GNSS data and mass changes from the GrIS, its peripheral glaciers, and glaciers in the Canadian Arctic to derive more accurate elastic VLM solutions that can improve constraints on GIA models. They find that discrepancies still remain between their elastic-corrected GNET data and predictions from multiple GIA models forced by deglacial ice change. Milne et al. 30 consider whether the influence of lateral variations in Earth structure could help explain poor fits to the present-day VLM observations, but conclude that the 3D viscosity structures considered in their study do not explain the large residuals.

Adhikari et al. 16 propose that this inconsistency is due to the exclusion of ongoing deformation caused by more recent ice-mass change that has occurred since the beginning of the LIA. They show that including a LIA ice history paired with a reduced mantle viscosity for the time period encompassing these mass changes allows for a better fit between their predicted present-day VLM and the observed values at the majority of GNET stations. They suggest that this lowering in mantle viscosity over the timeframe of ice-sheet changes since the beginning of the LIA hints at the presence of transient deformation beneath Greenland. Transient deformation is an extension of the classical Maxwell rheology. While the latter consists only of an instantaneous elastic response and a steady-state viscous creep, transient deformation, as is, for example, captured by a Burgers rheology, allows for additional dissipative mechanisms which act to lower the apparent viscosity over short periods. This deformation style is supported by high-temperature high-pressure laboratory experiments of rock deformation 34,35 and has been explored in other GIA studies e.g., 36–40. In the case of loading events that occur over centennial timescales (like the LIA), the viscosity of the mantle can effectively be lower than that inferred for longer-timescale events (like the LGM or mantle convection), allowing the mantle to deform more quickly and with a greater magnitude. Paxman et al. 41 use rheological laws paired with an inference of the mantle thermodynamic state from seismic tomography to show that one indeed would expect a lowering in apparent viscosity for ice-mass changes that occur over decadal-to-centennial timescales consistent with LIA mass loading. However, Pan et al. 42 caution that a sub-lithospheric low-viscosity layer with a regular Maxwell rheology would also be able to explain both deglacial RSL data and present-day VLM if the LIA ice history is included in the GIA model. Further work is needed to understand the possible extent of a sub-lithospheric low-viscosity layer around and below the cratonic keel of Greenland 29,43.

The goal of this study is to produce reliable projections of sea-level change around Greenland for this century. We build on previous work by accounting for both lateral variability in viscosity over long timescales and a reduced mantle viscosity over shorter timescales in our projections. We assimilate paleo RSL data from the Holocene as well as VLM data from GNET stations to reduce uncertainty in our knowledge of solid Earth deformation. Finally, we produce a complete estimate of future sea-level change around Greenland by projecting sea level using a GIA model forced with predictions of past and future ice-sheet behavior and combining it with additional sources of future sea-level change (e.g., global-ocean thermal expansion and melt of global glaciers) under a low and high emissions scenario.

In this work, we find that RSL will fall along Greenland’s coast, reaching a median value and likely range of −0.9 m [−1.2 m to −0.66 m] and −2.5 m [−3.8 m to −1.7 m] for Representative Concentration Pathways (RCP) 2.6 and 8.5, respectively, at 2100 CE relative to 2017 CE. The viscous contribution to projected RSL associated with future ice-sheet change is 20–40%, emphasizing the importance of modeling the full viscoelastic response to projected surface load changes over this century. A sea-level fall of this magnitude has implications for coastal communities and could affect ice-sheet dynamics at marine-based outlets.

Results
To make projections of 21st century sea-level change around Greenland we consider 3 contributions: (1) future sea-level change driven by historical ice-mass changes over the deglaciation (30 ka to 1 ka), (2) future sea-level change driven by historical and projected ice-mass changes from the LIA (beginning at 1000 CE) to 2100 CE, and (3) future sea-level change driven by processes unrelated to ice-sheet mass change and associated GIA (which we label the “other processes" contribution). For (1), we use a suite of 65 available 3D Maxwell GIA simulations (Supplementary Table 1) and weight them based on their fit to Holocene RSL data. For (2), we use a suite of 14,040 1D Maxwell GIA simulations. We first predict present-day VLM in Greenland driven by the LIA and compare these to observations of uplift measured at 55 GNET stations. In order to make a direct comparison, we correct the observed VLM for uplift associated with deglacial and coeval mass changes, the former of which is calculated in step (1) and the latter is taken from Berg et al. 20. We further account for potential differences between the reference frame used in the 1D GIA model and that used for the VLM observations (Supplementary Section 1). We perform a Bayesian inversion to identify the distribution of most likely ice and Earth parameters varied in the 1D GIA model and create a posterior prediction of present-day VLM driven by the LIA in Greenland. We then use these optimal parameters to create a suite of 1D Maxwell GIA runs forced by LIA-to-2100 CE ice-mass change, using an ensemble of projections of Greenland and Antarctic ice change for both RCP 2.6 and RCP 8.5 22,44. These represent low and high scenarios for radiative forcing to 2100 CE 45, which are considered in the Intergovernmental Panel on Climate Change (IPCC) Assessment Reports. Modeling contributions (1) and (2) separately allows for the consideration of different viscosities over different timescales to approximate a transient response. An assumption of 1D viscosity over shorter timescales (contribution 2) is made to improve computational efficiency and because the behavior of the mantle is inferred to be more laterally homogeneous over decadal-to-centennial timescales than over multi-millennial timescales 41. For (3), we use published estimates 1,46–48 associated with the IPCC AR6 for the contribution to sea-level change around Greenland from thermal expansion, ocean dynamics, land water storage, and global glaciers. Projections of RSL from these contributions are provided for SSP 1–2.6 and SSP 5-8.5, which combine socioeconomic narratives defined by the SSPs with emission pathways based on the RCPs.

Present-day VLM from deglacial mass change
Our calculations indicate that GIA from deglacial ice loss causes present-day uplift of up to 3.3 mm yr−1 along the northern and eastern coasts of Greenland and subsidence at rates less than 2.5 mm yr−1 in the south and west (Fig. 1). Our modeled rates of present-day VLM caused by the deglaciation are in general agreement with previous estimates 6,27,30,31. The spatial pattern is governed by multiple processes. First, uplift in response to overall ice-sheet retreat since the LGM dominates the northeastern margin, which began to retreat at 12 ka before remaining relatively stable around its present-day position since the mid-Holocene 6,24. Second, subsidence in response to late-Holocene regrowth of the GrIS during the Neoglacial period contributes to negative VLM along the southwest coast 6,49. And, third, subsidence due to Greenland’s position on the peripheral bulge of the Laurentide Ice Sheet results in larger subsidence rates along the western coast of Greenland compared to the east. Like previous studies 16,20,26,30,32, we find that the rates of VLM caused by deglacial ice loss alone are insufficient to fit the elastic-corrected observations at the majority of GNET stations. The largest discrepancy is seen along the southeastern coast, where the residual between the elastic-corrected observed and the modeled deglacial VLM rates reaches a maximum of 6.2 mm yr−1.

The weights for the deglacial GIA simulations (normalized so that they sum to one) fall within the range of 0.0081 to 0.0227, and do not deviate substantially from the value that would be assigned if all predictions were weighted equally (1/65 = 0.0154). Our highest-weighted simulation is that from Antwerpen et al. 49, which utilizes the Huy3 6 ice history and a combination of 2 seismic tomography models 50,51 to infer the 3D viscosity structure of the mantle. Comparing the highest-weighted simulation with the lowest-weighted simulation from Austermann et al. 52, in which the Earth structure is the same as that used in Antwerpen et al. 49 but the ice history is given by ICE-6G 53 instead of Huy3, suggests that the Huy3 ice reconstruction (rather than the Earth structure) is the key feature of the best-fitting prediction. This is expected, as additional and improved data constraints were used in the creation of Huy3 compared to the GrIS portion of ICE-6G.

Bayesian inversion of present-day VLM from LIA mass change
Incorporating VLM from LIA mass changes has been shown to reduce the residuals between modeled VLM and observed VLM within the period 2011–2017 CE 16. In this study, the contribution of LIA mass change to present-day VLM is determined through assimilating VLM observations into our modeling framework. Figure 2 shows this quantity’s prior estimate (blue markers, note that this shows LIA’, which includes the correction for a difference in the reference frame between model predictions and observations), which is the unweighted average taken across our entire suite of GIA model predictions, and posterior estimate (red markers, note this also shows LIA’), which is the weighted average taken after performing the Bayesian inversion. There is only a slight difference in mean between our prior and posterior solutions, as the model parameters used in this study (Table 1) were sampled from the posterior distribution of parameters from Adhikari et al. 16. Nonetheless, the posterior solutions have decreased variance at all stations and the predicted mean is closer to the observed mean at 34 out of the 55 stations when compared to the prior, which demonstrates learning from the data.

The posterior prediction of the LIA contribution to present-day VLM (LIA-VLM) shows positive uplift rates everywhere in Greenland, with the highest magnitude rates of up to 5.7 mm yr−1 seen along the southeast coast (Fig. 3). The standard deviation of the posterior solution has a similar pattern to the posterior mean and is 20–30% of its magnitude (averaged across GNET stations). Regions with outlet glacier systems whose mass has decreased significantly since the LIA, specifically the Kangerlussuaq Glacier and Helheim Glacier (station 23, Fig. 2) in the southeast 15–17, show higher magnitudes of uplift than other regions along the ice sheet margin. The average reference frame contribution to the observed uplift rates (Supplementary Fig. 1) is between −1.8 mm yr−1 and 0.45 mm yr−1, with the maximum magnitude contribution seen along the southeast coast. The posterior LIA-VLM rates calculated in this study are of similar magnitude and spatial pattern to those predicted by Adhikari et al. 16, suggesting that the prediction of present-day LIA-VLM is not heavily influenced by our use of different observed, elastic, and deglacial VLM rates or the difference in our methods for weighting the LIA-GIA model predictions.

Our solutions of LIA-VLM are sensitive to both mantle viscosity and lithospheric thickness and insensitive to LIA termination time and mass anomaly. Preferred values for mantle viscosity and lithospheric thickness are 6.3–7.0 × 10^19 Pas (probability = 0.33) and 200 km (probability = 0.17), respectively. The covariance between these parameters shows the expected trade-off where a thicker lithosphere requires a lower mantle viscosity (Supplementary Fig. 2). Our preferred mantle viscosity is significantly smaller than viscosities found by fitting deglacial data 6,9,25. However, it is in agreement with the inferred apparent viscosity of the mantle beneath Greenland for centennial timescales determined by the use of full-spectrum rheological models that examine the influence of transient deformation 41. Following this work and the interpretation of Adhikari et al. 16, our LIA viscosity estimate may be indicative of transient deformation in the mantle beneath Greenland. Our preferred lithospheric thickness is larger than the best-fitting estimates found by reducing the misfit of deglacial GIA models to observations 9,25, however, a higher apparent lithospheric thickness is expected for shorter LIA timescales 41.

Adhikari et al. 16 find a preferred viscosity of 6–11 × 10^19 Pas with little sensitivity to lithospheric thickness. This differs slightly from our results, which may be due to different choices in the data and their processing (between Adhikari et al. 16 and Berg et al. 20) or due to differences in the inference framework. To test these two options, we repeat our analysis with the observed and elastic VLM rates from Adhikari et al. 16 instead of Berg et al. 20. We find that this produces a preferred viscosity of 6–7 × 10^19 Pas (in line with our original inference) and a preferred lithospheric thickness of 60–80 km (different to both our original inference and the result by Adhikari et al. 16). This shows that both choices in the data and the inference scheme lead to differences in the inferred Earth structure.

Projections of RSL change
Projections of RSL change along coastal Greenland to 2100 CE are a combination of changes driven by deglacial GIA (Fig. 4a, b), LIA GIA (Fig. 4c, d), GIA associated with future ice-sheet change (Fig. 4e, h), and processes unrelated to ice-sheet change and associated GIA (i.e., the “other processes"; Fig. 4i–l). All projections reported in this study are given as the median and likely range (17th–83rd percentile) and are relative to 2017 CE. The deglacial contribution causes a slight rise in RSL, reaching 0.21 m [0.17 m to 0.26 m], in west and southeast Greenland due mainly to Laurentide forebulge collapse and a fall in RSL, reaching −0.17 m [−0.46 m to −0.10 m], in the northeast caused by ongoing uplift in response to ice loss since the LGM (Fig. 4a). The LIA contributes a RSL fall along the entire coastline and is a significant influence in the northwest, where it reaches −0.33 m [−0.40 m to −0.24 m], and in the southeast, where it reaches −0.41 m [−0.50 m to −0.30 m] (Fig. 4c). The contribution from future ice-sheet mass change is calculated by pairing an ensemble of ice-sheet projections for RCP 2.6 and RCP 8.5 22,44 with preferred parameters from our LIA Bayesian inversion (see “Methods”). This contribution produces falling sea levels everywhere around Greenland (driven by mass loss from the GrIS), with the largest changes, reaching −0.88 m [−1.1 m to −0.73 m] for RCP 2.6 and −2.5 m [−3.9 m to −1.9 m] for RCP 8.5, occurring adjacent to the GrIS margin where melting is projected to be greatest (Fig. 4e, g).

We performed an additional calculation to quantify the relative influence of solid Earth deformation versus gravitational effects on our projection of RSL change driven by LIA-2100 CE ice load variations. Gravitational effects include both the ice and the solid Earth’s pull on the ocean’s surface and a uniform gravitational shift caused by the input of meltwater into the ocean. We calculate this by running the GIA code from the LIA to 2100 CE using the ensemble of ice-sheet projections for RCP 2.6 and RCP 8.5 with one parameter combination that performs well in fitting observations (mantle viscosity = 7.9 × 10^19 Pas, lithospheric thickness = 200 km, LIA termination time = 1870 CE, LIA mass anomaly = 16,500 Gt) and separating the two components. We find that RSL at 2100 CE caused by ice-mass changes from the LIA to 2100 CE is dominantly driven by solid Earth deformation, with sea-level changes caused by gravitational effects constituting ~0–20% of the total for RCP 2.6 and ~10–30% of the total for RCP 8.5, depending on location around the coast (Supplementary Fig. 3).

The contribution to RSL at 2100 CE from thermal expansion and ocean dynamics, global glacier change, and land water storage 1,46–48 is positive along Greenland’s coasts except for a small area of RSL fall in the northwest (Fig. 4i, k). Thermal expansion and ocean dynamics are the primary drivers of RSL rise along the coast, causing up to 0.28 m [0.15 m to 0.42 m] and 0.62 m [0.36 m to 0.87 m] at 2100 CE for SSP 1–2.6 and SSP 5–8.5, respectively. The contribution of land water storage to RSL change in Greenland is also positive, but insignificant and does not reach above 0.03 m [0.02 m to 0.04 m] for either scenario. The fall in RSL in northwest Greenland is caused by projected ice loss from the peripheral glaciers of Arctic Canada and their resulting GIA signal, which was calculated using an elastic Earth structure 1,54. It reaches −0.38 m [−0.65 m to −0.10 m] for SSP 1-2.6 and −1.0 m [−1.3 m to −0.70 m] for SSP 5-8.5. Calculating the viscoelastic GIA response to future changes in mountain glaciers around Greenland and Arctic Canada using our preferred Earth structure for sub-centennial timescales may significantly increase the magnitude of this component.

Summing all contributions leads to a RSL change at 2100 CE (Fig. 5a, b) that, in pattern, resembles the contribution from GIA associated with the ice-sheet projections (Fig. 4e–h). At 2100 CE, the contributions to RSL change from GIA driven by past ice-mass variations and from the other processes have opposing signals around the coastline and tend to cancel out. The total projected fall in sea level reaches a maximum value of −0.90 m [−1.2 m to −0.66 m] near Sermeq Kujalleq (formerly known as Jakobshavn Isbrae) for RCP 2.6 and −2.5 m [−3.8 m to −1.7 m] near Helheim Glacier for RCP 8.5. Note that here and for the remainder of the text, we refer to the RSL projections as representing emissions scenarios RCP 2.6 and 8.5, which include the SSP 1-2.6 and SSP 5-8.5 projection components for the other processes.

The spatial pattern of the likely range of total RSL at 2100 CE demonstrates that the largest source of uncertainty in the projections comes from the different ways in which the future ice models 22 project marine-terminating glacier response to climate warming (Fig. 5c, d). This is illustrated further in Supplementary Fig. 4, which shows projected RSL change at 2100 CE under RCP 8.5 for each of the individual ice-sheet models considered in this study.

Time series of RSL from 2017–2100 CE for RCP 8.5 in Tasiilaq and Nuuk demonstrate that each contribution to RSL, and therefore the total RSL change, will vary for different settlements in Greenland (Fig. 6). In Tasiilaq, the most populous settlement along the eastern coast, made up of nearly 2000 inhabitants, RSL is projected to fall by −1.4 m [−1.8 m to −0.95 m] for RCP 8.5. Projected ice-sheet mass loss contributes −1.6 m [−1.9 m to −1.2 m] to this RSL fall, which is further exacerbated by the RSL fall driven by LIA mass loss and partially offset by the RSL rise from other processes. Nuuk, the capital of Greenland and the most populated area on the island with 19,611 inhabitants, is projected to experience a RSL fall of −0.83 m [−1.2 m to −0.54 m] by 2100 CE for RCP 8.5. Compared to Tasiilaq, the RSL fall driven by projected ice change and the LIA is smaller in Nuuk. This is combined with a larger RSL rise from deglacial ice change and other processes, leading to the smaller total RSL fall seen in Nuuk over this century. Though the magnitude of RSL change and the breakdown of its contributions will vary along Greenland’s coasts, RSL is projected to fall for both RCP 2.6 and RCP 8.5 at all of the six selected settlements discussed in this study (shown with the white stars in Fig. 5c and discussed further in Supplementary Section 2). For all locations considered here, it is noteworthy that projections for the two future emissions scenarios (RCP 2.6 vs. 8.5) only start to diverge around 2050 CE, which is consistent with prior work on sea-level projections globally [e.g. 1,55,56] (Supplementary Fig. 5).

The most recent Assessment Report (AR6) of the IPCC includes regional projections of sea-level change at 2100 CE relative to a 1995–2014 CE baseline 1. These differ from our Greenlandic projections in both magnitude and sign. For SSP 1-2.6, IPCC AR6 projects a RSL rise on the order of half a meter or less along most of Greenland’s coasts, except the northern and northeastern portions, where a RSL fall is projected. For SSP 5-8.5, projected RSL change is similar to that for SSP 1-2.6 in most areas, but RSL fall is now projected to occur in the northwest near Kullorsuaq and is stronger along the northern coast. Our projections differ from those in AR6 in three ways: (1) how GIA in response to past ice melt is calculated, (2) how GIA in response to future ice melt is calculated, and (3) the ice-sheet models that are used for future projections. Regarding the first point, the authors of AR6 estimate RSL caused by VLM from GIA driven by historical ice changes and additional processes with a Gaussian process regression method 57 that utilizes tide-gauge data. This dataset does not include tide-gauges in Greenland, making their inference for the long-term VLM component of RSL change in Greenland uncertain. The authors of AR6 recognize the limitations in this method and assign their projections of long-term VLM and GIA low-to-medium confidence, stating that for many regions, a regional analysis would result in projections of greater confidence 1. Regarding point (2), the authors of AR6 use a sea-level equation solver 54 with an elastic Earth structure. Viscous deformation associated with future ice change, which we estimate to be significant (see next section), is not included. Both factors likely cause AR6 projections to under-predict the amount of RSL fall in Greenland over the next century. Lastly, regarding point (3), we highlight that our RSL projections were made using future ice-sheet mass changes from ISMIP6 simulations that were forced by CMIP5 (and not CMIP6) climate output. The range of climate sensitivities represented in the ensemble of CMIP6 models has an upper tail that exceeds the range of climate sensitivities from the CMIP5 models used to create the ice-sheet projections considered in this study 58. Therefore, the ice-sheet projections used here represent a more conservative-to-moderate subset of those considered in AR6.

Implications of reduced mantle viscosity for sea-level projections
The ramifications of considering a lower mantle viscosity applied to sub-centennial timescales for projections of sea-level change around Greenland are quantified in Fig. 7, which shows the elastic and viscous contributions to RSL change at 2100 CE driven by projected ice-sheet mass loss. If the viscosity of the solid Earth is taken to be the steady-state viscosity used for long timescale loading events, it is common to calculate GIA associated with projections of ice-mass change using an elastic Earth structure [e.g., 1,54,59], as the viscosity is too high and timescale too short for there to be substantial viscous changes. However, if a shorter period forcing elicits a lowered mantle viscosity, the timescale in which the viscous response can become important is much shorter. Although the elastic contribution (Fig. 7a) is still dominant, the viscous signal (Fig. 7b) is responsible for 20–40% of the projected RSL change caused by future ice-sheet mass loss (Fig. 7c). Excluding this viscous response would lead to a significant under prediction of RSL fall around Greenland. This highlights the importance of calculating the full viscoelastic response of the solid Earth in response to future surface loading when projecting RSL, and the need for a continued focus on how to most appropriately model the seemingly time-dependent nature of mantle relaxation.

The fall in RSL by up to 2.5 m for RCP 8.5 that we project to occur along Greenland’s coasts by the end of this century could slow down the pace of grounding-line retreat for Greenland’s marine-terminating glaciers, as the position of the grounding line depends on sea level at the grounding line 60–62. Ice discharge from marine-terminating glaciers was responsible for approximately two-thirds of Greenland’s total contribution to global sea level over the period 1972–2018 CE 63,64. The response of the solid Earth and gravity field to this ice-mass loss can raise the bed topography and decrease RSL at the grounding line, which can lead to its stabilization and potential re-advance 61,62. This effect may be particularly impactful for large marine-terminating glaciers that rest on bedrock that deepens inland (like Kangerlussuaq and Sermeq Kujalleq 65) and are vulnerable to marine ice sheet instability feedbacks, where the flux of ice across the grounding line increases as it retreats 60,66–70. Adhikari et al. 67 find that uplift of 1–4 m by 2100 CE is sufficient to noticeably influence ice dynamics and promote the stabilization of marine segments of the Antarctic Ice Sheet. This amount of uplift is of similar magnitude to that projected to occur along Greenland’s coasts over a similar timeframe; however, a direct comparison between the dynamics of Greenland and Antarctica’s marine-terminating glaciers is difficult given differences in buttressing and other factors. This motivates the need for a detailed analysis that couples GIA and ice-sheet models on centennial timescales in Greenland.

Furthermore, a RSL fall of the magnitude projected in this study can impact the coastal communities of Greenland. RSL fall on the order of multiple meters has the potential to increase navigational hazards and limit the accessibility of harbors to deep draft boats by reducing the water depth beneath their keels. Marine channels that are already shallow today are vulnerable to becoming exposed at low tide, disconnecting navigational routes and water flow, which may impact marine habitats. Port infrastructure, such as docks and mooring sites, could become less useful as the falling sea level would make access to them more difficult or even impossible. A reliance on global sea-level projections can lead to the perspective that sea level will rise in Greenland due to climate change 71, which would create a very different set of hazards than those presented above. The regional RSL projections provided in this study can be used by communities to assess their site-specific risks to sea-level fall and prepare adaptation strategies accordingly. Subsequent interdisciplinary studies that consider these projections with regard to the national economy, coastal structure, and marine ecology are needed to better understand and communicate the impact that RSL fall could have on these communities.

Methods
Present-day VLM and projected RSL change driven by deglacial mass change
Ice changes over the deglaciation (30 ka-1 ka) 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-6G 53 and ANU 73, four seismic tomography models, which inform lateral variability in viscosity 50,51,74,75, two models for variations in lithospheric thickness 76,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 history 72. 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 & Lebedev 50 seismic tomography model for the mantle above the transition zone (<400 km depth) and the French & Romanowicz 51 model for the mantle below that. The Earth models explored 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 Huy3 6. 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 associated with the index and limiting data, according to:

WRSS nm = < (2r ε nm)^2 (0.5+0.5erf(q n r ε nm)), q n ≠ 0 (1)

where erf represents the error function, r nm 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 q n 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 WRSS nm, please refer to Creel et al. 79 and references within. The weight for each model is then calculated as the normalized inverse of WRSS nm 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 (1 ka–0 ka) 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 shorelines 3. 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 GrIS 17 and its peripheral glaciers 18 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 timeframe 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 × 10^19 Pas to 25 × 10^19 Pas (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 parameter 16. 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 Greenland 41, 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 mantle 42. 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 plume 16,26,82,83, which we don’t allow for in our calculations here. Observed uplift rates at each station (V OBS) are assumed to be the sum of VLM contributions due to ice-mass changes during the deglaciation (V DEG), ice-mass changes since the beginning of the LIA (V LIA), ice-mass changes during the GPS record (V ELA) and a reference frame correction (V RF; Supplementary Section 1):

V OBS = V DEG + V LIA + V ELA + V RF (2)

For the contribution of elastic deformation to present-day uplift rates, V ELA, 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 GrIS 84,85, Greenland’s peripheral glaciers 21, and Canada’s peripheral glaciers 20. 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 greatest 20. 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 Berliner 87, 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 & Wikle 88. Our posterior distributions are updates to those found in Adhikari et al. 16 since we used improved GNET data and estimated elastic contributions 20, 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 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 v = [v 1, ..., v N]^T according to:

x ~ N(v, Δ) (3)

where ~ means “is distributed as", N(a, b) is the multivariate normal distribution with mean vector a and covariance matrix b, and Δ = diag(δ^2 1, ..., δ^2 N) is the diagonal matrix of the data error variances where we define δ^2 i 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:

v ~ N(u, ϵ^2 I) (4)

where ϵ^2 is the unknown spatial variance of the local VLM process unrelated to GIA, I is the identity matrix, and u is the large-scale LIA’ VLM process, which we write as:

u ~ N(α1 + Gc, Ω) (5)

where G is the [K × N] matrix of K = 14,040 GIA model solutions at the N locations, 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:

Ω ij = ω^2 exp(−ρ|s i − s j|) (6)

where ω^2 is an unknown variance, ρ is an unknown inverse length scale, and |s i − s j| 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:

c ~ M(π) = πc k (7)

where M(π) is the multinomial distribution with probabilities π = [π 1, ..., π 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 α:

α ~ N(ηf α, ζf^2 α), ηf α = 3.7 mm yr−1, ζf^2 α = (11.5 mm yr−1)^2 (8)

where ηf α is the prior mean and ζf^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:

ϵ^2 ~ G^−1(ξe ϵ, χe ϵ), ξe ϵ = 0.5, χe ϵ = (0.7 mm yr−1)^2 (9)

ω^2 ~ G^−1(ξf ω, χf ω), ξf ω = 0.5, χf ω = (1.6 mm yr−1)^2 (10)

where ξ is the shape parameter and χ is the scale parameter of inverse-gamma distribution G^−1. For ρ, we use a log-normal prior:

log ρ ~ N(ηf ρ, ζf^2 ρ), ηf ρ = −6.9 log km−1, ζf^2 ρ = (0.35 log km−1)^2 (11)

where ηf ρ is the “mean" and ζf^2 ρ is the “variance" of the distribution. Finally, we use a symmetric Dirichlet prior for π:

π ~ D(μe) = Γ(μe K) / Γ(μe)^K * π^(μe−1) (12)

where D(μe) is the Dirichlet distribution with concentration parameter μe = 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(v, u, c, π, α, ω^2, ρ, ϵ^2 | x) ∝ p(α)p(ω^2)p(ρ)p(ϵ^2)p(π)p(c | π)p(u | α, ω^2, ρ, c)p(v | u, ϵ^2)p(x | 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 methods 89. 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 diagnostics 89 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 groups 22. 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 sheets 1. 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 SSP 1-2.6 (Shared Socioeconomic Pathway 1 paired with an approximate radiative forcing of 2.6 W/m2 by 2100 CE) and SSP 5-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 movement 1,46–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.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Reviewer #1 and #2 Evaluations:\nThe manuscript presents projections of 21st-century sea level fall along the coast of Greenland. The authors use available ice mass change estimates from the Last Glacial Maximum to 2100, derived from various ice models and reconstructions, to compute relative sea level change using the GIA model SEAKON for different Earth structures. They analyze separately the contributions of GIA-related processes associated with past deglaciation, the Little Ice Age, and future ice mass loss. In addition they use to non-GIA-related processes from the most recent IPCC WG1 report. The ensemble mean results indicate that sea level change in 2100 is primarily driven by future ice mass loss and non-GIA-related contributions, with projected maximum sea level fall of 1 ± 0.3 m for RCP 2.6 and 2.8 ± 0.8 m for RCP 8.5.\nGeneral Remarks:\nThe study is conducted thoroughly, and the manuscript is generally well written.\nHowever, the broader impact of the study is superficially addressed and remains somewhat unclear. While it is acknowledged that reliable sea level projections have significant implications for Greenlandic coastal communities, this aspect is only briefly mentioned and not elaborated upon in the manuscript.\nThe authors highlight the limitations of traditional approaches that consider only elastic deformation in projections. However, it is unclear whether they are referring to sea level projections based on GIA models or dynamical ice sheet models. Ice sheet models incorporate elastic and viscous deformation in different ways, often assuming a constant present-day rate. It is not clear in the manuscript how that compares to the results of this paper.\nThe results section assumed prior knowledge on linear, nonlinear and transient behavior of the Earth’s interior. Furthermore, different rheologies are used on different timescales, which is not explained in the manuscript. The introduction should include an explanation of these different types of deformation and its relevance on different timescales. Too much GIA expertise is presumed for the audience of Nat. Comm.\nThe spatial figures present ensemble means, but the standard deviation is relatively large due to variations among ice sheet models. Given this variability, displaying only the ensemble mean may not be the most informative approach. It would be beneficial to also present the mean per ice sheet reconstruction, as this likely determines the main pattern of deformation. In figures such as Figure 3, the relative error appears spatially consistent. It would be more informative to include a range of relative errors in the text and provide additional figures as suggested above.\nLine-by-Line Comments\nL12: Replace However by As a consequence\nL14: Predictions is used throughout the manuscript where projections are meant, please change\nL20: It is completely unclear what is meant with traditional approaches. Some applications have GIA as a constant offset. Does traditional approaches refer to GIA to ice sheet models or to IPCC style sea level projections\nL22: can these later two statements be more detailed. For a tenth of the change the sentence can still be made but the impact is negligible.\nL66-69: Provide the full range of upper mantle viscosity and lithospheric thickness values from references 5, 26, and 27 to give the reader a better overview of the uncertainty rather than highlighting a single result.\nL72: The reference to an AGU abstract (ref 29) appears redundant, given the availability of peer-reviewed sources supporting this statement.\nL71-74: Reference 33 states that lateral structure does not explain misfits to RSL data. Could you clarify how this compares to refs 5, 8, and 23? It is mentioned that the fit improves but it is not clear whether misfits remain after including lateral variations? If ref 33 also includes lateral variations, consider citing it together with refs 5, 8, and 23.\nL74-75: The phrase “observations beyond paleo RSL data” is somewhat vague. Could you specify what is meant here?\nL79: Define the time frame for “since the last deglaciation” and “contemporarily.” The period from the Last Glacial Maximum to the start of the Little Ice Age (~1000 CE) is defined in L114-115 but should be introduced earlier to improve readability. Could the authors clarify the exact time period they are referring to?\nL92: Do you mean during or since the LIA?\nL94: A non-GIA person will not understand why mantle viscosity would be smaller since the LIA, as it is suggested to be a material property. Please rephrase to make this clear for non-GIA specialists.\nL97: Is ref 39 the only reference discussing higher-order linear viscoelasticity and are refs 40 and 41 the only references to related to non-linear rheology? Please clarify the distinction.\nL102-103: The term “spatially extensive asthenosphere” is subjective. In Antarctica, an order of magnitude difference in mantle viscosity significantly affects the response, and similar variations can be expected between east and west Greenland. Please provide further clarification.\nL120: Between 1000 CE and 2100, VLM observations are corrected for coeval (and deglacial) ice mass change. However, coeval ice mass change is included in the GPS signal. Why is it corrected for separately?", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_3": "Review of \"Projections of 21st Century Sea-Level Fall Along Coastal Greenland\"\nThis manuscript aims at assessing the relative sea level along the coast of Greenland in the 21st century in response glacial\nisostatic adjustment and past and future ice-mass changes. My area of work is statistical methodology, and I am no expert in\nglaciology, so I will narrow my comments on the technical aspects of the methods.\nThe methods section is mostly a description of very basic approaches such as weighted sum of squares, decomposition of\nthe uplift in each station into a sum of different terms and linear regression, so there is not much to comment on, I assume the\nchoice of models was thorough and that the choice of these methods is based on well-established literature.\nThe interesting part is the Bayesian framework. This part is written unusually well for an applied work, so I commend the\nauthors for their effort. They propose a hierarchical Bayesian model for the GNSS observations, where the true VLM rates\nare a combination of 14,040 GIA model solutions.\nSome comments\n- No validation. The model makes many assumptions, but none are validated. For example, is x at least approximately\nGaussian? If one were to fit, say, a skew-Normal or some other members of the exponential family, how would the results\nchange? What about the exponential covariance function in (7), would it change if we had a Matérn covariance? Such\nassessment should have been done and their predictive skills measures with some cross-validation.\n- Lack of a continuous spatial process. While the model proposed is able to explain through the vector c which model better\nexplains the data, it does not assume any continuous spatial process. I am assuming the models have been interpolated on\nthe same grid as the data to avoid this issue. This, however, begs the question: to what extent is the interpolation impacting\nthe results? A much better solution would have been to provide a spatial map of a Multinomial distribution.\n- Lack of details on software and hardware. There is a wide variety of software able to perform Gibbs sampling and MCMC\n(e.g., STAN, JAGS, INLA, etc.). Which one has been used here and why? What computational resources have bene used?\nWhat was the computational time?\n- Lack of reproducibility. As of now, this work is not reproducible. No code has been provided, so I am not in the position of\ncommenting on it. Given the high standing of the journal, I would have expected, at the very least, a detailed commented\ncode on a GitHub repository for all methods used. This is in my opinion a considerable shortcoming for this submission.", "reviewer_4": "In this work, the authors present updated projections of 21st century sea-level change for Greenland, including updated VLM\nestimates generated by combining GIA models and VLM observations from permanent Global Navigation Satellite Systems\n(GNSS). My expertise is in sea level projections more broadly, and not in the specifics of GIA modeling/observations.\nTherefore, I have focused most of my comments below on the authors' generation of complete sea-level change projections,\nrather than the new VLM estimates.\nThe manuscript is generally well-written, and most figures are clear. However, while I think this manuscript has the potential\nto be a meaningful contribution to the sea level and GIA communities, I also have some concerns about the methods the\nauthors are currently using to generate complete sea-level change projections. If the authors can update their approach and\nthe explanation of their methods, I would look forward to revisiting the revised version of the study. My detailed comments\nare outlined below.\n~Andra J. Garner\n------\nMajor Comments\n--\nReading through the manuscript, I had quite a few questions about the methodology for the total RSL projections. The <10\nline explanation of these methods (Lines 523-531) generated more questions for me than it answered.\nLines 523-526 say, \"In addition to sea-level changes caused by GIA, other processes will contribute to the RSL experienced\nby coastal Greenland over this century. These processes include changes in global glaciers, sterodynamic changes, and\nchanges in terrestrial water storage, excluding glaciers and ice sheets. Note that the contribution from Antarctica is already\nincluded in the previous section.\"\nIf Antarctic ice loss is already included in the LIA to 2100 Ice Mass change, the first sentence of this section (and much of the\nresults) should be amended to reflect this. Because, in this case, this section is not \"Non-GIA related contribution to future\nsea-level change\", but \"Contributions to future sea-level change excluding GIA and Ice Sheet Mass Loss\". This fact was not\nclear to me until I reached line 526 of the manuscript.\nLines 526-528 state, \"Estimates for the contributions of each of these processes to RSL in Greenland to 2100 CE for RCP\n2.6 and RCP 8.5 were obtained from the IPCC AR6.\"\nThis is very confusing, since the IPCC AR6 did not make any explicit sea-level change projections for RCP scenarios, but\ninstead focused on SSP scenarios. Furthermore, to my knowledge, the IPCC AR6 did not provide RSL projections for\nGreenland in the full report--though the datasets made available for the IPCC AR6 projections (which should include\nadditional citations beyond AR6 Chapter 9) did. So, this description leaves me with a lot of questions about where the\ncontributions for each of these \"non-GIA\" processes actually came from.\nLines 528-529 say that, \"The mean of these estimates was made relative to 2017 CE and linearly interpolated at 5 year\nintervals between 2020 and 2100 CE.\" What does this mean for AR6 projections that are provided in terms of percentiles\n(median and likely range)? How did the authors construct a \"mean\" from these projections that are provided as a\ndistribution?\nIn addition, the IPCC AR6 projections are provided for a 1995-2014 baseline--so how were the authors' estimates made\nrelative to 2017? Furthermore, when the authors say that these estimates were linearly interpolated at 5 year intervals\nbetween 2020 and 2100, does that mean that the first interpolation was actually 3 years (going from 2017-2020), or should\nthe text say that the estimates were linearly interpolated at 5-year intervals from 2022-2097?\nLines 529-531 say that, \"For our total estimate of RSL around Greenland for the 21st century, we sum GIA-related\ncontributions associated with the deglaciation and the LIA to 2100 CE, and non-GIA related contributions.\" How does this\nwork with the uncertainty the authors provide? Are the authors also linearly summing uncertainty across components? This\nneeds explanation, since, again, the AR6 projections are provided in terms of percentiles, and unless a median estimate is\nbeing used, the uncertainty across percentiles cannot just be added together.\nIt seems to me that the correct approach for the analysis the authors present would be to use the\n\"IPCC AR6 Relative Sea Level Projections without Background Component\" dataset, which is available here: https://zenodo.org/records/5967269.\nThis dataset was specifically designed to allow scientists to use their own estimates of VLM to generate new RSL\nprojections in a manner consistent with IPCC methodology. Using this dataset, my suggestion would be to simply use the\nAR6 ice sheet projections (also from ISMIP, which the authors use for their LIA-2100 ice-mass change, and which are\nalready incorporated into the workflows as part of this dataset). The authors could still include the VLM component\ncalculated from their LIA-2100 ice-mass change for their GIA/VLM \"background\" component of RSL. But, I think from a\nstandpoint of correctly combining components within the AR6 framework, it's best to use the ice sheet estimates that\nhave already been included within the framework in this dataset.\nNote that when using this dataset, there are a number of citations/acknowledgements that should be included, as noted on\nthe Zenodo page:\n\"In order to document the impact of these sea-level rise projections, users of the projections are obligated to cite chapter 9 of\nWorking Group 1 contribution to the the IPCC Sixth Assessment Report, the Framework for Assessment of Changes To Sea-\nlevel (FACTS) model description paper, and the version of the data set used:\n• Fox-Kemper, B., H.T. Hewitt, C. Xiao, G. Aðalgeirsdóttir, S.S. Drijfhout, T.L. Edwards, N.R. Golledge, M. Hemer, R.E. Kopp,\nG. Krinner, A. Mix, D. Notz, S. Nowicki, I.S. Nurhati, L. Ruiz, J.-B. Sallée, A.B.A. Slangen, and Y. Yu, 2021: Ocean,\nCryosphere and Sea Level Change. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I\nto the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [Masson-Delmotte, V., P. Zhai, A. Pirani,\nS.L. Connors, C. Péan, S. Berger, N. Caud, Y. Chen, L. Goldfarb, M.I. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J.B.R.\nMatthews, T.K. Maycock, T. Waterfield, O. Yelekçi, R. Yu, and B. Zhou (eds.)]. Cambridge University Press, Cambridge,\nUnited Kingdom and New York, NY, USA, pp. 1211–1362, doi:10.1017/9781009157896.011.\n• Kopp, R. E., Garner, G. G., Hermans, T. H. J., Jha, S., Kumar, P., Reedy, A., Slangen, A. B. A., Turilli, M., Edwards, T. L.,\nGregory, J. M., Koubbe, G., Levermann, A., Merzky, A., Nowicki, S., Palmer, M. D., & Smith, C. (2023). The Framework for\nAssessing Changes To Sea-Level (FACTS) v1.0: A platform for characterizing parametric and structural uncertainty in future\nglobal, relative, and extreme sea-level change. Geoscientific Model Development, 16, 7461–7489.\nhttps://doi.org/10.5194/gmd-16-7461-2023\n• Garner, G. G., T. Hermans, R. E. Kopp, A. B. A. Slangen, T. L. Edwards, A. Levermann, S. Nowikci, M. D. Palmer, C. Smith,\nB. Fox-Kemper, H. T. Hewitt, C. Xiao, G. Aðalgeirsdóttir, S. S. Drijfhout, T. L. Edwards, N. R. Golledge, M. Hemer, G. Krinner,\nA. Mix, D. Notz, S. Nowicki, I. S. Nurhati, L. Ruiz, J-B. Sallée, Y. Yu, L. Hua, T. Palmer, B. Pearson, 2021. IPCC AR6 Sea\nLevel Projections. Version 20210809. Dataset accessed [YYYY-MM-DD] at https://doi.org/10.5281/zenodo.5914709.\nPlease also include in the acknowledgements of works citing these projections:\nWe thank the projection authors for developing and making the sea-level rise projections available, multiple funding\nagencies for supporting the development of the projections, and the NASA Sea-Level Change Team for developing and\nhosting the IPCC AR6 Sea-Level Projection Tool.\"\nTaking this approach for the total sea level projections (and making it clear that this was the approach) will greatly\nstrengthen the authors' discussion and comparison of their results to the AR6 regional results for Greenland. This approach\nwill ensure that all components (and scenarios) considered are the same as in the AR6 regional projections, allowing the\nauthors to more vividly emphasize the importance of their work on Greenland VLM when it comes to accurately assessing\nfuture SLR change for the region.\nFinally, the abstract says that the authors predict RSL fall of 2.8 +/- 0.8 m for RCP 8.5. Similarly, Lines 292-293 say, \"The fall\nin RSL by up to almost 3 m that we project to occur by the end of this century along much of the coast . . .\" However, looking\nat Figures 5-7, and considering the accompanying text, it's unclear how the authors conclude that their findings suggest\n\"almost 3 m\" of RSL fall \"along much of the coast\". It looks like up to ~2 m in many places, not 3. Please provide clarification\nof how the 3 m value is determined--right now I don't see data/results to support this conclusion.\n------\nMinor Comments\n--\nLines 259-261: \"For all locations considered here, it is noteworthy that the projections for the two future emissions scenarios\n(RCP2.6 vs. 8.5) only start to diverge around 2050 CE.\" This may be noteworthy, but it's probably also worth mentioning that\nthis is consistent with what many other studies have already shown (e.g., Garner et al., 2018, Horton et al., 2018, IPCC AR6\nChapter 9):\n• Garner, A. J., J. L. Weiss, R. M. Horton, A. Parris, R. E. Kopp, J. T. Overpeck, B. P. Horton (2018), Evolution of 21st Century\nSea-level Rise Projections. Earth’s Future. doi: 10.1029/2018EF000991.\n• Horton, B. P., R. E. Kopp, A. J. Garner, C. C. Hay, N. S. Khan, K. Roy, and T. A. Shaw (2018), Mapping sea-level change in\ntime, space, and probability. Annual Review of Environment and Resources. 43:1. doi:10.1146/annurev-environ-102017-\n025826.\n• Fox-Kemper, B., H.T. Hewitt, C. Xiao, G. Aðalgeirsdóttir, S.S. Drijfhout, T.L. Edwards, N.R. Golledge, M. Hemer, R.E. Kopp,\nG. Kopp,\nG. Krinner, A. Mix, D. Notz, S. Nowicki, I.S. Nurhati, L. Ruiz, J.-B. Sallée, A.B.A. Slangen, and Y. Yu, 2021: Ocean,\nCryosphere and Sea Level Change. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I\nto the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [Masson-Delmotte, V., P. Zhai, A. Pirani,\nS.L. Connors, C. Péan, S. Berger, N. Caud, Y. Chen, L. Goldfarb, M.I. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J.B.R.\nMatthews, T.K. Maycock, T. Waterfield, O. Yelekçi, R. Yu, and B. Zhou (eds.)]. Cambridge University Press, Cambridge,\nUnited Kingdom and New York, NY, USA, pp. 1211–1362, doi: 10.1017/9781009157896.011.\nPlease ensure that the baseline for total RSL projections is stated clearly at least once in the results section for Future\nProjections of RSL, and in the captions of all figures that show total RSL projections, to provide clarity to readers.\n------\nComments on Figures\n--\nThe maps generally look very good, however the colorbar on some of them becomes quite dark (e.g., Fig. 3b, Fig. 4b,d,f, and\nh), making portions of the map somewhat difficult to see. Would it be possible to adjust the colorbar to keep it lighter?\nThe station numbers on the inset map of Greenland in Figure 2 are not legible. I think this could be solved by making the\nmap of Greenland larger--perhaps make it a third panel to the right of panels a) and b)?\nThe white stars/black points on Figure 5d don't work well. The black points are particularly difficult to see. I'd suggest putting\nthe black points on a different panel than the white stars . . . And perhaps make them a different color and/or shape, so that\nthey don't blend into the map outlines.\nIn Figure 8, it seems odd to discuss panel c) before a) and b) in the caption. Perhaps consider reordering the figures, or\nrevisiting the caption."}, "author_response": "We thank the reviewers for their detailed and positive feedback. We have now revised the\nmanuscript to address the comments that have been raised. We describe these changes below by\ninterspersing the original feedback in black with our response in blue. We note that all our line\nnumbers refer to the revised manuscript. We also include a version of the manuscript that\nindicates all changes to easily track changes from the original manuscript.\nReviewer #1 and #2 Evaluations:\nThe manuscript presents projections of 21st-century sea level fall along the coast of Greenland.\nThe authors use available ice mass change estimates from the Last Glacial Maximum to 2100,\nderived from various ice models and reconstructions, to compute relative sea level change using\nthe GIA model SEAKON for different Earth structures. They analyze separately the\ncontributions of GIA-related processes associated with past deglaciation, the Little Ice Age,\nand future ice mass loss. In addition they use to non-GIA-related processes from the most recent\nIPCC WG1 report. The ensemble mean results indicate that sea level change in 2100 is primarily\ndriven by future ice mass loss and non-GIA-related contributions, with projected maximum sea\nlevel fall of 1 ± 0.3 m for RCP 2.6 and 2.8 ± 0.8 m for RCP 8.5.\nGeneral Remarks:\nThe study is conducted thoroughly, and the manuscript is generally well written.\nHowever, the broader impact of the study is superficially addressed and remains somewhat\nunclear. While it is acknowledged that reliable sea level projections have significant implications\nfor Greenlandic coastal communities, this aspect is only briefly mentioned and not elaborated\nupon in the manuscript.\nThe authors highlight the limitations of traditional approaches that consider only elastic\ndeformation in projections. However, it is unclear whether they are referring to sea level\nprojections based on GIA models or dynamical ice sheet models. Ice sheet models incorporate\nelastic and viscous deformation in different ways, often assuming a constant present-day rate.\nIt is not clear in the manuscript how that compares to the results of this paper.\nThe results section assumed prior knowledge on linear, nonlinear and transient behavior of the\nEarth’s interior. Furthermore, different rheologies are used on different timescales, which is not\nexplained in the manuscript. The introduction should include an explanation of these different\ntypes of deformation and its relevance on different timescales. Too much GIA expertise is\npresumed for the audience of Nat. Comm.\nThe spatial figures present ensemble means, but the standard deviation is relatively large due to\nvariations among ice sheet models. Given this variability, displaying only the ensemble mean\nmay not be the most informative approach. It would be beneficial to also present the mean per\nice sheet reconstruction, as this likely determines the main pattern of deformation. In figures\nsuch as Figure 3, the relative error appears spatially consistent. It would be more informative to\ninclude a range of relative errors in the text and provide additional figures as suggested above.\nWe thank the reviewer for their positive assessment and for their constructive feedback\nand summary of the main concerns here. Most of the points above are also reflected in\nthe Line-by-Line Comments below and we therefore address them as they appear in the\ncomments below. The only comment that is solely raised here relates to the variation in\nsea-level change that is driven by differences among ice-sheet models used in our\nprojections (last paragraph). We have included an additional figure in the supplement that\nshows spatial projections of sea-level change at 2100 CE for each ice sheet modelling\ngroup used in this study (Supplementary Fig. 4).\nLine-by-Line Comments\nL12: Replace However by As a consequence\nWe agree that both “However” and “As a consequence” could be used here. However,\nwe do believe that our use of “However” highlights the unexpectedness of projected\nsea-level fall in Greenland (considering the ice sheet is causing global sea level to rise)\nand therefore better fits the narrative that we are attempting to tell through this paper. We\ntherefore decided to keep the original wording.\nL14: Predictions is used throughout the manuscript where projections are meant, please change\nWe thank the reviewer for noticing this and have changed the wording from predictions\nto projections where relevant.\nL20: It is completely unclear what is meant with traditional approaches. Some applications have\nGIA as a constant offset. Does traditional approaches refer to GIA to ice sheet models or to IPCC\nstyle sea level projections\nWe are referring to approaches in which sea-level projections driven by projected\nice-mass change use an elastic Earth structure to calculate deformation (as is used, for\nexample, in refs. 1, 62, 63). This statement has been removed from the abstract to adhere\nto the word limit, but we have added these references to our discussion on this point\n(lines 305-316).\nL22: can these later two statements be more detailed. For a tenth of the change the sentence can\nstill be made but the impact is negligible.\nWe agree that this sentence is brief and have added some details that allow us to still\nadhere to the word limit of the abstract (lines 14-16):\n“A fall in relative sea level has implications for Greenlandic communities whose\neconomy, near-shore infrastructure, and food security are vulnerable to coastal\nchanges, and for the dynamics of marine-terminating glaciers”\nFor more details on these two statements, we refer to our results and discussion section,\nwhere we have added text to better describe the implications of our study (lines 317-346).\nL66-69: Provide the full range of upper mantle viscosity and lithospheric thickness values from\nreferences 5, 26, and 27 to give the reader a better overview of the uncertainty rather than\nhighlighting a single result.\nWe agree that including the preferred values of these parameters from relevant studies\nwould enhance the reader’s understanding of their uncertainty. We have added to the text\n(lines 65-72):\n“Lecavalier et al. [5] find that GIA models forced by a GrIS deglaciation history\n(termed Huy3) derived from a dynamical ice sheet model best fit RSL data in\nGreenland when using a 120 km thick elastic lithosphere and upper and lower\nmantle viscosities of 5 × 1020 Pa s and 2 × 1021 Pa s, respectively. However, these\nvalues are not consistent across all studies that tune their model parameters to fit\nGreenland RSL data [23, 24]. Tarasov & Peltier [23] find that their GrIS history\nperforms well when coupled with a lithosphere of 90 km. Fleming & Lambeck\n[24] also shows preference for a thinner lithosphere (80 km), in addition to higher\nvalues for the viscosities of the upper mantle (4 × 1021 Pa s) and lower mantle (10\n× 1021 Pa s).\nL72: The reference to an AGU abstract (ref 29) appears redundant, given the availability of\npeer-reviewed sources supporting this statement.\nWe appreciate the reviewer pointing this out. This was an oversight on our part and the\nreference has been removed from the revised manuscript.\nL71-74: Reference 33 states that lateral structure does not explain misfits to RSL data. Could you\nclarify how this compares to refs 5, 8, and 23? It is mentioned that the fit improves but it is not\nclear whether misfits remain after including lateral variations? If ref 33 also includes lateral\nvariations, consider citing it together with refs 5, 8, and 23.\nWe realized that this portion of the text was unclear and have added further explanation.\nReferences 5 and 8 found that using different 1D Earth parameters for the East and West\ncoasts of Greenland improved their ability to fit RSL data. Ref 33 (now 32) followed on\nthese studies by running a GIA model with a 3D varying Earth structure and found that it\nwas unable to explain all of the misfits to the RSL data that were pointed out in ref 5. We\nremoved former ref 23 from this portion of the text as it noted how a laterally varying\nEarth structure could improve model fit to GPS data, and in this section we are talking\nabout RSL data. The revised text now reads (lines 72-78):\n“A difference in the optimal 1D Earth parameters inferred for the East and West\ncoasts of Greenland suggest that improvements in the misfit between modeled\npredictions and RSL data could be achieved by considering a more complex,\nlaterally-varying (i.e., 3D) Earth structure [5, 8]. The influence of a 3D Earth\nstructure on GIA predictions is shown to be significant in numerous studies\nfocusing on various regions around the world [e.g. 28–32]. In Greenland,\nHolocene RSL predicted by 3D GIA models can differ by tens of meters when\ncompared to the 1D case, however it does not necessarily explain all\npreviously-indicated misfits to the RSL data [32].”\nL74-75: The phrase “observations beyond paleo RSL data” is somewhat vague. Could you\nspecify what is meant here?\nWe have changed the sentence to more explicitly describe what we intended (lines\n78-79). We use this sentence as a segway into the following paragraph, where we\ndescribe contemporary observations of GIA in more detail:\n“Incorporating observations of present-day deformation can be useful to further\nconstrain values for the properties of Earth’s interior.”\nL79: Define the time frame for “since the last deglaciation” and “contemporarily.” The period\nfrom the Last Glacial Maximum to the start of the Little Ice Age (~1000 CE) is defined in\nL114-115 but should be introduced earlier to improve readability. Could the authors clarify the\nexact time period they are referring to?\nWe agree that being specific about what is meant by these two time frames would be\nhelpful for the reader and have reworded it for clarity (lines 81-84):\n“Present-day VLM recorded at the 57 permanent GNSS stations that span coastal\nGreenland, called GNET, is caused by deformation driven by ice changes that\nhave occurred in the past (most importantly over the last deglaciation, i.e. last 26\nkyr) and that are occurring today (i.e. during the recorded period)[33].”\nL92: Do you mean during or since the LIA?\nWe mean during and since the LIA but agree that this a little confusing. We have\nreworded this to “since the beginning of the LIA”(line 97).\nL94: A non-GIA person will not understand why mantle viscosity would be smaller since the\nLIA, as it is suggested to be a material property. Please rephrase to make this clear for non-GIA\nspecialists.\nWe agree that our original explanation of transient rheology was lacking and could be\nfurther elaborated. We have added to the original paragraph to make the explanation of\ntransient rheology understandable for the broader audience (lines 100-109):\n“They suggest that this lowering in mantle viscosity for ice-sheet changes since\nthe beginning of the LIA hints at the presence of transient deformation beneath\nGreenland. Transient deformation is an extension of the classical Maxwell\nrheology. While the latter consists only of an instantaneous elastic response and a\nsteady-state viscous creep, transient deformation, as is, for example, captured by a\nBurgers rheology, allows for additional dissipative mechanisms which act to\nlower the apparent viscosity over short periods. This deformation style is\nsupported by high-temperature high-pressure laboratory experiments of rock\ndeformation [36, 37] and has been explored in other GIA studies [e.g. 38-42]. In\nthe case of loading events that occur over centennial timescales (like the LIA), the\nviscosity of the mantle can effectively be lower than that inferred for\nlonger-timescale events (like the LGM or mantle convection), allowing the mantle\nto deform more quickly and with a greater magnitude.”\nL97: Is ref 39 the only reference discussing higher-order linear viscoelasticity and are refs 40 and\n41 the only references to related to non-linear rheology? Please clarify the distinction.\nWe decided to rephrase this sentence in response to the comment above and now simply\nrefer to past studies that have explored transient rheology in GIA modeling (line 106).\nL102-103: The term “spatially extensive asthenosphere” is subjective. In Antarctica, an order of\nmagnitude difference in mantle viscosity significantly affects the response, and similar variations\ncan be expected between east and west Greenland. Please provide further clarification.\nWe agree with the reviewer that our wording made the point we were trying to make\nambiguous. In Pan et al. (2024), they insert a low-viscosity layer beneath a 75 km thick\nlithosphere and on top of an upper mantle of viscosity 5x1020 Pa s. They show that\nincluding this low-viscosity channel allows them to fit both present-day VLM and paleo\nRSL data in Greenland. This study shows that transient rheology is not required to\nreconcile these datasets. However, seismic tomography studies (cited in line 115) have\nshown high velocity anomalies under portions of the lithosphere in Greenland,\nhighlighting old cratonic regions, and contradicting the existence of a homogenous\nsub-lithosphere low-viscosity layer at one depth. The authors have changed our\nterminology to make our point more clear (lines 112-115):\n“However, Pan et al. [44] caution that a sub-lithospheric low-viscosity layer with\na regular Maxwell rheology would also be able to explain both RSL data and\npresent-day VLM if the LIA ice history is included in the GIA model. Further\nwork is needed to understand the possible extent of a sub-lithospheric\nlow-viscosity layer around and below the cratonic keel of Greenland [45, 46].\nL120: Between 1000 CE and 2100, VLM observations are corrected for coeval (and deglacial)\nice mass change. However, coeval ice mass change is included in the GPS signal. Why is it\ncorrected for separately?\nHere we consider the GPS signal to be a result of viscoelastic deformation in response to\nice mass change over three distinct time periods: the deglacial (30-1ka), the LIA (1ka to\nthe beginning of the GPS record), and the duration of the GPS record. We predict the\nsecond of these contributions and want to compare it to the observation. We therefore\nneed to remove the first and third contribution (deglacial and coeval) from the\nobservation. In response to the reviewers comment on L154-164 we have expanded on\nthis here.\nL131-136: Ref 44 is not mentioned in the introduction but should be placed in context in the\nintroduction to improve readability. Refs 8, 33, and 44 are not discussed further, but ref 5\nis—why? Additionally, why are refs 15 and 45 not mentioned alongside refs 5, 8, 33, and 44 in\nL132?\nWe agree that Ref 44 (now ref 33) is an important study that should be included in the\nintroduction and we have added it to line 84. We realized we cited the wrong study from\nSimpson et al. (ref 8 in original manuscript) by accident. We have now cited the correct\nstudy in the revised manuscript (ref 26). Reference 5 is included in the subsequent\ninterpretation of our present-day VLM caused by the deglaciation to support the\nstatements made about the ice history, as this reference provides a detailed chronology of\nice extent for Greenland from the LGM. Ref 15 isn’t mentioned in Line 155 because it\nuses deglacial uplift rates produced by ref 32 (former reference 33), so including it would\nbe redundant. Ref 45 (now reference 53) is not included in Line 155 as it does not\nproduce predictions of Greenland-wide uplift rates for us to compare our results to. We\ngeneralized text in line 157 by changing “6-ka” to “mid-Holocene” so that our statement\nwould be consistent with all the deglacial ice histories that we cite.\nL143: Add a brief explanation of what “normalized weights” mean at the start of the paragraph.\nWe have changed the sentence to clarify this (lines 166-168):\n“The weights for the deglacial GIA models (normalized so that they sum to one)\nfall within the range of 0.0081 to 0.0227, and do not deviate substantially from\nthe value that would be assigned if all models were weighted equally (1/65 =\n0.0154).”\nL143-152: Since ice history is a key factor, please consider adding a supplementary figure\nshowing average ice mass change from ICE6G, Huy3, and other reconstructions/models. Include\na brief description of these reconstructions.\nWe thank the reviewer for this suggestion and do agree that the ice history is a key factor\nwhen modelling GIA. However, we consider a more detailed description of the deglacial\nice histories to be beyond the scope of this (already long) paper since it is easily findable\nin the cited literature.\nL154-164: The method for removing elastic and GIA signals from observations is not entirely\nclear. Please include a brief explanation in the introduction of this section. Additionally, a short\nexplanation of the Bayesian approach in L120 would help in understanding Figure 2.\nWe have added details to the introduction of this section in order to make these methods\nmore clear (lines 130-140):\n“For (2), we use a new suite of 14,040 1D Maxwell GIA simulations. We first\npredict present-day VLM in Greenland driven by the LIA and compare these to\nobservations of uplift measured at 55 GNET stations. In order to make a direct\ncomparison, we correct the observed VLM for uplift associated with deglacial and\ncoeval mass changes, the former of which is calculated in step (1) and the latter is\ntaken from Berg et al. [19]. We further account for potential differences between\nthe reference frame used in the 1D GIA model and that used for the VLM\nobservations. We perform a Bayesian inversion to identify the distribution of most\nlikely ice and Earth parameters varied in the 1D GIA model and create a posterior\nprediction of present-day VLM driven by the LIA in Greenland. We then use\nthese optimal parameters to create a suite of 1D Maxwell GIA runs forced by\nLIA-to-2100 CE ice-mass change, using an ensemble of projections of Greenland\nand Antarctic ice change for both Representative Concentration Pathway (RCP)\n2.6 and RCP 8.5 [47, 48].”\nAd Figure 2: very difficult to see a difference between prior and posterior. Can you plot a panel a\nthe prior and as panel b the posterior-prior or maybe in one plot observed-prior and\nobserved-posterior. Alternatively, We suggest to show all three components\n(observations-ELA-GIA, Prior LIA and Posterior LIA) in one figure using overlapping boxplots\nin three different colors instead of shaded areas.\nWe agree that it is hard to see the difference between the prior and posterior solutions in\nthis figure and thank the reviewer for their suggestions on ways in which this figure could\nbe improved. We have replotted this figure so that the observations and predictions for\nboth the posterior and prior are all in one figure (Figure 2). We have also added to the text\nto describe why the posterior and prior means look similar (lines 184-186) :\n“There is only a slight difference in mean between our prior and posterior\nsolutions, as the model parameters used in this study (Table 1) were sampled from\nthe posterior distribution of parameters from Adhikari et al. [15].”\nAd Figure 3: How large is the relative standard deviation? 15%? seems an insightful quantitative\nvalue. It also looks like there is barely a spatial pattern in the standard deviation it is more or less\na relative error. So what is the added value of panel b, maybe you can leave it out and just\nmention the relative error.\nThe average relative error taken across GNET stations is 27%. We have chosen to leave\nthe original figure as is, but have included this number in the text (line 191-192):\n“The standard deviation of the posterior solution across GNET stations is 20-30%\nof the posterior mean.”\nL171: It is not clear whether the weighted average LIA vertical land motion of all simulations is\nsimilar to the average of Adhikari. Could you clarify this?\nWe now clarify in line 196 that we refer to the posterior (weighted) LIA vertical land\nmotion.\nL180-184: The finding of low upper mantle viscosity does not imply that the entire upper mantle\nhas this viscosity, as a 1D model is used. Rather, it suggests that the regions where the mantle is\nsensitive exhibit this viscosity, while other areas (insensitive to the ice load) may have different\nviscosities. Could you clarify how this leads to an indication of transient viscosity approaching\nsteady state? Why would the mantle viscosity be smaller since the LIA?\nThe reviewer is correct that our findings do not indicate that the entire upper mantle has\nthis viscosity and we changed the sentence to reflect that we are only referring to the\nmantle below Greenland. Inferring a lower viscosity for ice changes since the LIA\ncompared to the estimates of viscosity inferred in Greenland over the deglaciation (e.g.\n5x1020 Pa s, ref 5) supports the idea of transient rheology, as it suggests that the\nfrequency of loading (i.e. how quickly ice grows/melts) can impact the apparent viscosity of the\nmantle’s response to that loading. We say that the viscosity would approach that of the\nsteady-state viscosity on deglacial timescales because a higher viscosity is expected on\nlonger timescales within the transient regime of deformation and a higher viscosity for\nthe deglaciation has been inferred from previous studies (see lines 65-72). However,\nsince our study does not explicitly infer deglacial-timescale viscosity, we removed this\npart of the sentence. We also added an explanation to the introduction (lines 100-109) that\nexplains transient rheology better (and why the mantle viscosity would be smaller over\nthe LIA). Lastly we changed the wording in the sentence that is being referred to here\nthat our viscosity inference may be indicative of transient deformation (rather than is)\ngiven the results by Pan et al., (2023).\nL194 Do you intend to say effect the results only marginally?\nThe authors were intending to say that both the data and the inference methods that we\nuse lead to the minor differences in results between this study and Adhikari et al. (2021)\n(ref 15). We have changed “affect” to “leads to” in line 220 to make our meaning clearer.\nL198 Repeat that this is the forebulge effect of North America, if that is what explains it.\nWe have revised this sentence for clarify (lines 225-228):\n“The deglacial contribution is the smallest component but causes a slight rise (<\n0.21 m ± 0.05 m) in RSL in west and southeast Greenland due mainly to\nLaurentide forebulge collapse and a fall in RSL (∼ −0.1 m ± 0.09 m) in the\nnortheast (Fig. 4a) caused by ongoing uplift in response to ice loss since the\nLGM.”\nL204 Unclear from where you took the Antarctic values.\nWe now summarize in the beginning of the results and discussion section where the\nAntarctic contribution comes from (lines 137-140). It is also further described in the\nmethods section (lines 537-541).\nL207: Could you specify which gravitational effects are included—self-gravitation of ice on\nrelative sea level, self-gravitation of the Earth on deformation, or both?\nIt’s both. We now describe this explicitely lines 235-237):\n“Gravitational effects include both the ice and the solid Earth’s pull on the ocean’s\nsurface and a uniform gravitational shift caused by the input of meltwater into the\nocean.”\nL208: You refer often to Table S2, would it be possible to absorb it in the main text?\nTable S2 has been absorbed in the main text as Table 1.\nL211: Please include the spatial distribution of gravitational effects on sea level in the\nsupplementary materials. The range of ~20-30%—what does this depend on? Was only one\nsimulation conducted? How much does this effect amount to in meters? Please also show the\nspatial distribution in the supplementary material.\nThe range depends on location and RCP scenario. We only used one set of LIA ice\nhistory and Earth structure parameters in our GIA model, but ran the model to 2100 for\nRCP 2.6 and RCP 8.5 using the ensemble of projections of ice change from the 8 ice\nsheet modeling groups used in this study. We have clarified this in lines 237-243:\n“We calculate this by running the GIA code from the LIA to 2100 CE using the\nensemble of ice sheet projections for RCP 2.6 and RCP 8.5 with one parameter\ncombination (Table 1) that performs well in fitting observations, and then\nseparating the two components. We find that RSL at 2100 CE caused by ice-mass\nchanges from the LIA to 2100 CE is dominantly driven by solid Earth\ndeformation, with sea-level changes caused by gravitational effects constituting ∼\n0-20% of the total for RCP 2.6 and ∼ 10-30% of the total for RCP 8.5, depending\non location around the coast (Supplementary Figure 3).”\nWe also included a spatial plot of the percentage of RSL change at 2100 driven by\nLIA-to-2100 CE ice-sheet mass change that is made up from gravitational effects in the\nsupplementary material (Supplementary Figure 3). This effect amounts to a maximum\ncontribution of approximately -0.2 m for RCP 2.6 and approximately -0.7 m for RCP 8.5,\nwhich is represented by the percentage range described in the text.\nL212-218: Are these results derived from literature rather than modelled? If so, please add\nreferences.\nWe thank the reviewers for pointing this out. The relevant references have been added to\nline 245.\nL240: So separating does not add much. I think you can describe one and argue that qualitatively\nthe other locations are yielding comparable results.\nWe thank the reviewer for this suggestion. We have described projected sea-level change\nfor Tasilaaq and Nuuk and have moved the analyses for the other towns/settlements to the\nsupplement.\nFigure 4: unclear what the white areas in panel i to l are (same in Figure 5)\nProjections for the thermal expansion, ocean dynamics, global glacier change, and land\nwater storage contributions to RSL do not exist within these areas, so we have excluded\nthe regions from our projections at 2100 CE. An explanation of this has been added to the\nfigure captions for figures 4 and 5.\nL260: I am not sure whether that is noteworthy as this holds for most points along the world\ncoastline.\nThis is a great point that was also brought to our attention by another reviewer. We have\ndecided to keep this sentence in the text, but have added references for other studies that\nhave made a similar finding (lines 279-281):\n“For all locations considered here, it is noteworthy that the projections for the two\nfuture emissions scenarios (RCP 2.6 vs. 8.5) only start to diverge around 2050\nCE, which is consistent with prior work on sea-level projections globally [e.g. 1,\n59, 60] (Supplementary Figure 5).”\nAd Figure 6: try to homogenize the lay-out of the figure. Unclear what the figure adds. I believe\nit can be left out.\nWe have moved this figure to the supplement and have attempted to homogenize it better.\nAd Figure 7: What is the added value of a second location I think only Nuuk is sufficient?\nExplain the red curve in the caption as well.\nThe break-down of the different factors that contribute to local sea-level projections will\ndiffer from place to place. Especially the contribution from ice change since the\nbeginning of the LIA can vary significantly. We believe that showing the break-down for\ntwo locations highlights this.\nL294: Not sure whether this is true. It assumes the marine terminating glaciers dynamics are\ncontrolled by water depth. I think that is questionable, at least make clear this is what you\nassume and provide reference for this assumption.\nThough not a primary control of marine-terminating glacier dynamics, we do assume that\nthe impact of GIA on bed topography and water depth at the grounding line modulates\nglacier behavior, as illustrated in several papers [refs 65, 66]. We have changed the text to\nexplicitly state how a RSL fall of the projected magnitude could impact the glacier (by\nchanging the position of the grounding line) and have added additional text and\nreferences to support this (line 317-326):\n“The fall in RSL by up to 2.8 m for RCP 8.5 that we project to occur along\nGreenland’s coasts by the end of this century could act to slow down the pace of\ngrounding-line retreat for Greenland’s marine-terminating glaciers, as the position\nof the grounding line depends on sea level at the grounding line [64-66]. Ice\ndischarge from marine-terminating glaciers was responsible for approximately\ntwo-thirds of Greenland’s total contribution to global sea level over the period\n1972–2018 CE [67, 68]. The response of the solid Earth and gravity field to this\nice-mass loss can raise the bed topography and decrease RSL at the grounding\nline, which can lead to its stabilization and potential re-advance [65, 66]. This\neffect may be particularly impactful for large marine-terminating glaciers that rest\non bedrock that deepens inland (like Kangerlussuaq and Sermeq Kujalleq [69])\nand are vulnerable to marine ice sheet instability feedbacks, where the flux of ice\nacross the grounding line increases as it retreats [64, 70-74].”\nL305: Adhikari et al. studied Antarctica and the response of the Antarctic Ice Sheet and its ice\nshelves cannot be directly compared to that of the Greenland Ice Sheet and its outlet glaciers due\nto fundamental differences in their physical processes and environmental conditions. Please\ninclude this nuance in the manuscript.\nWe thank the reviewer for this pointing this out and agree. We have included this nuance in\nour manuscript and use it to motivate the need for a Greenland specific coupled GIA-ice\nsheet modeling study (lines 328-333):\n“This amount of uplift is of similar magnitude to that projected to occur along\nGreenland’s coasts over a similar time frame, however a direct comparison\nbetween the dynamics of Greenland and Antarctica's marine-terminating glaciers\ncannot be made. This limitation, along with the potential for the GIA-induced\nstabilization of these marine-terminating glaciers to reduce Greenland’s\ncontribution to global sea-level rise, motivates the need for a detailed analysis that\ncouples GIA and ice-sheet models on centennial timescales in Greenland.”\nL307 can you explain this in more detail. Why has a ~ meter an impact? in view of the hazards\nand threats communities are exposed to?\nThe authors agree with the reviewer that the discussion of the impact of projected\nsea-level fall on Greenlandic communities can be expanded. We have added to the\nparagraph to reflect this (lines 335-344):\n“RSL fall on the order of multiple meters has the potential to increase\nnavigational hazards and limit the accessibility of harbors to deep draft boats by\nreducing the water depth beneath their keels. Marine channels that are already\nshallow today are vulnerable to becoming exposed at low tide, disconnecting\nnavigational routes and waterflow, which may impact marine habitats. Port\ninfrastructure, such as docks and mooring sites, could become less useful as the\nfalling sea surface would make access to them more difficult or even impossible.\nA reliance on global sea-level projections has led to the perspective among many\nresidents that sea level will rise in Greenland due to climate change [75], which\nwould create a very different set of hazards than those presented above. The\nregional RSL projections provided in this study can be used by communities to\nassess their site-specific risks to sea-level fall and prepare adaptation strategies\naccordingly.”\nL311: In the context of publishing the manuscript in Nature Communications, it would be\nvaluable to elaborate on the interdisciplinary aspects of the study. Could you clarify what\ninterdisciplinary studies are meant?\nWe now specify what we mean by interdisciplinary studies in this case (lines 344-346):\n“Subsequent interdisciplinary studies that consider these projections with regard\nto the national economy, coastal structure, and marine ecology are needed to\nbetter understand and communicate the impact that RSL fall could have on these\ncommunities.”\nAd Figure 8: If understood correctly I think the text discusses the relative contribution for the\nviscous part why not show that or show both?\nWe have changed Figure 8 panel d to show the viscous percentage of RSL at 2100 CE.\nL312: clarify that you prescribe ice sheet history and what that implies. a. Not a consistent\nglaciological model, and b not a coupled solution with a glaciological model.\nWe added 2 sentences (lines 367-370) to clarify this:\n“While the Greenland and North American portions of Huy3 were created using\nglaciological models that account for ice dynamics, ICE6G, ANU, and the other\nglobal ice masses in Huy3 were not. Further, none of the reconstructions used\nhere were created using a coupled GIA-ice sheet model.”\nL334: Confusing phrasing. it is the same model but 65 different members with different physical\ncomponent or parameter settings.\nWe agree that this phrasing is confusing. We have revised it to (line 373-374):\n“We use output from the suite of 65 GIA model runs described above to calculate\nthe contribution of the deglaciation to present-day VLM and future sea-level\nchange.”\nL335: What is meant by the latter I presume the present-day VLM; unclear what that precisely is.\nThe present-day rate or the VLM caused by the future melt?\nWe have revised this (line 374) to explicitly say ‘Present-day VLM’ rather than ‘The\nlatter’.\nL349: observed RSL at location y add the last bit.\nThe authors thank the reviewer for catching this small detail. We have decided to remove\nthe superscript y, as it was added in to suggest we were calculating the residual between\nthe observed elevation of RSL and the predicted. However, we have realized that since\nwe aren’t calculating any other residual, it is unnecessary to specify this within the\nequation, as it is clarified in the text.\nLine 354: can you explain in words what the second line in equation 1 secures. I guess it has to\ndo with the limiting aspect but unclear how that works.\nThe reviewer is correct that the second line in the equation treats the limiting data points.\nFor each RSL data point, 𝑑, there exists a model prediction of RSL at that location and\n𝑖\n( )\ntime, 𝑚. We wanted to define a misfit function, ρ 𝑑, 𝑚 , so that the total error, 𝐸 , for\n𝑖 𝑖 𝑖 𝑖\none model run is equal to the sum of the misfits between that model run’s predictions of\nRSL and the data at every location/time. Following Menke, 2018 (Geophysical Data\nAnalysis, Eq. 9.6), this can be achieved by relating the misfit to the conditional\nprobability of the data given the model prediction:\nρ = − 2𝑙𝑛𝑝(𝑑|𝑚)\n𝑖 𝑖 𝑖\nSo that:\n𝐸 = − 2∑𝑙𝑛𝑝(𝑑|𝑚)\n𝑖 𝑖\n𝑖\n( )\nIf we take the terrestrial limiting case as an example, we would want 𝑝 𝑚 = 0 when\n𝑖\n𝑑 ≪ 𝑚 (the model predicts that RSL is higher than that inferred from the limiting data\n𝑖 𝑖\n( )\npoint , which is inconsistent) and 𝑝 𝑚 = 1 when 𝑑 ≫ 𝑚 (the model predicts RSL\n𝑖 𝑖 𝑖\nlower than that of the limiting data point, which is consistent). This is consistent with a\nHeaveside function, with the step situated at the model RSL prediction, 𝑚. If we assume\n𝑖\nthat the data points are represented by normal distributions, when 𝑑 ≈ 𝑚, we want\n𝑖 𝑖\n( )\n𝑝 𝑚 to be proportional to the amount of the data’s normal probability distribution that is\n𝑖\n( )\nsituated to the right of this step. This leads us to the following equation for ρ 𝑑, 𝑚 :\n𝑖 𝑖 𝑖\n[ ]\nρ ( 𝑑, 𝑚 ) = 1 + 1 𝑒𝑟𝑓 (𝑞 𝑑 𝑖 −𝑚 𝑖 )\n𝑖 𝑖 𝑖 2 2 𝑛 σ\n𝑖\nFor the terrestrial limiting case, 𝑞 = 1:\n𝑛\n( [ 𝑑−𝑚 ])\n- If 𝑑 ≫ 𝑚, 𝑑 − 𝑚 is large and positive, so 𝑞 𝑖 𝑖 = 1 and\n𝑖 𝑖 𝑖 𝑖 𝑛 σ\n𝑖\n( )\nρ 𝑑, 𝑚 = 1\n𝑖 𝑖 𝑖\n( [ 𝑑−𝑚 ])\n- If 𝑑 ≪ 𝑚, 𝑑 − 𝑚 is large and negative, so 𝑞 𝑖 𝑖 =− 1 and\n𝑖 𝑖 𝑖 𝑖 𝑛 σ\n𝑖\n( )\nρ 𝑑, 𝑚 = 0\n𝑖 𝑖 𝑖\n( [ 𝑑−𝑚 ])\n- If 𝑑 ≈ 𝑚, 𝑑 − 𝑚 is small and 𝑞 𝑖 𝑖 would lie somewhere between 0\n𝑖 𝑖 𝑖 𝑖 𝑛 σ\n𝑖\nand 1 (the closer together that 𝑑 and 𝑚 are to each other, the smaller the value of\n𝑖 𝑖\nthe error function). The probability then becomes dependent on the proportion of\nthe data’s normal distribution lies to the left or right of 𝑚. For example, if the\n𝑖\nmean of the normal distribution is equal to 𝑚, there is an equal proportion of the\n𝑖\nnormal distribution to the left and right of 𝑚 and the probability would be equal\n𝑖\nto ½.\nWe believe that this is too much detail to include in the text but now point the interested\nreader to reference 83 (and references within) for a detailed explanation of this equation\n(line 390-392).\nLine 378: but that those not exclude that lateral variation in the structure are important. You have\nto comment on this despite you use a 1D model to provide context on the limitations of your 1D\nset-up\nWe thank the reviewer for pointing out this missing detail. We have added text and a\nreference to describe our reasoning for excluding both lateral and radial variations in\nviscosity in our LIA GIA models (lines 424-427):\n“We assume a laterally and radially homogeneous viscosity structure because the\ninferred mantle viscosity over centennial timescales varies minimally for the\nupper mantle below Greenland [43] and it is unlikely for small-scale load\nchanges, like those experienced in Greenland since the beginning of the LIA, to\ninvoke deformation within the lower mantle [44].”\nIn addition to justifying our model choice we have also added the caveat that 3D\nvariations may be relevant (lines 427-429):\n“However, we also acknowledge that omitting 3D variations in Earth structure,\nespecially those driven by lateral variations in lithospheric thickness, may have\nsome impact on our results.”\nL395: why is that not affecting the entire south(western) part?\nWe believe that the reviewer meant to ask why the plume isn't affecting the entire\nsoutheastern portion of Greenland. Published estimates of the hotspot track through\nGreenland have the northwest crossing over the plume at approximately 80 Ma and then\ncentral east Greenland passing over it at around 40 Ma (see Fig 1 in ref 25 and Fig 2 in\nref 86). The two stations in central East Greenland that we removed have anomalously\nhigh uplift rates, suggesting that the mantle below them is still weak from passing over\nthe hotspot, a process that we do not account for in our GIA models (which is why we\nremove them from our analysis). The weakening of mantle structure driven by the hotspot\nis thought to be relatively narrow (ref. 87), so we do not believe the entire southeastern\npart of Greenland to be affected.\nL422: I would expect V_RF to be independent of time, but you suggest you use each of the\n14,040 predictions?\nIt is correct that the correction is independent of time. V_RF is the contribution to VLM\ndriven by differences between the reference frame of the GIA model (center of mass of\nthe solid Earth reference frame) and that of the observations. Because each of the 14,040\nGIA simulations that we use distributes mass within the solid Earth slightly differently\n(depending on Earth structure and ice history parameters), the center of mass of the solid\nEarth that the prediction is in reference to will be slightly different between each of the\nmodel runs, and thus we must make this correction for each.\nWe thank the reviewers again for this very detailed and constructive review. Our appreciation is\nreflected in the acknowledgements."}, {"round": 2, "reviewer_comments": {"reviewer_1": "We appreciate the improvements made by the authors and respect their choices, nevertheless we add an annotated\nmanuscript with what we consider minor suggestions for further improvement. We don't need to review them any further.", "reviewer_2": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature\nCommunications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career\nResearchers who co-review manuscripts.", "reviewer_3": "I am overall ok with the review and I don't have further comments", "reviewer_4": "I have reviewed the authors' updates to the manuscript, and commend them for considerable improvements, especially to\nthe readability of some figures that were previously difficult to distinguish, and for the inclusion of important clarifications in\nthe text.\nHowever, I do still have one major concern, which again relates to the way that sea-level projections from AR6 are being\nincorporated.\nIn response to my original comment detailing concerns with the authors' approach, the authors said that \"We have plotted\nthe percentiles provided by the AR6 projections and those associated with our approximated mean and standard deviation\nfor all processes at select times and locations, and find that they fit well.\" They also provided a plot of both AR6 data for the\nocean dynamics component of sea-level rise and their own analysis, assuming a normal distribution.\nHowever, for the ocean dynamics component for the AR6 projections, a t-distribution with the mean/co-variance of a multi-\nmodel ensemble constructed from CMIP6 models was used. A normal distribution fit to the ocean dynamic component of the\nprojections would therefore be expected to align very closely, with the greatest discrepancies being in the tail, owing to the\nrelatively few CMIP6 models available (~20), for which the t-distribution accounts for with a degree-of-freedom term. (see\nFox-Kemper et al., [Ref. 1 in the manuscript], Kopp et al., 2023, [Ref. 50 in the manuscript] Kopp et al., 2014 [\nhttps://agupubs.onlinelibrary.wiley.com/doi/10.1002/2014EF000239 ]).\nThis does not mean that it is appropriate to assume a normal distribution for all components or for the full SLR projections,\nwhich are well-known to have long (and important) upper tails. The AR6 datasets that the authors use are provided as full\ndistributions for a reason. The appropriate methodology would require working with the full distribution and framing the\nauthors' results in terms of that distribution. The current assumption of a normal distribution for sea-level rise projections is\ninappropriate--particularly given that sea-level rise projections are well-known not to conform to a normal distribution.\nBuilding upon my concerns about this methodology, I continue to be concerned about the validity of the uncertainty that the\nauthors provide for their total estimate of RSL around Greenland for the 21st century. As noted above, the correct approach\nfor using the AR6 projections is not to treat them as normal distributions, but to work with the full distribution. If an approach\nthat uses the full distribution is adopted (as it should be), the analysis of uncertainty will also need to be adjusted, since, as\nnoted in my original response, uncertainty across percentiles cannot just be added together (unless it is the median)."}, "author_response": "We thank the reviewers for their positive assessment. We have addressed all minor suggestions\nbelow.\nLine-by-Line Comments\nL21: why not immediately mention SSPs here?\nThe contributions to RSL around Greenland from thermal expansion and ocean dynamics, land\nwater storage, and global glacier change are based on SSP-RCP scenarios (i.e SSP1-2.6),\nhowever, the ice sheet projections we use are based solely on the emissions scenarios RCP 2.6\nand RCP 8.5. We refer to the RCP emission scenario when we talk about total projected RSL, as\nonly some of the contributions are based on SSP scenarios. We specify this in the revised text\n(lines 265-267):\n“Note that here and for the remainder of the text, we refer to the RSL projections as\nrepresenting emissions scenarios RCP 2.6 and 8.5, which include the SSP1-2.6 and SSP5-8.5\nprojection components for the other processes.”\nL22: not clear from the rebuttal where the change from 30-40 to 20-40 comes from?\nWe misread the result of the relative contribution of viscous deformation to RSL at 2100 CE (Fig\n8d) in the original submission, which led to this correction in the revised manuscript.\nL66-87: This paragraph assumes that misfits can only arise from GIA, but uncertainties in ice\ndynamics and climate forcing maybe equally well be considered to be responsible for misfits.\nWe thank the reviewers for pointing out this nuance. We have included a sentence that says\ndata-model residuals can come from uncertainties in both the ice reconstruction and the Earth\nstructure (lines 71-73):\n“Each of these studies finds that their predictions of RSL are unable to match all of the\nobserved RSL data, a consequence of uncertainties in both the ice reconstruction and\nprescribed Earth parameters.\nL79: The previous sentences contain a very nice overview of the different viscosities that other\nstudies found, but it doesn’t state anything about the misfits of those studies. Please add\ninformation on the misfits of other studies before stating in the highlighted part that improvement\nof the misfit is needed.\nWe have added a sentence that mentions how each of the studies we described that tunes their\nmodel simulations to fit paleo RSL data has remaining misfits to that data (line 71-73, see\nprevious comment). We believe a more in-depth description of the residuals that exist for each\nstudy would dilute the main point of the paper. We highlight that we use published model runs\nfor this component and it only contributes minorly to future sea level change around Greenland.\nWhile important to include, this component is not the main focus of the paper.\nL85: Not clear where in the text these were indicated before.\nAs described in the response to the review for L66-87, we have now added a sentence where\nthese misfits are explicitly mentioned (lines 71-73).\nL85: This document does not contain a reference list, but according to the reference list in the\nnew manuscript without tracked changes, this refers to Bevis et al., PNAS, 2012. I cannot find a\ncomparison of 1D and 3D GIA models in this paper, is this reference correct?\nWe thank the authors for pointing out this reference error. It has been corrected in the revised\nmanuscript.\nL86: We have commented about this sentence before and the sentence is improved by replacing\n\"beyond paleo RSL data\", but it is still not connected to the rest of the paragraph. Please discuss\nbefore this sentence what the previous studies did incorporate to constrain the mantle viscosity\nand whether or not they incorporated present-day observations. Only after this information you\ncould state that incorporating present-day observations could further constrain viscosity.\nIn lines 66 and 68-69, we specify that the study from Lecavelier et al. (2014), Tarasov &\nPeltier (2002), and Fleming et al. (2004) tune their GIA models to paleo RSL data to constrain solid\nEarth parameters. However, we do agree that the last sentence did seem disconnected from the\nrest of the paragraph and have removed it from the revised manuscript.\nL117: Reference error.\nThis has been corrected in the revised manuscript.\nL146: explain better why you use 1D and 3D. Why is 1D sufficient for changes since LIA,\nwhy not use 3D as well?\nWe already provided a justification for using a 1D GIA model for the LIA instead of a 3D one in\nour last version, which is given in lines 143-146 of the revised manuscript.\nL158: Why RCP and not SSP?\nThe ice-sheet projections considered in this study were taken from ISMIP6 simulations forced by\nCMIP5 climate models, which only represent the RCP scenarios and not the combined SSP-RCP\nscenarios. Our reasoning for using the CMIP5 forced ISMIP6 simulations is detailed in the\nresponse to the comment made on L382 below.\nL185: Unclear what the basis is for the weighing presumably misfit with the observational record\nbut you don't mention this.\nIn lines 128-130, where we give a brief description of how we calculated the deglacial\ncontribution to both present-day VLM and future sea-level change, we write that we weight the\nmodels based on their fit to Holocene RSL data. This is then further detailed in the methods\n(lines 383-398).\nL196: whereas ICE6G is .....\nWe have changed the sentence to provide a clear comparison between ICE-6G and Huy3 (lines\n173-175):\n“This is expected, as additional and improved data constraints were used in the creation\nof of Huy3 compared to the GrIS portion of ICE-6G.”\nFigure 2: Include in the caption over which time period the mean in taken.\nWe believe that including this in the caption would over-complicate the figure caption as the time\nperiod over which the mean is calculated differs for each site and also slightly for each\ncontribution (i.e., the elastic and observed, the deglacial GIA, and the LIA contributions).\nL352: A reference to new supplemental Fig. 4 is missing. The caption of that figure also misses a\nreference to where the figure was taken from I think.\nSupplementary Figure 4 is referenced in lines 270-271. We have now included the reference\nwhere we obtained the individual ice sheet projections from in the caption.\nL382: is there CMIP6 forced ISMIP, if so why did you not include that?\nIn order to account for the uncertainty in LIA ice history and earth structure, we could only run a\nsubset of the future ice projections considered in the IPCC AR6. Most of the ensemble of CMIP6\nmodels have Effective Climate Sensitivities (ECS; how much global mean temperature changes\nwhen carbon dioxide concentration in the atmosphere is doubled) that range from 2.3°C-3.2°C or\n4.6°C-5.6°C. The CMIP5 models used in Goelzer et al. (2020), from which we obtained our\nice-sheet projections, have ECS that are either close to the median (2.7°C-2.9°C) or are from the\nupper end (4.1°C-4.7°C) of ECS from the ensemble of CMIP5 models. So the ice-sheet\nprojections used in our study were forced by CMIP5 models with similar ECS as that in many of\nthe CMIP6 models. However, we don't represent the higher end climate sensitivity from CMIP6\nmodels (ECS > 4.7°C), so our estimates of RSL at 2100 CE might be more conservative, as it is\nshown that those CMIP6 models do produce a greater GMSL contribution from Greenland\n(Payne et al., 2021). We have added a sentence to the manuscript to describe this limitation (lines\n306-310):\n“The range of climate sensitivities represented in the ensemble of CMIP6 models has an\nupper tail that exceeds the range of climate sensitivities from the CMIP5 models used to\ncreate the ice-sheet projections considered in this study. Therefore, the ice-sheet\nprojections used here represent a more conservative-to-moderate subset of those\nconsidered in AR6.”\nL404: layout mistake\nWe thank the reviewer for pointing this out and have corrected it in the revised manuscript.\nReferences\nPayne, Antony J., et al. \"Future sea level change under coupled model intercomparison project\nphase and phase 6 scenarios from the Greenland and Antarctic ice sheets.\"\nGeophysical Research Letters 48.16 (2021): e2020GL091741.\nReviewer #4 Evaluations:\nI have reviewed the authors' updates to the manuscript, and commend them for considerable\nimprovements, especially to the readability of some figures that were previously difficult to\ndistinguish, and for the inclusion of important clarifications in the text.\nWe thank the reviewer for this positive assessment.\nHowever, I do still have one major concern, which again relates to the way that sea-level projections from AR6 are being\nincorporated.\nIn response to my original comment detailing concerns with the authors' approach, the authors\nsaid that \"We have plotted the percentiles provided by the AR6 projections and those associated\nwith our approximated mean and standard deviation for all processes at select times and\nlocations, and find that they fit well.\" They also provided a plot of both AR6 data for the ocean\ndynamics component of sea-level rise and their own analysis, assuming a normal distribution.\nHowever, for the ocean dynamics component for the AR6 projections, a t-distribution with the\nmean/co-variance of a multi-model ensemble constructed from CMIP6 models was used. A\nnormal distribution fit to the ocean dynamic component of the projections would therefore be\nexpected to align very closely, with the greatest discrepancies being in the tail, owing to the\nrelatively few CMIP6 models available (~20), for which the t-distribution accounts for with a\ndegree-of-freedom term. (see Fox-Kemper et al., [Ref. 1 in the manuscript], Kopp et al., 2023,\n[Ref. 50 in the manuscript] Kopp et al., 2014 [\nhttps://agupubs.onlinelibrary.wiley.com/doi/10.1002/2014EF000239 ]).\nThis does not mean that it is appropriate to assume a normal distribution for all components or\nfor the full SLR projections, which are well-known to have long (and important) upper tails. The\nAR6 datasets that the authors use are provided as full distributions for a reason. The appropriate\nmethodology would require working with the full distribution and framing the authors' results in\nterms of that distribution. The current assumption of a normal distribution for sea-level rise\nprojections is inappropriate--particularly given that sea-level rise projections are well-known not\nto conform to a normal distribution.\nBuilding upon my concerns about this methodology, I continue to be concerned about the\nvalidity of the uncertainty that the authors provide for their total estimate of RSL around\nGreenland for the 21st century. As noted above, the correct approach for using the AR6\nprojections is not to treat them as normal distributions, but to work with the full distribution. If\nan approach that uses the full distribution is adopted (as it should be), the analysis of uncertainty\nwill also need to be adjusted, since, as noted in my original response, uncertainty across\npercentiles cannot just be added together (unless it is the median).\nWe thank the author for this detailed comment and agree that using the full distributions instead\nof approximate normal distributions would improve the accuracy of our results. We have now\nchanged our approach to capture the full distribution of each contribution. We calculate the total\nsea level projection by sampling the distribution of each contribution (at each location and\ntimestep) and summing it together. This results in a distribution for the total future sea level\nchange of which we report the median and likely range (which we define as the difference\nbetween the 17th-83rd percentiles, following the IPCC AR6). We also plot the 5th-95th ranges of\nthe RSL distributions for each contribution and the total at 2100 CE for Tasiilaq and Nuuk to\nshow the spread between the tails of our distribution (Fig. 6e-f).\nThese changes are propagated throughout the revised manuscript (e.g. see Figs. 4-6) and the\nchanges in our methods are described in lines 398-399 and 571-582. We believe that our\nprojections now incorporate the IPCC projections sufficiently well and analyze uncertainty as\nsuggested, and again, thank the reviewer for their guidance on this methodology."}, {"round": 3, "reviewer_comments": {"reviewer_4": "I appreciate and commend the effort the authors have put into revising their methodology for the sea-level rise projections. I\nhave no further comments at this time, and would recommend that it move forward for publication in the present form.\n~Andra Garner"}, "author_response": ""}]

Decision letter