Greenland ice sheet climate disequilibrium and committed sea-level rise

Ice loss from the Greenland ice sheet is one of the largest sources of contemporary sea-level rise (SLR). While process-based models place timescales on Greenland’s deglaciation, their confidence is obscured by model shortcomings including imprecise atmospheric and oceanic couplings. Here, we present a complementary approach resolving ice sheet disequilibrium with climate constrained by satellite-derived bare-ice extent, tidewater sector ice flow discharge and surface mass balance data. We find that Greenland ice imbalance with the recent (2000–2019) climate commits at least 274 ± 68 mm SLR from 59 ± 15 × 103 km2 ice retreat, equivalent to 3.3 ± 0.9% volume loss, regardless of twenty-first-century climate pathways. This is a result of increasing mass turnover from precipitation, ice flow discharge and meltwater run-off. The high-melt year of 2012 applied in perpetuity yields an ice loss commitment of 782 ± 135 mm SLR, serving as an ominous prognosis for Greenland’s trajectory through a twenty-first century of warming. Greenland ice sheet melt is currently the largest single contributor to sea-level rise. This work combines observations and theory to show that Greenland ice sheet imbalance with recent climate (2000–2019) has already committed at least 3.3% ice volume loss, equivalent to 274 mm of global sea-level rise.

G reenland's ice budget deficit emerged after the 1980s from increases in surface meltwater run-off 1,2 and ice flow discharge from its tidewater sectors 3,4 . Yet, despite its importance for future sea-level rise (SLR), our capacity to accurately predict Greenland's response to climate change is hampered by process limitations in ice sheet models and their imprecise coupling to land, atmosphere and ocean boundaries [5][6][7] . Given these constraints, we pursue a complementary approach to obtain Greenland's thus far lower-bound-committed SLR contribution.
Our approach does not directly solve the transient ice flow equations but rather calculates the committed areal and volumetric changes incurred by the up-to-present ice geometrical disequilibrium with climate 8,9 . The method determines the ice extent and thickness perturbations required to bring the current ice sheet into equilibrium with surface mass balance (SMB). Changes in flow dynamics are implicitly accounted for by a glaciological power law scaling function that relates imposed areal changes in ice extent to an implied ice volume 10 . To account for marine-terminating sectors and tidewater outlets where the ablation area is truncated by iceberg calving, we introduce an effective ablation area treatment. For application to Greenland-including its peripheral ice masses-the essential empirical requirements are met with new multi-year inventories of: (1) tidewater glacier discharge 4 ; (2) SMB (that is, snowfall accumulation minus run-off) from observational reanalysis and regional climate data 11 ; and (3) the accumulation area ratio (AAR): the glacierized area with net annual mass gain divided by its total area, readily retrieved from optical satellite imagery 12,13 .
For grounded ice masses with an ablation area, the maximum snow line elevation at the end of each melt season marks the transition between the lower-elevation dark bare ice and the bright upperelevation snow accumulation areas. This equilibrium line and its corollary, the AAR, conveniently integrate the competing effects of surface mass loss from meltwater run-off and mass gain from snow accumulation. Minimum AAR each year demarcates hydrological years on a sector basis (Extended Data Fig. 1). By regressing annual AAR and mass balance, we obtain the statistical property of AAR in the condition of mass balance equilibrium (AAR 0 ) that is necessary for the current ice surface morphology to be in dynamic equilibrium with climate (Fig. 1). The ratio of the observed AAR to AAR 0 yields the fractional imbalance (α) that quantifies the area perturbation required for the ice mass to equilibrate its shape to an imposed climate shift away from that associated with AAR 0 (ref. 8 ). This disequilibrium approach exploits how climatically driven SMB perturbations are at least an order of magnitude faster than the associated dynamic adjustment of the ice mass 14 . The resulting derivation for the adjustment in ice volume (ΔV) and committed eustatic SLR follows glaciological scaling theory relating the glacierized area change to ice volume perturbation using a power law function 10 ( Fig. 1) with exponent (γ) (Methods).
While under the most up-to-date ice thickness and subglacial topography mapping 15 , Greenland's current ice sheet configuration implies an area-volume scaling exponent of γ = 1.24 that closely abides the theoretically derived value of 1.25 (ref. 10 ), we apply a linear exponent of 1 to avoid the mathematically intractable regional case in which some ice flux between adjacent flow sectors is inevitable. The choice of a linear exponent represents an absolute minimum committed loss, encompassing flow interaction between adjacent ice sheet sectors, since it accounts for how scaling techniques are best applied to ensembles of many ice masses 10,16 , which we accomplish by summing the volumes from 473 subregions of the ice sheet. While it is possible to scale the entire ice sheet with an exponent of 1.25, yielding volume changes roughly 20% higher than our results, this is an ill-posed inversion with potentially large random errors that grow exponentially with the size of the ice mass, a problem shared by both numerical models and scaling theory 17 . Our aggregate sums from many regions exploit the law of large numbers to dramatically reduce these errors 10 . As such, while an exponent of 1 underestimates the volume change for each region, such an approach guarantees a mathematically sound lower bound of the ice volume loss along with the resulting SLR while minimizing methodological errors.
While the method has previously only been applied to assess mountain glacier and ice cap disequilibria 8,9,18 , the theoretical basis and derivation of ice area-volume scaling analysis can be applied independently of size and area 8,9,18 . Hence, we apply the method to Greenland's entire glacierized area through summing over sectors to obtain lower-bound estimates of its committed mass loss and SLR resulting from its imbalance with recent climate. We calculate ice disequilibrium for three reference climate scenarios applied in perpetuity: 'recent' (2000-2019 average) climate to determine Greenland's already committed ice loss and then take the respective high-and low-melt years of 2012 and 2018 to assess potential future area and volume changes under extreme end-member climate states, with the proviso that no long-term reversal in climate warming trends are anticipated this century.

Results
Application of the average 2000-2019, hereafter 'recent' , climatology to Greenland's entire glacierized area of 1,783,090 km 2 gives an AAR/AAR 0 (α) disequilibrium with the current ice configuration corresponding with a 3.3 ± 0.8% committed area and volume loss. Taken in perpetuity, this imbalance with recent climate results in 59 ± 15 × 10 3 km 2 of committed retreat of Greenland's ice area, equivalent to 110 ± 27 × 10 3 km 3 of the ice sheet volume or 274 ± 68 mm of global eustatic SLR.
Over the 2000-2019 period, we found that surface ablation through meltwater run-off is the primary control on the trend and interannual variability of the Greenland ice sheet mass budget. The magnitude of the SMB trend is twice that of tidewater ice flow discharge (D) and the variability in SMB is an order of magnitude greater (Fig. 2). SMB variability is dominated by that from meltwater run-off (SMB versus R correlation: −0.855, 1 − P > 0.999), with the SMB versus precipitation (P) correlation = 0.528, 1 − P > 0.983. Despite run-off (R) averaging less than the ice flow discharge, its interannual variability exceeds that of D by a factor of ten. Indicative of an intensifying hydrological system, the ice sheet has been undergoing increased mass turnover, while paradoxically in a state of retreat and thinning 19,20 . The AAR to mass balance correlation holds independently for both tidewater and land-terminating sectors, demonstrating that the impact of increasing surface ablation and run-off is not limited to land-terminating sectors but prevails across the entire ice sheet given that 84% of its volume is marine-terminating (Extended Data Fig. 2).
The high-melt year of 2012 was driven by a persistent negative extreme in the North Atlantic Oscillation (NAO index = −1.6), delivering excess atmospheric heat to Greenland 21,22 , which caused enhanced and sustained surface ablation and meltwater run-off, for example 3,23 . Taking this high-melt year as an analogue for a . The ice sheet cross-sections illustrate the present and future ice area and volume adjustments necessary for dynamic equilibrium with shifting ELA and the equilibrium parameter ELA 0 , which is linearly proportional to AAR 0 . The associated area to volume scaling exponent γ, in this study taken conservatively as 1.0, is discussed in the Methods. projected late-twenty-first-century sustained warmer climate (Tedstone et al. 24 ), we found that the year 2012 AAR applied in perpetuity to Greenland commits a minimum of 10 ± 2% ice area (169 ± 29 × 10 3 km 2 ) decline equivalent to 314 ± 54 × 10 3 km 3 of the ice sheet volume or 782 ± 135 mm eustatic SLR (Table 1).
Conversely, summer 2018 was dominated by an extreme positive NAO (index = +1.5) with cold polar air drawn along Greenland's broad western flank, suppressing surface melt and resulting in the lowest surface mass loss recorded so far this century 3 (Fig. 1). Across virtually all sectors, this low surface ablation year had an AAR ≈ AAR 0 , indicating that the ice sheet was temporarily in a state of mass budget equilibrium. Western Greenland, which includes the large tidewater outlet glaciers of Jakobshavn Isbrae, Rink, Upernavik and Store (Mouginot et al. 19 ) still had an implicit adjustment with net mass loss despite predominantly cool surface climatology but which was insufficient to offset the mass gains from snowfall and suppressed surface melting.
The respectively high and low-melt years of 2012 and 2018 are useful to represent future extremes that are becoming the hallmark of Arctic amplified climate change [25][26][27] . Indeed, Greenland climate variability is observed to be increasing 28,29 , yet it is not well captured by global climate models 30 .
Given no area dependence on specific mass balance across the 473 sectors (Extended Data Fig. 3), we may obtain ice mass closure for the missing 1% of ice area arising from areas smaller than 30 km 2 not considered in this study due to a 1-km resolution constrain on the Moderate Resolution Imaging Spectroradiometer (MODIS) albedo grids by applying the mean specific imbalance for land-terminating sectors of −446 ± 309 kg m −2 y −1 . That residual area has a recent climate (years 2000-2019) imbalance of 3.4 ± 2.7 mm sea level equivalent.
Regional imbalance. Strong regional patterns of Greenland ice disequilibrium and committed SLR from the 2000 to 2019 climate are evident with the present-day imbalance of western Greenland being four times that of eastern sectors (Table 1). Such regional contrast arises from strong east-west gradients in snowfall accumulation 31 , coupled with high eastern topographic relief 15 and proximity to high south-eastern (SE) snowfall North Atlantic storms 31 that altogether result in a 30% steeper SE ice surface slope compared to the western ice sheet. Hence, the gentler western sectors of the ice sheet are generally more susceptible to recent surface warming and yield larger melt volumes amplified by the flatter elevation profile, thinner seasonal snow cover and a darker ablation area 32,33 . The year 2012 high-melt case resulted in a fourfold increase in disequilibrium from predominantly land-terminating south-western (SW) sectors and a 186 ± 30 mm SLR commitment, which is the highest of all Greenland sectors. Under the anomalously low-melt year 2018 surface climate, the SW incurs a mass budget surplus of 22 ± 13 mm equivalent sea-level drawdown. Although of a smaller magnitude than the sea-level drawdown from the north-eastern (NE) and central eastern (CE) regions, the magnitude of the difference between the 2012 and 2018 scenarios is greatest for the SW, highlighting the SW as a surface-climate-sensitivity hotspot region.
The predominantly marine-terminating north-western (NW) sector is most out of balance with the 2000-2019 climate (Table 1 and Fig. 3), committing 79 ± 12 mm equivalent SLR. The region is sensitive to climate warming 34 and meltwater run-off 35 compounded by its marine-terminating margin with a relatively high regional ice flow discharge. Under a perpetual high-melt year 2012, this region's mass loss commitment increases to 169 ± 28 mm equivalent SLR; under a perpetual 2018 climate, the sector is in mass budget equilibrium. In contrast to the NW, the relatively small SE sector, with just 8% of Greenland's ice volume, yields just 4% of the recent ice imbalance (Table 1) despite this region having the largest ice flow discharge; 20% of the total 36 . However, under the 2012 climate applied in perpetuity, the mass loss commitment from the SE increases substantially to 67 ± 8 mm equivalent SLR but has, like for the NW, a stable mass budget under the low-melt 2018 climate.
The relatively high elevation, cold and dry CE region is in mass budget equilibrium with the 2000-2019 climate (Table 1), a finding confirmed by recent satellite altimetry 3,20 . Despite rising air temperatures during the past two decades, a negative correspondence of mass balance with NAO across the region has delivered increased snowfall while limiting surface melt 37 . However, the CE region is vulnerable to future sustained mass loss as represented by the high-melt year of 2012. The CE region likewise exemplifies considerable climate sensitivity with an ice mass budget surplus under the cool year 2018 climate, yielding a substantial 44 ± 8 mm hypothetical sea-level drawdown.
Greenland's NE sector, despite including major outlet glaciers and ice streams that are undergoing considerable dynamic changes 19,38 , yields a relatively modest 24 ± 30 mm equivalent SLR under Greenland's recent climate but which rises to 65 ± 8 mm under a perpetual 2012 climate. Similarly, the northern (NO) sector with 14% of Greenland's total ice area, is also mainly marine-terminating but contributes just 3% to the recent mass imbalance due to much lower surface ablation across this region, increasing by half to 72 ± 30 mm SLR under the high-melt year 2012 climate.
Although land-terminating sectors comprise around a quarter of Greenland's glacierized area, they disproportionately account for 41% of Greenland's recent ice imbalance (Table 1) and demonstrate the dominance of surface ablation in currently committed  and grid-resolution dependence 45 . However, for application to Greenland, these have not yet been fully integrated with sufficiently realistic atmosphere/land/ocean forcing data nor with sufficient treatment of amplifying processes to yield a sensitive sub-centennial response to intensifying climate change 6,46 . The Greenland ice loss commitment presented in this study is complementary to process-based models, yet it is grounded in observations and theory, albeit lacking a response timescale. Our assessed Greenland's ice disequilibrium with 2000-2019 surface climate that commits the equivalent of 274 ± 68 mm SLR aligns with expert judgement predicting a mean committed Greenland SLR contribution of 330-490 mm by 2100 (ref. 47 ). Moreover, application of the anomalous 2012 melt year in perpetuity taken as a representative analogue for sustained later-this-century climate yields a Greenland SLR commitment of 782 ± 135 mm, which also resonates with expert judgement of metre-scale SLR under the unabated SSP58 future climate scenario 47 .
The conservative lower bound and error budget (Methods) encompass how our estimates of committed SLR neglect potential future changes to ice discharge through fast-flowing outlet glaciers that are not the result of geometric adjustment to climate but result from changes in fjord circulation, oceanic heat delivery, calving front dynamics and/or subglacial conditions. By demonstrating how the low (2018) and high (2012) melt years encompass a state of Greenland ice sheet stability on the one hand, or extreme ice loss amplified by multiple feedbacks on the other, we underscore the high sensitivity and susceptibility of Greenland's terrestrial ice to extremes in atmospheric circulation that may not be implicit in global climate model outputs, for example, Coupled Model Intercomparison Project Phase 6 (ref. 48 ). Furthermore, the large interannual Greenland surface climate variability represented by 2012 and 2018 is obscured when presentation of mass balance accounting is made as multi-method aggregates with a time-invariant uncertainty envelope 3 . Rather, given that observations demonstrate increases in interannual variability 28,29 , we contend that accounting for time-varying magnitudes of extremes may be essential in obtaining a realistic envelope of future ice response states.
Newly emerging processes are driving rapid ice sheet response: tidewater glacier acceleration and destabilization by submarine melting [49][50][51] ; loss of floating ice shelves 52 ; accelerating interior motion from increased melt and rainfall 53 ; enhanced basal thawing due to hydraulically released latent heat and viscous warming 54 ; amplified surface melt run-off due to bio-albedo darkening 55 ; and impermeable firn layers 56 amplified by ice sheet surface hypsometry 32,33 . Given the breadth and potency of those processes, we contend that known physical mechanisms can deliver most of the committed ice volume loss from Greenland's disequilibrium with its recent climate within this century. Nevertheless, we underscore that a SLR of at least 274 ± 68 mm is already committed, regardless of future climate warming scenarios.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41558-022-01441-2.

Methods
AAR and imbalance. The equilibrium line altitude (ELA; Fig. 1) separates an ice sheet's higher elevation accumulation area from its lower elevation ablation area. The AAR is the fraction of the ice sheet's total area that lies within the accumulation area. In high-melt and/or low snowfall years, the ELA moves to higher elevations and the AAR decreases. The opposite occurs in years with reduced melting and/or increased snowfall.
The AAR is retrieved using optical satellite imagery, which distinguishes the high brightness of snow cover from darker exposed glacier ice 12,57 . From the observed year-to-year variation of AAR, we exploited how the ratio between two AAR states measures an ice sheet's committed area and volume losses 10 , using scaling techniques elaborated in the following sections. One AAR state is the observed recent climate (2000-2019), established by averaging the observed time-varying annual AAR values (AAR 2000-2019 ). The second AAR state is that which is necessary for an ice sheet to maintain equilibrium (AAR 0 ); it is estimated from a regression of annual AAR and ice sheet mass balance values 9 . The observed average AAR being less than AAR 0 indicates that the accumulation area has shrunk (Fig. 1). Thus, the ice sheet is overextended and must retreat to reach a new equilibrium with the recent climate. This retreat from overextension occurs after changes in the SMB that happen quickly relative to changes in the area of the glacier, which depends on slower ice dynamic adjustment. Typical e-folding response times for glacier and ice sheet flow range from 10 to 1,000 s of years depending on the ice volume, among other factors 14,58,59 . The ratio α = AAR/AAR 0 is directly proportional to the area change needed for an ice mass to adjust to a new equilibrium 10,60 . The corresponding volume change (ΔV = (α γ − 1)V) involves the volume-area scaling exponent γ. In this study, we took γ equal to 1, appropriate for ice caps and ice sheets and allowing flow interactions between adjacent sectors 10 . The change in eustatic sea level is ΔV/β where β is 362 billion tonnes of water equivalent per mm global sea-level change.
Area-volume scaling. While we found that an exponent of 1.24 is evident from the observations of area versus volume of the ice sheet 15 , closely matching the theoretically derived value of 1.25 (ref. 10 ), we chose a linear exponent of 1 to avoid the mathematically intractable regional case in which flow interaction between adjacent flow sectors is inevitable. The choice of a linear exponent represents an absolute minimum committed loss that encompasses any flow interaction between adjacent ice sheet flow sectors. The linearization has two significant benefits. First, combining the aggregate volumes of the ice sheet's subregions in this manner is not possible with a non-linear exponent 10 but is possible with a linear exponent. Second, scaling techniques are best applied to ensembles of many ice masses 10,16 , which we accomplished by summing the volumes from many subregions of the ice sheet. While it is possible to scale the entire ice sheet with an exponent of 1.25 (yielding volumes roughly 20% higher than our results), this is an ill-posed inversion with potentially large random errors that grow exponentially with the size of the ice mass, a problem shared by both numerical models and scaling 17 ; our aggregate sums from a number of smaller regions exploit the law of large numbers to dramatically reduce such errors 10 . As such, while an exponent of 1 underestimates the volume of each region by roughly 20%, our approach guarantees a mathematically sound lower bound of the volume (mass) lost and the resulting SLR, while minimizing errors.

Separation of bare ice and snow-covered areas.
Over the Greenland ice sheet and peripheral ice caps, with areas greater than 30 km 2 , AARs are measured daily via a bare-ice classification made for the full 20-year (2000-2019) catalogue of gap-filled 13 v.6 NASA Terra satellite MODIS MOD10A1 albedo observations. The representative bare-ice albedo value (0.565 ± 0.109) computed from 105 station years of the Programme for the Monitoring of the Greenland Ice Sheet (PROMICE) climate station data was used to separate the ablation and accumulation areas 61 .
The AARs incorporate the BedMachine v3 (ref. 15 ) ice mask. The 473 individual flow sectors were defined by Mouginot et al. 19 .
Derivation of Greenland ice imbalance. The parameter α = AAR/AAR 0 is an accurate measure of glaciated catchment imbalance 9 that is independent of volume (V)-area (A) scaling 10 . It can be applied to any glaciated catchment at any scale, including the entire Greenland ice sheet, specific sectors of the ice sheet, partial regions or individual catchments and partial ice caps. Conversely, the derivation of volume change (∆V = (α γ − 1)V) depends on the well-established scaling relationship V = cA γ where c is a scaling parameter that is unused and unnecessary in this analysis. For ice caps and ice sheets, γ = 1, a linearized solution, is the appropriate value in V = cA γ to encompass the situation where flow interaction between adjacent sectors can occur. This linearization gives a lower-bound estimate of equilibrium volume change and hence committed future SLR, without requiring dynamic independence. Thus, ∆V = (α − 1)V.

Committed area loss.
We assessed committed area loss using the definition 8 of α as: where ΔA is the committed area loss and A is the current ice area. Therefore, Substituting our values, gives: ΔA = 1, 783, 090 km 2 (1.035 ± 0.012 − 1) = 1, 783, 090 (0.035 ± 0.012) = 62, 408 ± 21, 397 km 2 The value of α is as in Figs. 2 and 3 but for all Greenland ice, that is, not regionally. The Greenland ice area is from BedMachine v.3 data 15 . The ± values result from the error propagation described below.
The linear area-volume scaling exponent used in this study has no numerical influence.
Effective ablation area. Tidewater sector ablation areas are truncated by calving, a process not accounted for by previous theory 10 . Without calving, tidewater sectors would have larger ablation areas and smaller AAR 0 . Likewise, without calving, land-terminating sectors and adjacent tidewater sectors would have similar AAR 0 on average because they are exposed to similar surface climates. In this study, to unify land-terminating and tidewater sectors, additional ablation area is added to each tidewater sector so that each one's equilibrium AAR 0 is equivalent to the observed average AAR 0 = 0.780 ± 0.02 for all Greenland land-terminating sectors. This effective ablation area increases the assessed imbalance from tidewater glaciers by 25 ± 4% on average among regional basins.

Derivation of committed volume loss.
A detailed derivation of the committed volume loss ∆V can be found in Bahr et al. 8 and is abbreviated in this article. Let P A = ΔA/A and P V = ΔV/V be the respective fractional changes in ice area and volume. From volume-area scaling: Substituting the definition of α (above) gives, From scaling theory, the exponent γ is a constant that cannot be changed by boundary conditions such as rapid sliding, retrograde bed slopes, and so on 10 . Only the multiplicative scaling parameter c is changed by boundary conditions, but in the above derivation, c has factored away and plays no role in this study's analysis.
Using our expected errors for each of the Greenland regions, the average of the dimensionless P V has an expected error of 0.0126 (or 1.3%), that is, the ±104 mm in the all-region SLR total (Table 1).
Multiplying the mean of P V (376 mm × 362 Gt mm −1 eustatic SLR) by Greenland's total ice volume gives the predicted change in Greenland's ice volume of 136,112 km 3 .
Non-use of the scaling parameter c. The scaling parameter c can be eliminated altogether when estimating global (or regional) changes in aggregate ice volume 8,62,63 . After derivation 10 , the total regional ice volume V T is estimated from each sector's current volume V before volume changes are projected into the future. In general, a scaling analysis applies to a population of glaciers, ice caps and ice sheets. With the linear volume-area scaling exponent discussed above, this translates to a population of regions or 'fractional parts' of an ice sheet. Section 8.5 in Bahr et al. 10 includes the specifics. SMB model data. SMB is defined in this study as precipitation minus sublimation minus run-off (note that Cogley et al. 64 refer to this as climatic mass balance). SMB is prescribed daily at catchment scale after two regional climate models (RCMs) 11 : (1) Modèle Atmosphérique Régional (MAR) v.3.10; and (2) Regional Atmospheric Climate Model (RACMO) v.2.3p2.
MAR with 24 vertical levels was forced by ERA5 until December 2019 at a MAR horizontal resolution of 7.5 km including a snow and ice surface scheme based on the CROCUS snow model. The model allows meltwater refreezing and snow metamorphosis, influencing the transformation of snow to ice and the surface albedo, which are calibrated to agree with the satellite-derived melt extent over 1979-2009. The MAR outputs are downscaled to a 1 km grid 11 .
With 40 vertical-atmospheric-level RACMO2 is run at 5.5 km horizontal resolution and forced by a combination of climate reanalyses including ERA- 40 (1958-1978), ERA-Interim (1979-2018) and ERA5 (2019) 57 . The model incorporates a multilayer snow module that simulates melt, liquid water percolation and retention, refreezing and run-off, and accounts for dry snow densification. RACMO2 also implements an albedo scheme that calculates snow albedo based on prognostic snow grain size, cloud optical thickness, solar zenith angle and soot concentration in snow. Ice albedo is prescribed from a 500 m MODIS 16 d Albedo product (MCD43A3), as the 5% lowest surface albedo records for the period 2000-2015, clipped between 0.30 for bare ice and 0.55 for bright ice covered by perennial firn in the accumulation area. The snowpack was initialized from September 1957 using temperature and density profiles derived from the offline firn densification model (IMAU-FDM). The SMB and individual components are statistically downscaled to a 1 km horizontal resolution 11 .
Sector scale mass balance. For land-terminating sectors, the total mass balance equals the sum of the SMB. The SMB for each individual flow sector or region 19 (Fig. 2) is from the two regional climate models. (The SMB model details are shown above.) The SMB data sources are applied independently to generate twice as many end-members in our final ensemble (see error propagation below). Ice flow D (ref. 4 ) data enable defining total mass balance for tidewater sectors as SMB-D in gigatonne (Gt) units. The SMB and ice flow D data are totalled for hydrological years according to dates between annual minimum AAR (Extended Data Fig. 1). Sector scale ice flow D data are averaged for the same hydrological years as for the SMB.
Scaling for full ice mass closure. Tidewater SMB is consistently positive, as it should be for all the basins that are not extremely far out of mass balance, with total mass balance occurring after iceberg calving and underwater melting (Fig. 1). Total specific mass balance (mass balance divided by area) is negative, which is indicative of observed mass budget deficit and sea-level contribution 3 . Negative average SMB for land-terminating sectors was also evident in observations. Error propagation. To yield uncertainty estimates, a range in values for mass balance, snowline and AAR parameters was propagated through the mass balance imbalance calculation (see Fig. 1) for 2ach catchment and hydrological year. The results among 473 sectors were averaged and their s.d. was taken as the uncertainty. The selection of parameters was after the following cases with 3 × 3 × 2 × 3 yielding a 54-member ensemble: (1) adjustment of equilibrium AAR (AAR 0 ) for tidewater sectors to a representative land-terminating sector value (from 0.76 to 0.80), that is, 0.78 ± 0.02, three cases ( Fig. 1 and Extended Data Fig. 2); (2) albedo at firnline, three cases. By varying the albedo threshold (0.565 ± 0.109) used to separate snow from bare ice in the NASA MODIS albedo imagery 62 , a mass balance variability is introduced for each sector parameter by ± half its uncertainty (0.054). This introduced range has the effect of varying the dates used for obtaining hydrological years for SMB and ice flow D; (3) SMB, two cases. The SMB for each catchment and hydrological year varies depending on whether it comes from MAR v.3.10 or RACMO2.3p2 SMB output downscaled to 1 km (ref. 11 ); (4) ice flow D, three cases. The parameter D for each tidewater sector and hydrological year is varied after the stated ±10% uncertainty 4 .

Mass balance validation.
To validate our imbalance results, monthly observations of changes in the Earth's gravity field from the Gravity Recovery and Climate Experiment (GRACE) (2002-2017) and GRACE Follow-On (GRACE-FO) (2018 to present) satellites were used to derive annual mass balances for the Greenland ice sheet (including peripheral glaciers and tundra) 65 . When considering short timescales, these variations are directly related to mass redistribution at the Earth's surface. Using an ensemble of several level 2 GRACE spherical harmonics, monthly mass changes of the Greenland ice sheet were inverted with an error-weighted ensemble average to reduce uncorrelated noise and mass changes due to glacial isostatic adjustment corrected for using a spherical Earth model. Annual mass balances were derived by differencing successive September values, with missing September months filled using spline interpolation. For the 2016/2017 and 2017/2018 mass balances, we used June-June differences because the one-year gap between the GRACE and GRACE-FO missions would result in excessively large uncertainties in the spline interpolation.
A minor amount of calving from ice caps, not resolved in this study, means that the actual Greenland land ice mass balance may be lower. Visual inspection reveals that 12 of the 20 largest ice caps have tidewater glaciers. However, total mass balance from this study during the period of independent satellite gravimetry (2003-2019), that is, from GRACE and GRACE-FO (−270 ± 86 billion tonnes per year), is equivalent to the average (−272 ± 45 billion tonnes per year) from this study (correlation = 0.640, P = 0.005).

Data availability
The scripts and output data generated during and/or analysed as part of the submitted study are available from https://doi.org/10.22008/FK2/D5JEZ0 66 . The MAR and RACMO regional climate model output is available from the respective lead scientists.

Code availability
The scripts and output data generated during and/or analysed as part of the submitted study are available from https://doi.org/10.22008/FK2/D5JEZ0 66 .