Significance We reconstruct the mass balance of the Greenland Ice Sheet for the past 46 years by comparing glacier ice discharge into the ocean with interior accumulation of snowfall from regional atmospheric climate models over 260 drainage basins. The mass balance started to deviate from its natural range of variability in the 1980s. The mass loss has increased sixfold since the 1980s. Greenland has raised sea level by 13.7 mm since 1972, half during the last 8 years.

Abstract We reconstruct the mass balance of the Greenland Ice Sheet using a comprehensive survey of thickness, surface elevation, velocity, and surface mass balance (SMB) of 260 glaciers from 1972 to 2018. We calculate mass discharge, D, into the ocean directly for 107 glaciers (85% of D) and indirectly for 110 glaciers (15%) using velocity-scaled reference fluxes. The decadal mass balance switched from a mass gain of +47 ± 21 Gt/y in 1972–1980 to a loss of 51 ± 17 Gt/y in 1980–1990. The mass loss increased from 41 ± 17 Gt/y in 1990–2000, to 187 ± 17 Gt/y in 2000–2010, to 286 ± 20 Gt/y in 2010–2018, or sixfold since the 1980s, or 80 ± 6 Gt/y per decade, on average. The acceleration in mass loss switched from positive in 2000–2010 to negative in 2010–2018 due to a series of cold summers, which illustrates the difficulty of extrapolating short records into longer-term trends. Cumulated since 1972, the largest contributions to global sea level rise are from northwest (4.4 ± 0.2 mm), southeast (3.0 ± 0.3 mm), and central west (2.0 ± 0.2 mm) Greenland, with a total 13.7 ± 1.1 mm for the ice sheet. The mass loss is controlled at 66 ± 8% by glacier dynamics (9.1 mm) and 34 ± 8% by SMB (4.6 mm). Even in years of high SMB, enhanced glacier discharge has remained sufficiently high above equilibrium to maintain an annual mass loss every year since 1998.

In the last several decades, the Greenland Ice Sheet (GIS) has lost mass to the ocean (1⇓⇓⇓–5). The mass loss has been quantified by three independent techniques using changes in ice volume (6, 7), time-variable gravity (8), and input versus output fluxes or mass budget method (1, 2, 4, 9⇓–11), for the time period 1992–2016 or 2002–2016. The mass budget method is the only one that provides information about the physical processes controlling the mass loss, i.e., the partitioning between surface mass balance (SMB) processes (accumulation minus runoff and other forms of ablation) and glacier dynamics (ice mass flux into the ocean), which is important to inform numerical models. A drawback of this method is that it requires comprehensive, precise glacier fluxes into the ocean and reconstruction of SMB over the ice sheet, i.e., the differencing of two large numbers. The gravity method does not extend before 2002. The ice volume method does not extend before 1992 with satellite data. Aerial photos were used to quantify ice volume changes in coastal areas from 1900 to 1980s using a single Digital Elevation Model (DEM) (5).

Here, we extend the mass budget method to the start of the Landsat historical archive in 1972, 20 years longer than with altimetry, and 30 years longer than with gravity. The revision benefits from a more complete time series of ice velocity (12⇓⇓⇓–16), improvements in ice thickness from NASA’s Operation IceBridge (OIB) (17, 18), bathymetric surveys from NASA’s Ocean Melting Greenland (OMG), and gravity surveys from OIB and the Gordon and Betty Moore Foundation (19⇓–21). A mass-conservation method constrained with a high-resolution velocity vector map yielded a high-resolution (350 m) ice thickness and bed elevation map of Greenland, named “BedMachine” (22, 23), based on physical principles instead of an interpolation (24). We use BedMachine Version 3 combined with new gravity inversion data in southeast Greenland. We benefit from major improvements in surface topography mapping. We use a 30-m-spacing Greenland Ice Mapping Project DEM for 2007–2008 (25), 8-m-spacing time series of WorldView DEMs (Polar Geospatial Center, University of Minnesota) for 2011–2018, and a 30-m-spacing historic DEM for 1980s (26). When combined with ice thickness and bed elevation at the date of the radar surveys, OIB and pre-OIB airborne laser altimetry from 1993–2017, the DEMs yield a time series of glacier thickness to calculate glacier fluxes with precision since 1972. Finally, regional atmospheric climate models used to reconstruct SMB have improved in spatial resolution (5.5 km downscaled to 1 km instead of 11 km) to match, in size, the typical width of outlet glaciers, which improves the fidelity of reconstruction of ice melt at low elevation (27, 28). We present the methodology; discuss the history of Greenland mass balance over the last 46 years until the most recent data, glacier by glacier, region by region, for the entire ice sheet; compare the estimates with prior work; and conclude on the recent and near-term contributions of the GIS to sea level rise.

Results We divide Greenland, including its peripheral glaciers and ice caps, into 260 basins (Fig. 1A and Dataset S1; ref. 29) grouped in seven regions: (i) southwest (SW), (ii) central west (CW), (iii) northwest (NW), (iv) north (NO), (v) northeast (NE), (vi) central east (CE), and (vii) southeast (SE). These regions are selected based on ice flow regimes (30), climate (31), and the need to partition the ice sheet into zones comparable in size (200,000 km2 to 400,000 k m 2 ) and ice production (50 Gt/y to 100 Gt/y, or billion tons per year). Out of the 260 surveyed glaciers, 217 are marine-terminating, i.e., calving into icebergs and melting in contact with ocean waters, and 43 are land-terminating, i.e., with zero discharge at the terminus (Dataset S2). The actual number of land-terminating glaciers is far larger than 43, but we lump them into larger units for simplification because we only need the total SMB to conduct the assessment. Peripheral glaciers and ice caps are assumed to be in balance at the beginning of our survey, and only SMB processes are considered afterward (32). We calculate the mass balance as SMB over the drainage basin minus ice discharge, D, at the glacier grounding line, or at the ice front if the glacier does not develop a floating section. Results are added by regions and for the entire ice sheet (Dataset S2). We calculate mass balance on a decadal time scale to reduce errors. For SMB, we use the output products from the recent Regional Atmospheric Climate Model v2.3p2 downscaled at 1 km (28). We use BedMachine to calculate the eustatic sea level equivalent (SLE) of each basin and region. Fig. 1. (A) Glacier catchments/basins for the GIS and seven regions overlaid on a composite map of ice speed (12). (B–D) For 1972–2018, the percentage (B) thickness change, (C) acceleration in ice flux from each basin, and (D) cumulative loss per basin. The surface area of each circle is proportional to the change in ice discharge caused by a change in (B) thickness or (C) speed; the (blue/red) color indicates the (positive/negative) sign of the change in (B) thickness, (C) speed, and (D) mass. SW holds a 74-cm SLE over an area of 216,207 k m 2 controlled at 72% by land-terminating glaciers. Twelve tidewater glaciers discharge, on average, 30 ± 5 Gt/y in 1972–2018. Changes in D are controlled by Qajuutap Sermia (4.5 ± 1.1 Gt/y in 1987), Ukaasorsuaq (6.4 ± 1.9 Gt/y), and Kangiata Nunaata Sermia (6.4 ± 1.9 Gt/y). We find no trend in D during 1972–2018. In contrast, SMB decreased from 59 ± 4 Gt/y in 1972–1980, 40 ± 4 Gt/y in 1990–2000 (SI Appendix, Fig. S1A and Dataset S2), and −13 ± 2 Gt/y in 2010–2018. The total mass balance decreased from +29 ± 6 Gt/y in the 1970s to −43 ± 4 Gt/y in 2010–2018 (Fig. 2). Overall, SW gained 426 ± 80 Gt between 1972 and 2001 (Fig. 3A) and lost 486 ± 17 Gt between 2001 and 2018, for a net loss of 63 ± 98 Gt or 0.2 ± 0.3 mm eustatic SLR, the smallest contributor in Greenland. Fig. 2. Partitioning of the mass loss between anomalies in SMB, dSMB, and ice discharge, dD, by regions for the time period 1972–2018 in gigatons, or 10 12 kg; (A) 1972–1980, (B) 1980–1990, (C) 1990–2000, (D) 2000–2010, and (E) 2010–2018. Shown are SMB in light tone (red for loss, blue for gain), and D in dark tone (red for loss, blue for gain). The size of the circle is proportional to the absolute magnitude of the change in SMB or D. Fig. 3. Cumulative anomalies in SMB (blue), discharge (D, red), and mass (M, purple) in gigatons (gigaton = 10 12 kg) for the time period 1972–2018 for the seven regions of Greenland and the entire ice sheet component: (A) SW, (B) CW, (C) NW, (D) NO, (E) NE, (F) DE, (G) SE, and (H) GIS. CW holds a 134-cm SLE over an area of 236,648 k m 2 drained at 91% by 15 tidewater glaciers (SI Appendix, Fig. S1B). Jakobshavn Isbræ controls 45% of D (30.0 ± 7.7 Gt/y in 1975), followed by 32% from Store Gletscher (8.6 ± 1.4 Gt/y in 1975) and Rink Isbræ (11.0 ± 1.4 Gt/y in 1984). D decreased from 68 ± 4 Gt/y in 1972–1980 to 65 ± 1 Gt/y in 1990–2000 due to a slowdown of Jakobshavn Isbræ (33), jumped to 93 ± 8 Gt/y in 2013 due to the speed up of Jakobshavn Isbræ, and decreased to 78 ± 6 Gt/y in 2018. Three other glaciers, Eqip Sermia, Kangilerngata Sermia, and Sermeq Silarleq, increased their ice flux by almost 100% from 1984–1998 to 2018. Store Gletscher and Rink Isbræ fluctuated little. SMB dropped from 62 Gt/y in 1972–1989 to 38 ± 3 Gt/y in 2010–2018, or 39%. SMB and D were in balance between 1972 and 2000, but combined to form a large loss (679 ± 15 Gt) from 2000 to 2018 (Fig. 2). In total, CW lost 738 ± 75 Gt, or 2.0 ± 0.2 mm SLR since 1972 (Fig. 3B). NW holds a 127-cm SLE over 283,654 k m 2 drained by 64 tidewater glaciers (Fig. 1A). D decreased from 90 ± 3 Gt/y in 1972–1980 to 87 ± 1 Gt/y in 1990—2000, and increased to 112 ± 1 Gt/y by 2010–2018, or 45% (SI Appendix, Fig. S1C). The largest changes occurred for Upernavik Isstrøm C (+73 ± 16% for 1993–2018), Upernavik Isstrøm N (+141 ± 12% for 1996–2013), Kakivfaat Sermiat (+136 ± 11% for 1995–2015), Alison (+169 ± 10% for 1995–2018), Kjer (+372 ± 30% for 2005–2018), Steenstrup-Dietrichson (+93 ± 13% for 1985–2012), and Sverdrup (+112 ± 22% for 1986–2018). A number of glaciers with no speedup were already out of balance in the 1970s: Hayes Gletscher M and SS, and Upernavik S. SMB averaged 78 ± 2 Gt/y for 1972–1980 but dropped to 45 ± 2 Gt/y in 2010–2018. NW changed from near balance in the 1970s to a small loss in the 1980s, then equilibrium between 1995 and 2000, before a rapid loss from 2000 to present. The cumulative loss is 1,578 ± 56 Gt, or 4.4 ± 0.2 mm SLR, the largest contributor in Greenland (Fig. 3C). NO holds a 93-cm SLE over 263,534 k m 2 drained at 82% by 12 tidewater glaciers (Fig. 1A). D varied from 23 ± 1 Gt/y in 1972–1980, to 22 ± 1 Gt/y in 1990–2000, to 24 ± 1 Gt/y in 2010–2018. The largest discharge is from Humboldt (4.8 ± 1.1 Gt/y in 1975), Petermann (10.5 ± 1.4 Gt/y in 1975), Ryder (3.0 ± 0.2 Gt/y in 1985), and Ostenfeld (1.7 ± 0.2 Gt/y in 1984). Humboldt discharge increased from 4.8 ± 1.1 Gt/y in 1972 to 6.4 ± 1 Gt/y in 2018. Petermann remained relatively stable until 2010, when it experienced major calving events (34) and increased D by 10%. SMB decreased from 25.4 ± 1.5 Gt/y in 1972–1980 to −5 ± 2 Gt/y in 2010–2018. The cumulative mass loss is 474 ± 30 Gt, or 1.3 ± 0.1 mm SLR, the largest loss per unit discharge. NE holds a 180-cm SLE over 425,250 k m 2 controlled at 87% by 14 tidewater basins, mostly Nioghalfjerfjorden (60-cm SLE), Zachariae Isstrøm (56-cm SLE), and Storstrømmen (26-cm SLE) (Fig. 1A). D for Nioghalfjerfjorden has increased by 10% over the entire period. Zachariæ Isstrøm sped up after 2012 (35, 36), doubling D from 10 ± 1 Gt/y in 1972 to 18 ± 2 Gt/y in 2018. Storstrømmen is a surge-type glacier with 1.2 ± 1.0 Gt/y discharge in 1972, 20.6 ± 3.4 Gt/y in 1979, and 0 Gt/y after 1990 as the next surge builds up (37, 38). Overall, D increased from 36.2 ± 1 Gt/y in 1972–1980, to 40 ± 1 Gt/y in 1980–1990, back to 31 ± 1 Gt/y in 1990–2000, and up to 35 ± 1 Gt/y in 2010–2018. SMB dropped from 22 ± 1.3 Gt/y in 1961–1990 to 13 ± 2 Gt/y in 2010–2018. The cumulative loss is 532 ± 52 Gt since 1972, or 1.5 ± 0.1 mm SLR. CE holds a 72-cm SLE over 218,628 k m 2 drained by 41 tidewater glaciers (Fig. 1A). The largest D is from Kangerlussuaq (21.1 ± 2.0 Gt/y in 1984), Unnamed Deception ø CN and CS (8.2 ± 0.7 Gt/y), and Daugaard-Jensen (8.9 ± 0.7 Gt/y in 1980s). D increased moderately from 75 ± 5 Gt/y in 1972–1980 to 78 ± 2 Gt/y in 1990–2000, and 87 ± 2 Gt/y in 2010–2018. In single years, D increased 20 ± 9 Gt/y in 2005, 10 ± 7 Gt/y in 1991, and 7 ± 8 Gt/y in 1996 due to the speedup of one to two glaciers. SMB decreased from 73 ± 1.5 Gt/y in 1960–1989 to 60 ± 3 Gt/y for 2010–2018. The mass loss since 1972 totals 508 ± 83 Gt, or 1.4 ± 0.3 mm SLR. SE holds a 55-cm SLE over 165,349 k m 2 drained by 59 tidewater glaciers. Its mass budget has been difficult to assess due to uncertainties in glacier thickness and speed in an area of high rates of snowfall with liquid water in the firn (39). Gravimetric surveys (21), combined with OMG bathymetry in the fjords, solved the problem of calculating D at the ice front. SE has the largest D (136 ± 6 Gt/y in 1972–1980) in Greenland. D increased to 158 ± 2 Gt/y in 2000–2010 and 160 ± 2 Gt/y in 2010–2018. The rapid increase in 2000–2004 was driven by Helheimgletscher and less well-known glaciers Køge Bugt C and S, Umiivik Fjord, A. P. Bernstoff, Tingmiarmiut Fjord, and Anorituup Kangerlua. SMB decreased from 131 ± 1.6 Gt/y in 1960–1989 to 111 ± 3 Gt/y in 2010–2018. SE was near balance before the 1990s but has lost 1,089 ± 91 Gt since 1972, or 3.0 ± 0.3 mm SLR. Greenland’s SMB averaged 422 ± 10 Gt/y in 1961–1989 (SI Appendix, Fig. S1H). It decreased from 506 ± 18 Gt/y in the 1970s to 410 ± 17 Gt/y in the 1980s and 1990s, 251 ± 20 Gt/y in 2010–2018, and a minimum at 145 ± 55 Gt/y in 2012. In 2018, SMB was above equilibrium at 449 ± 55 Gt, but the ice sheet still lost 105 ± 55 Gt, because D is well above equilibrium and 15 Gt higher than in 2017. In 1972–2000, D averaged 456 ± 1 Gt/y, near balance, to peak at 555 ± 12 Gt/y in 2018. In total, the mass loss increased to 286 ± 20 Gt/y in 2010–2018 due to an 18 ± 1% increase in D and a 48 ± 9% decrease in SMB. The ice sheet gained 47 ± 21 Gt/y in 1972–1980, and lost 50 ± 17 Gt/y in the 1980s, 41 ± 17 Gt/y in the 1990s, 187 ± 17 Gt/y in the 2000s, and 286 ± 20 Gt/y in 2010–2018 (Fig. 2). Since 1972, the ice sheet lost 4,976 ± 400 Gt, or 13.7 ± 1.1 mm SLR.

Discussion Our assessment extends prior records in time by 20 to 30 years and in quality by more than 20%. We have fewer velocity data in 1972–1992 compared with 1992–2018, but the glacier variations in speed are also less in 1972–1992 (11, 18, 40). Our long time record makes it possible to examine decadal estimates of mass balance instead of yearly estimates, which reduces our errors by a factor of 3. We constrain 85% of D with precision thickness and 15% from velocity-scaled reference fluxes (see Materials and Methods). For precision thickness, we use BedMachine at the ice front of glaciers contributing to 47% of D, contrary to ref. 11, who use it for all glaciers. For the glaciers that we exclude, the uncertainty in BedMachine exceeds our requirements (±100 m) and yields errors in flux up to 100%. For the remaining 38% of D, we use gravity-derived thickness (13%) and direct radar-derived thickness upstream of ice fronts (25%). The gravity-derived thickness used for SE (21) changes its balance flux from 47.2 Gt/y to 64.2 Gt/y, or +37%, thereby enhancing the role of SE in the total budget. For the glaciers using a velocity-scaled reference flux, a 10% uncertainty in reference flux yields a 1.5% error in total mass balance, or 4 Gt/y, which is negligible. Our study also uses systematic corrections for ice thickness. The uncertainty in ice thickness correction is less than 1%. Without it, D would be 10% too high in 2018, or 55 Gt/y. Radar altimetry indicates a loss of 269 ± 51 Gt/y in 2011–2014 (7), and laser altimetry indicates a loss 243 ± 18 Gt/y for 2003–2009 (6). We find 323 ± 28 Gt/y and 220 ± 21 Gt/y, respectively, for the same periods; that is, our estimates agree within errors, especially laser altimetry. The loss is lower with radar altimetry because of the difficulty for radar altimeters to sample coastal sectors with steep slopes. With the Gravity Recovery and Climate Experiment (GRACE), the loss of 280 ± 58 Gt/y for 2003–2013, with an acceleration of 254 ± 12 Gt/y per decade (8), matches our 274 ± 17 Gt/y estimate, with an acceleration of 249 ± 54 Gt/y per decade. For 1992–2011, the 142 ± 49 Gt/y mass loss from multiple methods (41) is within errors of our 144 ± 13 Gt/y estimate. Few estimates exist before the 2000s (42⇓⇓–45), because of mission duration (GRACE launched in 2002) or limitations over steep slopes (radar altimeters). Our assessment indicates with certainty that the ice sheet was near balance (−13 ± 14 Gt/y) in 1972–1990. Within that time period, SMB was above balance in 1972–1980 and slightly below balance in 1980–1990, whereas D was slightly above equilibrium the entire time. This result contradicts an earlier study that concluded the ice sheet was losing mass the entire time (5) based on differences in ice elevation between two epochs 83 years apart, with no data in the ice sheet interior. We posit that balance conditions in 1972–1990 were compensated by periods of high loss, e.g., in the 1920s through 1940s (46), so that the mass loss varied significantly between the two epochs. Secondly, the same study suggested that the contribution from glacier dynamics to the mass loss was constant over the entire period 1900–2016. Our annual to subannual time series show that D increased continuously during 1972–2018 and overall has had a larger impact on mass balance than SMB. SMB domination of the mass balance is only recent. Over the past 46 years, ice dynamics contributed 66 ± 8% of the mass loss versus 34 ± 8% for SMB. Hence, ice dynamics, i.e., changes in glacier flow, play a major role in the mass budget. In 2000–2012, half of the cumulative loss was from four glaciers, (i) Jakobshavn Isbræ, (ii) Kangerlussuaq, (iii) Køge Bugt, and (iv) Ikertivaq S (4), but, between 1972 and 2018, Ikertivaq S gained 26 ± 15 Gt. Over 46 years, the conclusions are different. The largest losses are from (i) Jakobshavn Isbræ (327 ± 40 Gt), (ii) Steenstrup-Dietrichson in NW (219 ± 11 Gt), (iii) Kangerlussuaq in CE (158 ± 51 Gt), (iv) Humboldt in NO (152 ± 7 Gt), (v) Midgårdgletscher in SE (138 ± 5 Gt), and (vi) Køge Bugt C in SE (119 ± 37 Gt), hence highlighting glaciers that are seldom mentioned in the literature. Steenstrup-Dietrichson, Humboldt, and Midgårdgletscher contributed to the mass loss during the entire period, versus only after year 2000 for Jakobshavn, Kangerlussuaq, and Helheimgletscher. This result illustrates the risk of summarizing the ice sheet loss on the basis of the fate of a few glaciers. Some glaciers gained mass during the survey: Saqqap, Majorqaq, and Russell glaciers in SW gained 282 ± 18 Gt in 1972–2006. This mass gain is consistent with the glacier advance in SW in the 1970s to 1980s (47, 48, 49) and satellite-derived ice sheet growth (50, 51), which improves confidence in the SMB reconstruction. Conversely, a few large glaciers did not undergo large dynamic changes (<10%): Rink Isbræ, Hayes N and NN, Petermann, Nioghalfjerdfjorden, and Daugaard-Jensen glaciers. We attribute their stability to the bed configuration. Daugaard-Jensen calves on a stabilizing ridge. Nioghalfjerdfjorden retreats along a prograde bed. Rink and Petermann are protected by a buttressing ice shelf. Petermann and Nioghalfjerdfjorden sped 10% since 2010 and 2006, respectively, following a weakening of their buttressing ice shelves (34, 36). In terms of partitioning between discharge and SMB (Fig. 2), the loss from D anomalies (deviations in D from a reference state) decreased from 47 ± 19 Gt/y in the 1970s to 41 ± 8 Gt/y in the 1990s before reaching 127 ± 9 Gt/y in 2010–2018, for a cumulative 3,312 ± 124 Gt since 1972, or 9.1 ± 0.3 mm SLE (Fig. 3H). SMB anomalies (deviations in SMB from a reference state) were positive in 1970s (+95 ± 20 Gt/y) and near equilibrium in the 1980s and 1990s (−1 ± 17 Gt/y and 0 ± 17 Gt/y, respectively) before becoming negative in 2000–2018 (−99 ± 17 Gt/y and −160 ± 20 Gt/y for 2000–2010 and 2010–2018, respectively), for a cumulative loss of 1,670 ± 379 Gt since 1972, or 4.6 ± 1.0 mm SLE. Hence, glacier dynamics has played a stronger role in the mass loss (66 ± 8%) than SMB (34 ± 10%) during the last 46 years. SMB dominated (55 ± 5%) only in the last two decades. It is important to note, however, that D increased by 18% during that time, versus a 38% decrease in SMB; that is, D had a larger impact on the mass loss, because it was already above equilibrium conditions in the 1970s. This result has several implications. First, even in years of high SMB, e.g., 2018, the ice sheet loses mass because D is well above equilibrium. In fact, the ice sheet has lost mass every year since 1998. Second, this implies that changes in glacier flow remain of primary importance in driving the mass loss of the ice sheet, even though SMB played a larger role in the last 20 years. In the future, SW will remain controlled by SMB processes, whereas SE, CW, and NW will be controlled by the fate of their tidewater glaciers, i.e., ice dynamics (64 ± 12%, 65 ± 14%, and 86 ± 4%, respectively). We expect that NW, SE, and CW will continue to dominate the ice sheet mass budget, especially in NW where D increased 18 ± 0.1 Gt/y per decade from 1998 to 2018, with no sign of deceleration. In terms of long-term contribution to sea level rise, the NO and NE sectors are of greatest importance. The loss in that region is equally partitioned between D and SMB at present (62 ± 11% and 60 ± 15%, respectively). The glaciers do not produce a high D (25.9 ± 1.9 Gt/y and 39.5 ± 2.7 Gt/y, respectively, in 2018), but the glacier speeds are low compared with those in SE or NW; the potential for a large increase in D exists if the glaciers lose their buttressing ice shelves and start flowing as fast as their SE and NW counterparts, which would tap the largest potential SLE (273 cm) in Greenland. The evolution of NO and NE glaciers in the coming decades is therefore of greatest relevance to future sea level change as the ice shelves are weakened by climate change.

Conclusions Using improved records of ice thickness, surface elevation, ice velocity, and SMB, we present a 46-years reconstruction of glacier changes across Greenland that reveals a dominance from ice discharge over the entire record and a sixfold increase in mass loss from the 1980s. The largest mass loss is from NW, SE, and CW, which are controlled by tidewater glaciers. Several glaciers, including Humboldt, Steenstrup-Dietrichson, and Køge Bugt C. North Greenland (NO, NE) have played a stronger role in the total mass loss than reported previously, hence illustrating the value of an extensive time series of mass balance that includes all of the large glaciers. We also find that the ice sheet as a whole was near balance over the time period 1972–1990. In the future, we expect the mass changes in the northern part of Greenland to become of greatest importance to sea level rise, because of the large reserve of ice above sea level and the potential for manyfold increase in ice discharge.

Materials and Methods Ice Discharge. We combine ice thickness (19⇓⇓–22) and ice velocity (13⇓⇓–16) to calculate ice flux, D (SI Appendix). We assume no internal deformation; that is, surface ice speed equals the depth-averaged ice speed. We apply no firn correction. At the flux gate, D in cubic kilometers per year is the integral (trapezoidal rule) of ice velocity v i in the direction normal to the gate times ice thickness H i . Speed and ice thickness are sampled every 200 m. The associated error, σ D , is the sum of systematic and random errors as σ D = ∑ i v i σ H i + ∑ ( H i σ v i ) 2 , where σ H i is the error in ice thickness and σ v i is the error in speed. The volume flux is converted into mass using an ice density of 917.2 kg/ m 3 . Flux gates are placed at the most-retreated grounding line or ice front position. With gridded data sets from mass conservation (MC) and gravimetric inversion (GRA), we calculate the flux within 1 km of the most-retreated grounding line or ice front position. With direct radar sounder measurements (GAT), the measured flux at the flux gate is converted to a reference ice front flux, D r e f , as in ref. 2. Namely, we correct the flux obtained at the flux gate with the earliest velocity data by adding the average SMB for 1961–1989 (or S M B r e f ) between the flux gate and the ice front; that is, we assume that the discharge was in equilibrium at that time. A 2% error in that initial value (D increased by 18% in 1972–2018 on average) would introduce a 0.7% error in D r e f . We use 61, 16, and 20 flux gates for MC, GRA, and GAT, respectively. These 107 glaciers control 85% of D. For the other glaciers, we use D r e f equal to S M B r e f . The error in D combines the error in S M B r e f and error in the relationship between mean SMB and initial D. The linear regression between the GAT, MC, and GRA values of D before the 1990s and S M B r e f reveals a significant correlation ( R 2 = 0.93; SI Appendix, Fig. S2). Using a linear least-squares fit between these variables, the unexplained variance between mean SMB and D is ±7%. Detailed results for each glacier, gate, and area used for SMB correction are in SI Appendix. Annual D is calculated by scaling D r e f with the speed and thickness changes established previously. The associated error, σ D , is the combination of the systematic error in reference flux and independent errors on the scaling as σ D = α v α h σ D r e f + ( α v D r e f σ α h ) 2 + ( α h D r e f σ α v ) 2 , where α v is the scale factor for speed, α h is the scale factor for thickness, and D r e f is the reference flux. For years with no velocity data, D is linearly interpolated, or kept constant if at the end or beginning of the time series. If D is extended at the end or beginning of the time series, we assume that the relative error increases by 5% per year beyond the closest available estimate, or 100% after 14 years. When D is interpolated, we assume a relative error increasing at 2.5% per year from the closest available estimate. We assume that errors in D are independent between individual glaciers, so the error in D is ∑ σ D 2 . Mass Loss. Mass balance, dM/dt, is SMB minus ice discharge D. The error in dM/dt combines the errors in SMB and in D (SI Appendix, Fig. S2 and Dataset S2). We estimate dM/dt yearly between 1972 and 2018 (Fig. 2). For each basin, we compute the partitioning between anomalies in D (D – S M B r e f ) and anomalies in SMB ( S M B r e f – SMB). The anomaly in D is spread spatially using the ratio between flux density (speed times thickness) and distance to the ice front as a proxy. The anomaly in SMB is naturally spread over the ice sheet. Mass (in gigatons) is converted to eustatic sea level rise using 362 Gt = 1 mm SLE.

Acknowledgments This work was performed at the University of California, Irvine, and at California Institute of Technology’s Jet Propulsion Laboratory under Contract NNX13AI84 A with the NASA Cryosphere Science Program. J.M. acknowledges funding from the French Centre National d’Etudes Spatiales (CNES) and French L’Agence Nationale de la Recherche (ANR) (Grant ANR-15-IDEX-02). M.v.d.B. and B.N. acknowledge funding from the Polar Program of the Netherlands Organization for Scientific Research and the Netherlands Earth System Science Center. A.A.B. acknowledges support from the Carlsberg Foundation (Grant CF17-0529). We acknowledge the European Space Agency (ESA), Canadian Space Agency (CSA), Japanese Space Agency (JAXA), Italian Space Agency (ASI), and German Space Agency (DLR) for the use of SAR data, NASA and US Geological Survey for the Landsat data, NASA’s OIB and IceSAT missions for the use of ice thickness and elevation data. We thank the Space Task Group and Polar Space Task Group for coordinating satellite data acquisition efforts during IPY and post-IPY, respectively. Worldview DEMs were provided by the Polar Geospatial Center under National Science Foundation OPP Awards 1043681, 1559691, and 1542736.

Footnotes Author contributions: J.M. and E.R. designed research; J.M. performed research; J.M., A.A.B., M.v.d.B., R.M., M.M., B.N., B.S., and M.W. analyzed data; and J.M. and E.R. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Data deposition: Ice velocity data and basin boundaries used to perform the analysis have been deposited in University of California, Irvine Dash (UCI Dash) (DOIs: 10.7280/D1MM37, 10.7280/D1GW91, 10.7280/D1595V, 10.7280/D1595V, and 10.7280/D11H3X), the ice thickness data are available at National Snow and Ice Data Center (NSIDC) and UCI Dash (DOIs: 10.5067/2CIX82HUV88Y, 10.5067/GDQ0CUCVTE2Q, and 10.7280/D1NW9K), and the surface elevation are available at NSIDC (DOIs: 10.5067/ICESAT/GLAS/DATA209 and 10.5067/CPRXXK3F39RV) and the Polar Geospatial Center (https://www.pgc.umn.edu/data/arcticdem/).

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1904242116/-/DCSupplemental.