Introduction

The Greenland Ice Sheet (GrIS) has been losing mass for several decades1 due to both increased surface meltwater runoff and ablation of marine-terminating outlet glaciers via calving and submarine melting, termed ice discharge. Total GrIS mass loss over the 1992–2018 period was due to approximately equal contributions from both terms1 but with greater contribution from increased melt runoff after 2000, when mass losses accelerated2,3,4. Estimates of ice sheet discharge over multiple decades and at annual, or finer, resolution provide insight into the ice sheet’s response to long-term climate forcing and ongoing change4,5. Seasonal and interannual variability in ice sheet dynamics are challenging to resolve prior to the year 2000 due to temporal and spatial data gaps. Here, we apply the rigorous methodology in ref. 3 to gain improved constraints on estimates of ice sheet discharge over three decades, including the period leading up to the onset of rapid glacier retreat and acceleration. Rates of Greenland glacier retreat have accelerated6 and previous work has identified relationships between glacier speed and retreat7 and glacier area8 for smaller subsets of Greenland glaciers. We also then combine these data with high-resolution observations of time-varying calving front position changes and perform a GrIS-wide analysis of how these two variables relate on individual, regional, and ice sheet-wide spatial scales over the multi-decadal record. We discuss the timing of changes in retreat, thinning, and acceleration across the ice sheet, quantify the sensitivity of ice discharge to retreat, and describe the roles of both long-term changes in ice dynamics and surface mass balance in preconditioning regions of the ice sheet for rapid retreat, thinning, and accelerated discharge.

We find that GrIS-wide discharge is now ~14% greater than the rate observed during 1985–1999, following an observed step-increase during the early 2000’s. Widespread glacier retreat explains nearly all (>90%) of the observed multi-decadal variability in discharge, with a observed increase in discharge of 4–5% per every weighted mean kilometer of retreat. We find that this sensitivity is proportionally consistent across different regions of the ice sheet, despite highly variable long-term trends in discharge. Following the step-increase in discharge, GrIS-wide totals have remained relatively stable at rates near 495–500 Gt yr−1, reflecting an increase that was sufficient to effectively shift the ice sheet to a state of persistent mass loss.

Results

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.

Fig. 1: Net ice sheet discharge at contributions to total mass balance.
figure 1

a Total GrIS discharge (black curve) with early (blue) and late-thinning (red dashes) scenarios between 1985 and 1999, and discharge estimates derived from wintertime InSAR velocities from Mouginot et al.4 in green. Gray shading is 2-σ uncertainty. b Net annual mass balance (Ba, black curve), annual surface mass balance (SMBa, blue dashes with circles) and total annual discharge (Da, red dashes with crosses). c Total shading shows the cumulative mass anomaly (Bc) of Ba since 1985 broken into cumulative contributions from discharge (Dc, red) and surface mass balance (SMBc, blue). Bc including SMB from peripheral ice caps, to match the GRACE domain, is plotted by dark blue line. GRACE estimates, updated from Wouters et al.16 for the 2002–2018 time period, are also plotted on this axis as a moving 12-month average (dotted black), and shifted along the y-axis to begin at the 2002 Bc value.

Full size image

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.

Fig. 2: GrIS-wide relationships between discharge and front position.
figure 2

a Annual total D versus net front position change for the 128 glaciers with continuous measurements between 1985 and 2018, with the color scale and label corresponding to the measurement year. The black line is the least-squares best fit with a slope of 14 Gt yr−1 per km of retreat and intercept of 312 Gt yr−1. Negative front change indicates net retreat from the glacier fronts’ initial position over the calendar year. Note the reversed x-axis to highlight the inverse relationship. b Annual changes in D and front position with labels indicating the first year (e.g., ‘00 is the change from 1999 to 2000). The (black line) least-squares fit has a slope of 14 Gt yr−1 per km of annual retreat and intercept of −0.7 Gt yr−1.

Full size image

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.

Fig. 3: Mean regional trends in ice thickness and front position.
figure 3

Net regional trends in (black curves) ice thickness and (colored curves) front position, calculated from discharge-weighted individual glacier time series in the northwest (a), central east (b), central west (c), and southeast (d), with the colors corresponding to the locations mapped in the left plane. Negative values on the left y-axis indicate net cumulative retreat over the study period. Cumulative thickness changes are shown on the right y-axes, and are calculated from the mean of all regional glaciers and expressed as fractional changes from Ho. Lighter, dotted gray lines between the two starred points indicate linear interpolation between two sparse data points, and the apparent trend during this period should be treated with caution. Please note that the limits along the left y-axes vary between panels so that time series details are more visible.

Full size image

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.

Discussion

We estimate over three decades of GrIS-wide ice discharge at high temporal resolution and with well-defined uncertainties, capturing a rapid, step-increase between two distinct states of discharge, and reaching an annual maximum of ~500 Gt yr−1 in 2017 and 2018. This transition corresponds to widespread retreat of glacier calving fronts at the turn of the 21st century after a long and relatively stable period of static ice fronts on the west coast and slow retreat in the east, possibly triggered by substantial ice thinning between 1985 and 2000. In all regions except the SE, outlet glacier thinning was due, at least in part, to glacier discharge exceeding the long-term balance flux, indicating dynamic disequilibrium associated with a longer-term retreat. This is consistent with the retreat of these glaciers earlier in the 20th century26. In the SE region, thinning prior to 2000 can be primarily attributed to increased surface melt, indicating a shorter-term response to climate forcing. This agrees with the observation that the largest SE glacier, Helheim, was close to mass balance before and after its rapid retreat and temporary acceleration27. As a result of the new, semi-static rate of GrIS-wide D, that is >60 Gt yr−1 higher than in 1985, annual SMB greater than two standard deviations above the 19-year (2000–2018) mean is required for the ice sheet to gain mass and is thus likely in a long-term state of persistent loss.

The overall increasing trend in ice discharge is the combination of substantially different regional, and even individual glacier, behavior that is likely the result of differences in ocean and atmospheric forcing17,28,29. Within the past 3 years, this variability ranges from declining D in the CW due to slowing of Jakobshavn Isbræ, to variable but trendless D in the SE, to episodically increasing D in the CE due largely to Kangerdlugssuaq, to continued acceleration in the NW. Thus, rather than indicating stabilization, recent changes in GrIS discharge are the result of partially offsetting regional responses to past and current forcing. Despite this regional variability, changes in discharge and front position show a remarkably consistent relationship, with every region showing a 4–5% increase in mean discharge per km of weighted mean retreat. This indicates relative uniformity in the processes that regulate outlet glacier response to changes at the calving front and provides some constraint on future mass loss. The step-increase in discharge in the early 2000’s was driven primarily by synchronous, rapid retreat in the SE and secondarily in the CE, where further regional-scale rapid retreat is unlikely due to shallowing glacier beds. The potential for another such episode of regional, rapid retreat is in the NW where retreat and discharge continue to accelerate. Finally, we expect continued episodic retreat of the large glaciers such as Jakobshavn Isbræ and Kangerdlugssuaq, that currently dominate net GrIS discharge since they are positioned on retrograde slopes22,23,29. Ultimately, predictions of future change will require improved understanding of the ice/ocean boundary and controls on glacier calving.

Methods

We have extended the ice discharge estimates of King et al.3 by two decades by integrating historical velocity data (LANDSAT 4 and 5), derived using the MIMC2 feature tracking algorithm30, with the AeroDEM digital elevation model31. We also append the extended time series, now covering an ~34-year period from 1985 to 2018, with annually resolved SAR (Synthetic Aperture Radar) velocity mosaics from ref. 4 between 1990 and 1999, which were computed mostly from wintertime observations. We use this series to assess decadal trends in ice discharge from the GrIS, and where possible, resolve the seasonal amplitude through time. We combine the dynamic ice loss component with 1-km downscaled estimates of surface mass balance (SMB) from RACMO2.3p2 over the full 1985–2018 time period15. We sum the two terms to derive more than three decades [1985–2018] of total ice mass balance that is independent of the total ice sheet mass balance derived from elevation change and gravimetry methods.

We derive ice thicknesses every 250 m along the flux gates by differencing time-varying surface elevation data, beginning with the earliest AeroDEM data through ArcticDEM32, generated by SETSM (Surface Extraction with TIN-based Search-space Minimization) algorithm33 with bed topography from BedMachine v334. We replace BedMachine v3 data with gravity inversion solutions from ref. 29 at select regions along the southeast margin, which improved bathymetry constraints for ten fjord systems using multibeam echo sounding data combined with airborne gravity data from Operation IceBridge. We combine spatially variable errors in bed topography, treated as systematic errors derived from the BedMachine v3 product, with an estimated random error of ±5 m from surface elevation data when computing total ice thickness errors. Following ref. 5, we apply a correction based on typical velocity-thickness relationships to a subset of 39 glaciers with ice thicknesses that are too small for flowing ice, located mostly within the southeast and central east regions of the ice sheet. Using the remaining thickness data, we derive unique solutions from the best linear fit between centerline velocity and flux-gate thickness. We find that thickness increases more gradually with increasing velocity at wider flux gates. Therefore, we discretize glaciers into subgroups based on glacier width along the flux gate (Supplementary Fig. 9) rather than using a single linear fit for all glaciers and apply an adjusted ice thickness term based on the discretized solutions. This binning method improves our baseline thickness adjustments at the individual glacier scale and increases net ice sheet discharge by a time series average of 13 ± 1 Gt yr−1.

We extend the 2000–2018 velocity record back in time to 1985 using LANDSAT 4 and 5, obtaining measurements for 133 of the 234 gates for the full record. These glaciers accounted for 78% of the net GrIS discharge between 2000 and 2018. Rather than assume constant flow and extrapolate the earliest available velocity data backward in time to account for incomplete coverage4, we account for the post-2000 acceleration of many glaciers by scaling the observed total annual D by the proportion of the mean net D accounted for by glaciers with observations before 2000. The mean contributions of individual glaciers to total D are calculated from their 2000–2005 period values, although results are closely similar if a 2000–2018 average is used instead. We fill missing monthly velocities using the procedure of ref. 3, where a Kalman filter is applied to each glacier time series using the glaciers’ typical seasonal variation. Missing monthly data are more prevalent in wintertime months prior to when TerraSAR/TanDEM-X InSAR data first became available in 2007. We do not apply the Kalman filter for gap filling prior to 1998 due to multi-year gaps. Therefore, the temporal resolution of estimates in D varies from approximately monthly in the summers of 1985–1997, when optical data are available (including 1985–1990 and 1994/1995), to continuous monthly estimates after 1998.

Surface elevation observations, required for estimating ice thickness, are limited before 2000. The nominal dates of AeroDEM data vary spatially between 1978 and 1987, with a nominal date of 198331, and therefore, temporal gaps range from 15 to 22 years between AeroDEM and available Advanced Spaceborne Thermal Emission and Reflection radiometer (ASTER) DEM data. 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. This range in the timing of thinning yields an average uncertainty in D of 13 Gt yr−1 during the 1985–1999 period (Fig. 1a).

We also combine velocity and ice thickness data with an extended time series of changes in calving front position at all 234 glaciers. Centerline front position data are manually digitized using optical imagery from ASTER and LANDSAT 4–8 following methods described in ref. 35, and filtered and weighted by glacier flux, when applicable, as described in ref. 3. Last, we compute net annual surface mass balance (SMB) and reference period GrIS-wide SMB using RACMO2.3p215, and combine this term with D to derive net ice sheet mass balance for 1985–2018. Front position and discharge data from this study are publically available36.