A 21st Century Warming Threshold for Sustained Greenland Ice Sheet Mass Loss

Under anticipated future warming, the Greenland ice sheet (GrIS) will pass a threshold when meltwater runoff exceeds the accumulation of snow, resulting in a negative surface mass balance (SMB < 0) and sustained mass loss. Here, we dynamically and statistically downscale the outputs of an Earth system model to 1 km resolution to infer that a Greenland near‐surface atmospheric warming of 4.5 ± 0.3°C—relative to preindustrial—is required for GrIS SMB to become persistently negative. Climate models from CMIP5 and CMIP6 translate this regional temperature change to a global warming threshold of 2.7 ± 0.2°C. Under a high‐end warming scenario, this threshold may be reached around 2055, while for a strong mitigation scenario it will likely not be passed. Depending on the emissions scenario taken, our method estimates 6–13 cm sea level rise from GrIS SMB by the year 2100.

significantly exceeding total snow accumulation (e.g., 2012 and 2019), SMB has remained positive (Tedesco and Fettweis, 2020) because a large fraction (∼45%) of the meltwater is retained in the firn, the up to 100 m thick layer of compressed, perennial snow covering ∼90% of the ice sheet. In the firn layer, meltwater either refreezes as ice lenses/slabs (Machguth et al., 2016;MacFerrin et al., 2019) or is retained as liquid water in perennial firn aquifers (Forster et al., 2013).
Simply extrapolating the recent negative trend into the future predicts the SMB = 0 threshold to be reached within ∼20 years ( Van den Broeke et al., 2016). However, such predictions are uncertain given the relative brevity of the time series (Hanna et al., 2020) and the large interannual SMB variability (Wouters et al., 2013). Previous studies have relied on semi-empirical schemes such as positive-degree-day models (Gregory & Huybrechts, 2006) and idealized energy-balance models (Robinson et al., 2012) to project surface runoff and SMB into the future. Here, we present a novel physically based approach, using explicitly resolved surface energy balance from a state-of-the-art regional climate model (RACMO2.3p2) to dynamically downscale a CMIP6 projection from the Community Earth System Model (CESM2) under a high-end climate warming, further statistically downscaled to 1 km resolution. Uniquely, without any bias correction in the forcing, this physically based model setup reproduces a realistic present-day GrIS SMB , and shows high correlation with the native CESM2 forcing, enabling us to project meaningful GrIS SMB under various CMIP6 warming scenarios after applying a small linear correction. Since GrIS geometrical changes and the associated (elevation-temperature) feedback are expected to remain small during the 21st century (Aschwanden et al., 2019;Edwards et al., 2014;Le clec'h et al., 2019), first-order inferences about the SMB = 0 threshold can be made assuming constant ice sheet geometry.

Methods
Our main GrIS SMB product was obtained by dynamically downscaling the Community Earth System Model version 2.1 (CESM2 henceforth; Danabasoglu et al., 2020) over the period 1950-2100 using the regional climate model RACMO2.3p2 (RACMO2 henceforth; Noël et al., 2018). The most extreme warming scenario SSP5-8.5 was selected to obtain maximum leverage for correlations between SMB and regional temperature T GrIS (see supporting information). To better resolve the high ablation rates over low-lying, narrow outlet glaciers that are primary contributors to meltwater runoff, RACMO2 SMB components were further statistically downscaled from 11 to 1 km resolution and the resulting product is referred to as "RACMO2 SSP5-8.5." Statistical downscaling corrects daily melt and runoff for elevation bias on the relatively coarse 11 km grid using elevation gradients, and for underestimated bare ice albedo using remote sensing records (see supporting information). Table A1 provides an overview of the various data products, including a benchmark reanalysis-based RACMO2 simulation (RACMO2 ERA henceforth) used for model evaluation and discussed in Noël et al. (2019).
Both RACMO2 (native resolution ∼11 km; Noël et al., 2018) and CESM2 (native resolution ∼111 km; van Kampenhout et al., 2020) have been extensively evaluated over the GrIS and compare well to in-situ measurements, although a better agreement is generally found for the high-resolution RACMO2 than for the relatively coarse CESM2 Van Kampenhout et al., 2020). GrIS-integrated SMB, precipitation and runoff correlate well between the two models (R 2 = 0.96, Figure S1). The limited but persistent deviations are mostly resolution-related. For instance, too gentle coastal ice sheet slopes in the low-resolution model topography of CESM2 allow moist air masses to propagate too far inland, which results in overestimated precipitation in the GrIS interior (Van Kampenhout et al., 2019). The inferred regression statistics ( Figure S1) are then used to correct SMB, precipitation, and runoff fields for all available CESM2 ensemble members (see Table A1, column "Corrected SMB"), enabling us to reconstruct RACMO2-like SMB time series for other warming scenarios (see supporting information).

Results
The resulting time series are presented in Figure 1a for SMB and in Figure S2 for precipitation and runoff. Note that CESM2 is only forced with greenhouse gas and aerosol concentrations and does not assimilate observations, so that period averages, average interannual variability and trends can be compared to NOËL ET AL. 10.1029/2020GL090471 2 of 9 observations but not individual high/low SMB years. In spite of this, the historical ensemble members (gray; 1950-2014) agree well with the RACMO2 ERA SMB product (green line, see Table A1), including the relatively stable precipitation ( Figure S2a), the recent increase in runoff ( Figure S2b), and the subsequent decrease in SMB (Figure 1a), resulting in an almost seamless transition into the subsequent scenario runs.
In all scenarios, the overall SMB signal is dominated by greater runoff ( Figure S2c) leading to an SMB decrease and surface-driven GrIS mass loss ( Figure 1a). Most runoff (∼85% of the total) originates from the narrow and widening marginal ablation zone (Noël et al., 2019). In combination with the hypsometry of the ice sheet, which flattens inland, this leads to a stronger-than-linear increase in runoff with increasing temperature ( Figure S2b) and, since precipitation increases only linearly or not at all ( Figure S2a), leads to a nonlinear decrease in SMB ( Figure 1a). We find SSP5-8.5 as the only scenario in which precipitation increases significantly after 2015 (1.1-2.3 Gt yr −2 , p < 0.05, Figure S2a). The increase in rainfall (1.4-1.6 Gt yr −2 , p < 0.05) comes at the expense of snowfall (−0.3 to −0.6 Gt yr −2 , insignificant p > 0.05) with the rainfall fraction increasing at a rate of ∼0.15% yr −2 .
We consider the SMB = 0 threshold to be exceeded when the five-year moving average SMB becomes negative. Figure 1b displays five-year-averaged SMB as a function of the near-surface Greenland temperature anomaly (T GrIS ) with respect to preindustrial conditions (see supporting information). The black points representing RACMO2 SSP5-8.5 roughly fall onto a single quadratic curve (R 2 = 0.90) described by, allowing us to determine a regional warming threshold for which GrIS SMB becomes negative: 4.5 ± 0.3°C (vertical red line in Figure 1b). The uncertainty of 0.3°C is based on the model SMB uncertainty of 48 Gt yr −1 derived from a comparison between SMB from in situ measurements and from the historical CESM2-forced RACMO2 simulation at 1 km ) (see supporting information). The statistical uncertainty of the quadratic fit is comparatively small (<0.1°C calculated from a 68% confidence band) and therefore neglected. Furthermore, these results are robust under inclusion of the reconstructed CESM2 scenarios (colored dots in Figure 1b), which suggests that the nonlinear dependence of SMB on T GrIS is independent of the actual scenario or ensemble member used (see supporting information). Note that the apparent discrepancy between RACMO2 and the CESM2 members under SSP5-8.5 in 2095-2100 (red and black dots in Figures 1a and 1b) is due to internal variability in CESM2. Specifically, the CESM2 "parent" member forcing RACMO2 (not shown) has lower SMB than the two independent members in the figure.
(a) (b) Figure 2 shows a map of 5-year average specific (local) surface mass balance (units kg m −2 yr −1 or mm w.e. (water equivalent) yr −1 ) that represents the situation at the 4.5°C warming threshold, that is, when SMB ≈ 0. Ablation rates at the low-lying GrIS margins are enhanced by high atmospheric temperatures and the low albedo of bare ice that is exposed in the ablation zone during summer (Box et al., 2012), and in magnitude exceed the interior accumulation rates (Figure 2). The ablation zone area increases from an average of 8 % of the contiguous GrIS area in the period 1961-1990 to 21 % in 2051-2055 ( Figure S3), when the SMB = 0 threshold is reached. The strongest expansion occurs in the relatively dry southwest, from 13% to 32% and in the northwest GrIS, from 7% to 19%, in line with recent observations of bare ice zone expansion (Noël et al., 2019).
Under an SSP5-8.5 warming scenario, the equilibrium line altitude rapidly migrates inland, expanding the GrIS ablation zone area in a nonlinear fashion throughout the 21st century ( Figure S3), in line with the RACMO2 ERA simulation (green line). Exposure of bare ice in growing ablation zones triggers a strong melt-albedo feedback, amplifying runoff ( Figure S2b) and driving the nonlinear SMB decline (Figure 1a) as atmospheric temperature increases (Figure 1b).

Discussion
Previous work by Gregory and Huybrechts (2006) used contemporary climate model output and a degree-day model for ablation and found a regional warming threshold of 4.5 ± 0.9°C, similar to our results, and a global warming threshold of 3.1 ± 0.8°C. Using a linear regression between SMB from the Modèle Atmosphérique Régional (MAR) and annual mean T GrIS from multiple General Circulation Models, Rae et al. (2012) found that a 2.3 ± 0.4°C warming relative to 1980-1989 is required for the GrIS to reach the SMB = 0 threshold. In addition, Fettweis et al. (2013) used a cubic fit between MAR SMB and annual mean global temperature to show that the GrIS would pass the SMB = 0 threshold for a 3°C warming relative to 1980-1989. Here we refine these estimates by combining recently developed state-of-the-art modeling tools. However, the temperature threshold presented in this paper should be considered an upper bound for two reasons. First, SMB = 0 does not consider marginal retreat of the ice sheet, a response that operates at time scales greater than the ∼100 years considered here. By coupling an intermediate complexity climate model with a dynamical ice sheet model, Robinson et al. (2012) found that the GrIS may already tip towards deglaciation at global temperature perturbations in the range 0.8-3.2°C, with SMB initially above zero. Second, the ice sheet geometry in RACMO2 is held constant throughout the simulation, thereby ignoring dynamical thinning and the positive melt-elevation feedback. Fettweis et al. (2013) estimated that this feedback could cause an additional surface mass loss of 8 ± 5% towards the end of the 21st century, implying a lower temperature threshold for SMB = 0 even on short time scales.
It is unclear by how much T GrIS has already increased since the preindustrial period. Observations from (mainly coastal) meteorological stations indicate a significant warming since the mid-1990s of ∼5°C in winter and ∼2°C in summer (Hanna et al., 2012). Here we consider ensemble-averaged T GrIS anomalies from the CMIP5 (Figure 3a) and CMIP6 (Figure 3b) climate model archives, with the uncertainty band representing the intermodel standard deviation (see Figure S4). T GrIS from a hybrid ECMWF reanalysis product (Table A1) is scaled to the same baseline period (see supporting information) and is used as a proxy for observations, with the remark that reanalysis temperature products are not well constrained over the GrIS because of poorly represented snow/ice processes in the used land model and the scarcity of near-surface observations to assimilate. NOËL ET AL.  , separating the accumulation and ablation zones.
In agreement with state-of-the-art GrIS SMB estimates, the SMB = 0 warming threshold (indicated by the horizontal red bar, Figure 3) has been reached in neither the reanalysis nor CMIP models up until today. Looking towards the future, the CMIP5 RCP8.5 ensemble-mean crosses the warming threshold in 2056, with a 66% likely range of  with similar estimates for CMIP6 (mean 2060, likely range ). However, for the medium warming scenario (RCP4.5) the CMIP5 and CMIP6 archives disagree on the year of crossing. Whereas the CMIP5 ensemble mean crosses the threshold in 2125, the CMIP6 mean crosses four decades earlier in 2084, which could be attributed to the higher climate sensitivities of some CMIP6 models for the same amount of radiative forcing (Zelinka et al., 2020). For both ensembles, a likely range could not be calculated due to the flatness of the curves (Figure 3). Finally, in the high-mitigation scenario RCP2.6 the ensemble-mean remains below the warming threshold for the entire time period considered (Figure 3). An overview of estimated year of crossing for all scenarios is provided in Table S1.
Using the CMIP ensembles, we can translate the regional threshold of 4.5 ± 0.3°C into a global temperature warming threshold. Figure 4 indicates near perfect correlation (R 2 > 0.99) between T GrIS and global average surface temperature, T global , in the ensemble mean. The slope of the correlation yields a Greenland warming amplification of ∼1.6 in both CMIP5 and CMIP6, in line with a previous CMIP3 estimate of 1.5 ± 0.4 K K −1 (Gregory & Huybrechts, 2006). Results by Frieler et al. (2011) suggest that these results are robust when individual models are used, rather than the ensemble mean. From the linear relationship in Figure 4 we NOËL ET AL.
10.1029/2020GL090471 5 of 9 Figure 3. Multimodel estimates of regional Greenland warming under multiple scenario projections. Five-year simple moving average of the anomaly (relative to 1850-1899) of near-surface temperature averaged over Greenland (T GrIS ). Presented are the model ensemble means (solid lines) and model ensemble standard deviation (shading) for all available CMIP5 (a) and CMIP6 (b) models, and three emission scenarios. Number of models included is listed in the legend as n. The green solid line represents a hybrid reanalysis temperature product for comparison, which has been aligned to the CMIP data over the period 1958-1977 (see supporting information) and the inferred GrIS warming threshold is shown as a horizontal red bar.

(a) (b)
derive a T global warming threshold of 2.7 ± 0.2 C, with respect to the preindustrial, as a precursor for GrIS deglaciation under sustained warming.
It is important to note that the presented method is only valid for a fixed ice sheet geometry, an assumption that becomes questionable beyond 2100 (Fettweis et al., 2013). On these longer time scales, projections require a coupled atmosphere-ocean-ice sheet model system (Le clec'h et al., 2019). For the short term, however, we can apply the quadratic SMB − T GrIS anomaly relation (Equation 1) to the CMIP5 and CMIP6 model ensembles to reconstruct 21st century SMB ( Figure S5a) and global mean sea level rise due to GrIS surface mass loss ( Figure S5b). Obviously, the largest sea level rise contributions are predicted for the strongest warming scenarios (∼13 ± 4 cm in 2100 for RCP8.5, Figure S5b). Substantial contributions from GrIS surface mass loss are still expected for the strong mitigation scenarios (∼6 ± 3 cm in 2100 for RCP2.6) since persistent GrIS mass imbalance remains even when the warming threshold is not reached. For 2081-2100, Table S2 summarizes our global mean sea level estimates for the different scenarios and compares them against existing IPCC values. Our median estimates are consistently higher than those of IPCC, but the likely ranges have a large overlap. In summary, our study uses stateof-the-art SMB modeling tools for the GrIS to show that the temperature threshold for SMB = 0 is likely to occur in all but the strong mitigation scenarios, and is particularly likely to occur in the upper-end warming scenarios.

Conclusions
In this study we explored possible pathways of future GrIS SMB by combining Earth system model output with dynamical and statistical downscaling techniques. For SMB = 0, we found a regional warming threshold of 4.5 ± 0.3°C compared to preindustrial, which translates to a global warming threshold of 2.7 ± 0.2°C using climate model output from CMIP5 and CMIP6. Under high-end warming scenario RCP 8.5, SMB is projected to become persistently negative around 2055, while for the strong mitigation scenario RCP 2.6 the threshold for sustained mass loss is likely not reached. Our method estimates 6-13 cm of sea level rise by 2100, depending on the emissions scenario. Assuming that future sea level rise is to be curbed, our results underscore the need for a rapid stabilization of greenhouse gas concentrations, especially since our estimates are conservative. NOËL ET AL.
10.1029/2020GL090471 6 of 9  Yes, see Figure S1 van Yes, see Figure S1 Idem CMIP5 RCP2.6 1850-2100 ∼1 degree CMIP5 multimodel mean 27 models, see Table S3 --  Table S4 - Note. Not all CMIP models start exactly in 1850, which means a slightly different reference period may have been used for calculating anomalies (see supporting information).