Long-term changes in ice discharge and comparison with other studies
We find a step-increase in decadal-scale ice discharge (Fig. 1a), with a ~60 Gt yr−1, or 14%, increase between 1985–1999 and 2007–2018 means. After reaching a temporally local maximum in 2005, annual D then temporarily decreased for 3 years. Following the temporary decline, discharge accelerated again at a slower pace of 2 Gt yr−2 during 2008–2018, reaching a peak annual value of 502 ± 9 Gt yr−1 in 2017 and 2018, or 17% above the 1980’s average. The increase in mean annual D since 2008 has been mostly due to a steady increase in seasonal minimum values increasing with a trend of 3 Gt yr−2 since 2007, indicating greater wintertime velocities relative to summertime maxima, most evident in the northwest (Supplementary Fig. 1) and in the most recent 3 years of the central west. The seasonal amplitude in D has also changed, increasing by nearly 50%, from a 1985–1990 average of 17 ± 6 to 25 ± 6 Gt yr−1 for 2000–2018. To account for the uncertainty in D due to this temporal gap in ice thickness observations, we estimate D assuming the end member-cases of (1) all thickness change occurring in the first year, which maximizes the impact of thinning at the start of the period, and (2) all thinning occurring in the last year, which minimizes the impact until ~2000. We find that during the 1985–1999 period, estimates of D can vary by an average 13 Gt yr−1 (Fig. 1a) depending on when thinning occurred between temporally sparse elevation data (AeroDEM, ~1985 and ASTER, nominally ~2000), described in more details in the Methods section.
Discharge has increased in every region (defined in Supplementary Fig. 2) of the ice sheet since the 1980’s, but with contrasting temporal variability. In the west, discharge slightly declined from 1985 to the late 1990’s prior to larger increases after 2000. The central west (CW; Supplementary Fig. 1c) is dominated by Jakobshavn Isbræ, which accounted for an average of 30% of the regional discharge. Within the CW, discharge steadily increased by 27% after the initial retreat and acceleration of Jakobshavn Isbræ in 19999 to a peak in 2014 before declining 10% through 2018 due to Jakobshavn’s slowdown10,11. Jakobshavn’s observed slowdown was responsible for the majority (>75%) of the net decrease in D in the CW region, with only two other glaciers showing notable deceleration during this time. In the northwest (NW; Supplementary Fig. 1a), gradual increases in discharge began in the early 2000’s and have since accelerated, increasing 36% by 2018. In the east, discharge was relatively stable during the 1980’s and 1990’s, with increases beginning in the southeast (SE; Supplementary Fig. 1d) in 2000 with the synchronous retreat and acceleration of multiple glaciers, including Helheim12,13. After increasing ~18% between 2000 and 2005, discharge in the SE then suddenly declined to within 5% above its initial [1980s mean] values, varying by 10% over the next decade with no clear trend. Pronounced acceleration did not begin in the central east (CE; Supplementary Fig. 1b) until the 2004 retreat and acceleration of its largest glacier, Kangerdlugssuaq, that caused a nearly 20% increase in annual D between 2003 and 2005. Following a decline and relative stability, another large retreat event at Kangerdlugssuaq in 201614 caused D to again rise near its 2005 maximum. We exclude discussion of regional patterns observed in the northern and southwestern regions due to small glacier sample size, but include regional discharge totals (with each of the two regions contributing, on average, ~10 Gt yr−1) in Supplementary Fig. 1e, f.
We compare our results with those presented in ref. 5 which is the only other study, to our knowledge, that estimates ice sheet-wide discharge at sub-annual resolution for a similar time period shown here. Our results closely agree with those of ref. 5 with a mean difference <4 Gt yr−1 (~1%) and absolute biases of <8 Gt yr−1 (Supplementary Fig. 3), and summertime estimates within 2 Gt yr−1. Larger wintertime differences are within the uncertainty expected from differences in interpolation methods (i.e., Kalman filter model used in this study and linear interpolation done in ref. 5). While GrIS-wide biases are reduced by regionally compensating errors, we find that the two studies still closely agree on the regional scale, with mean temporal biases of <10% for all major drainage basins. This level of agreement is noteworthy given that the two studies use independently positioned flux gates and different velocity and elevation datasets. Our annual D estimates are also in general agreement with those from ref. 4 prior to 2002, but later diverge, with the estimates from ref. 4 exceeding ours, and those in ref. 5, by >40 Gt yr−1 (8%) in 2018, excluding peripheral glaciers and ice caps (Supplementary Fig. 3). We attribute this difference to the treatment of glaciers with missing bed topography data and how ice velocity is sampled, as described in the Supplementary Discussion (see also Supplementary Fig. 4).
Decadal variability in components of ice sheet mass balance
Combining our estimates of D with SMB from RACMO2.3p215, including the peripheral glaciers and ice caps, gives a cumulative net loss of 4200 Gt between 1985 and 2018 (Fig. 1c). The average loss rate of −240 Gt yr−1 between 2003 and 2018 is close to the approximate −235 Gt yr−1 rate estimated from GRACE for the same period (update of ref. 16). From 1985 to 1999, annual D consistently exceeded the mean 1961–1990 GrIS-wide SMB, typically assumed near equilibrium17,18, both with (424 Gt yr−1) and without (409 Gt yr−1) the inclusion of peripheral glaciers and ice caps. During this period SMB fluctuated around equilibrium, so that discharge contributed ~70% to the cumulative ice sheet mass loss of ~400 Gt prior to 2000 (Fig. 1c). After 1997, however, SMB anomalies became increasingly negative at a rate of nearly −20 Gt yr−2 for the next 15 years, reaching approximately −300 Gt yr−1 in 2012. This resulted in discharge contributing only approximately one-third of net mass losses between 2000 and 2012. Discharge has contributed to >50% of cumulative mass losses since 2013, however, as SMB anomalies have been less negative.
For the ice sheet to gain mass, the total annual SMB must exceed the total D. Assuming the mean 1985–2000 D of 426 ± 6 Gt yr−1, the mass gain from SMB offset the loss due to D, resulting in ice sheet growth, nearly 40% of the years between 1985 and 1999 (Fig. 1b), which is similar to estimated mass gain estimated in one-half of the years during the 1992–1999 period in ref. 1. Increasing D to the current, post-retreat average of ~495 Gt yr−1 cuts the percentage of years of net gain to one in five, demonstrating that the increase in discharge alone was sufficient to transition the ice sheet to a new state of persistent negative mass balance. However, SMB has also markedly declined and, assuming the post-2000 average SMB, the annual probability of net ice sheet mass gain is now only ~1%. Thus, the combined climate and dynamic states of the ice sheet are such that even a single year of mass gain is highly unlikely.
Long-term relationships between retreat and ice discharge
King et al.3 found interannual variability in D after 2000 to be primarily associated with changes in ice front position, likely due to acceleration caused by the perturbation to the glacier stress regime during retreat19. Here we assess changes in front position and discharge since 1985, capturing the period prior to widespread retreat and acceleration. In order to compare ice front retreat and discharge, we use a subset of glaciers (n = 128) that have continuous records of discharge and front position change. In order to assess the impact of mean front position change on the regional and ice sheet discharge and minimize the impact of small glaciers, we weight the front position of each glacier by its discharge and denote the weighted mean front position change as Fw. Figure 2a shows that GrIS D increased at an average rate of 14 ± 1 Gt yr−1 per km of cumulative retreat over the 1985–2018 period, with the cumulative change in front position explaining over 93% of the change in D. This plot reveals two distinct clusters of values, before the year 2000 and after the year 2009, separated by roughly 4 km of retreat and 50 Gt yr−1. D increased most rapidly during widespread acceleration and retreat between 2001 and 2005, as glaciers in the SE rapidly retreated off bathymetric highs13, at a consistent rate of 16 Gt yr−1 per km of cumulative retreat. D increased at a significantly lower rate of 8 Gt yr−1 per km of retreat following the readvance in 2006 and 2007. Since 2014, there has been little additional change in Fw or D for the ice sheet as a whole. Examining year-to-year changes in the Fw and D in Fig. 2b, we find a remarkably consistent relationship with annual changes in discharge and front position, averaging 14 ± 6 Gt yr−1, or ~4%, per km of retreat, with annual variability in the front accounting for 41% of the annual variability in D. Further, we find that the change in D per km of Fw is also consistently 4–5% from region to region (Supplementary Fig. 5), indicating a similar overall sensitivity of discharge to retreat.
Since changes in discharge are dependent on changes in both ice flow speed and thickness, we would expect that, following front retreat, discharge would (i) initially increase with acceleration and then decrease due to dynamic thinning if the front restabilizes or (ii) become less sensitive to front position change if retreat continues up a prograde bed slope. While such a response is not clearly apparent for the GrIS-wide results (Fig. 2), regionally we find temporary declines in D following retreat. This is particularly evident in the SE and CE regions, where D has declined and stabilized during the recent period of relative front stability (Supplementary Fig. 5). In the NW, we find an increasing sensitivity of discharge to ubiquitous retreat, likely due to several of the actively retreating glaciers doing so on increasingly retrograde slopes, resulting in a greater increase in ice flux with retreat20,21,22,23.
The ice sheet-wide and regional sensitivities of discharge to ice front change are consistent on the scale of individual glaciers. Out of the 128 glaciers with continuous records between 1985 and 2018, nearly 70% display a significant (p ≤ 0.05) relationship between changes in discharge and front position. These glaciers account for ~75% of the total D. Only 24 (19%) glaciers, accounting for 6% of D, exhibit a significant negative correlation (i.e., decreasing D with retreat) between discharge and front position. These glaciers tend to have more static fronts, varying only ~1 km between 1985 and 2018, compared with 4 km for those where positive correlations are observed, and include quiescent phase surge-type glaciers and those that have retreated along prograde bed slopes24 (i.e., Upernavik Isstrøm) during the observation period.
Retreat-driven acceleration initiates at the ice front within a distance determined by the stress coupling length (e.g., ref. 13) and then propagates up-glacier through diffusion. Therefore, the timing between retreat and changes in D may be sensitive to the location of the flux gates with respect to the ice fronts. To test this sensitivity, we repeated the above analysis using a subset of the 15 largest glaciers with dynamic flux gates positioned 2 km from the changing ice front and obtained nearly identical results, indicating that our fixed flux gates are close enough to the ice front to capture nearly instantaneous changes in front-driven dynamics. More detail regarding this analysis and an example is provided in the Supplementary Methods (see also Supplementary Fig. 6).
Temporal patterns of thinning and retreat
Retreat occurs through increased calving due to melting at the ice front and/or thinning of the glacier terminus toward flotation. To assess which mechanism triggered the onset GrIS-wide retreat of the late 1990’s/early 2000’s, we investigate the timing and magnitudes of regional retreat and thinning, weighted by discharge, as above, and plotted in Fig. 3. Between 1985 and 2000 (i.e., between the AeroDEM and start of the ASTER surface elevation datasets), glaciers thinned across their sampled flux gates by an average of 2–4% in all regions except the central west, where a 6% thinning was primarily due to the initial acceleration of Jakobshavn Isbræ in 1998–1999. During most of this initial period of thinning, ice fronts slightly advanced in the west and retreated steadily in the east since the early 1990’s. The onset of rapid retreat around the year 2003/2004 was accompanied by a brief 3–4-fold increase in thinning rates in the east, abruptly ending when fronts temporarily stabilized. In contrast, thinning rates were more constant in the west, with a doubling in thinning in the NW and little change in the CW (Fig. 3) after the onset of rapid retreat around the year 2000. We note that in the NW, the disintegration of floating ice tongues at large glaciers, as observed at Kong Oscar Glacier25, likely resulted in a delayed dynamic increase of discharge, beginning around 2003, in response to retreat (Supplementary Fig. 5). More constant thinning rates along the western margin mirrored long-term trends in retreat, which slightly decelerated after 2005 in both the NW and CW. Regional plots similar to those shown in Fig. 3 using unweighted means are shown in Supplementary Fig. 7.
Thus, in all regions, substantial thinning occurred prior to the onset of rapid retreat and acceleration in the early 2000’s. To assess the potential contribution of increased surface melt to thinning, we examined changes in surface mass balance (SMB) estimated by RACMO2.3p2 relative to the 1961–1990 baseline average. We estimate thinning due to SMB anomalies along glacier flux gates, located in the ablation zone, by assuming a density of 910 kg m−3. Since thinning is the sum of declining SMB and increased dynamic thinning due to ice stretching, we can estimate their relative contributions to the observed thinning for the same 1961–1990 mean reference mass balance as used above. Discharge-weighted, regional averages of thickness change due to SMB are plotted in Supplementary Fig. 8, along with the total observed thickness change and thickness change due to ice dynamics, estimated as the difference between thickness change due to SMB and the total.
Thickness changes due to SMB between 1985 and 2000 were much smaller than the observed thickness changes for all regions except the SE, where an average thinning of 8 m due to declining SMB is very close to the total observed thinning. Conversely, ice dynamics accounted for nearly all of the rapid thinning, averaging ~20 m, during the episode of retreat and acceleration in the SE between 2002 and 2005. Since 2008, however, the rates of thinning due to SMB and dynamics have been nearly equal, each accounting for ~1 m yr−1.
Thinning due to declining SMB in the other regions began later than in the SE; the decline in SMB started in the mid-1990s in the CE, and the late 1990s in the NW and CW. In these regions, the contribution of SMB to total thinning has increased through time, reaching 30–35% in the CW and NW by 2018. Less thinning was due to SMB in the CE, accounting for <10% of the total thinning in that region by 2018.
