Articles | Volume 13, issue 10
Earth Syst. Sci. Data, 13, 5001–5025, 2021

Special issue: Extreme environment datasets for the three poles

Earth Syst. Sci. Data, 13, 5001–5025, 2021

Data description paper 29 Oct 2021

Data description paper | 29 Oct 2021

Greenland ice sheet mass balance from 1840 through next week

Greenland ice sheet mass balance from 1840 through next week
Kenneth D. Mankoff1, Xavier Fettweis2, Peter L. Langen3, Martin Stendel4, Kristian K. Kjeldsen1, Nanna B. Karlsson1, Brice Noël5, Michiel R. van den Broeke5, Anne Solgaard1, William Colgan1, Jason E. Box1, Sebastian B. Simonsen6, Michalea D. King7, Andreas P. Ahlstrøm1, Signe Bech Andersen1, and Robert S. Fausto1 Kenneth D. Mankoff et al.
  • 1Department of Glaciology and Climate, Geological Survey of Denmark and Greenland (GEUS), Copenhagen, Denmark
  • 2SPHERES research unit, Department of Geography, University of Liège, Liège, Belgium
  • 3Department of Environmental Science, iClimate, Aarhus University, Roskilde, Denmark
  • 4Danish Meteorological Institute (DMI), Copenhagen, Denmark
  • 5Institute for Marine and Atmospheric Research, Utrecht University, the Netherlands
  • 6Geodesy and Earth Observation, DTU Space, Technical University of Denmark, Lyngby, Denmark
  • 7Polar Science Center, University of Washington, Seattle, WA, USA

Correspondence: Ken Mankoff (


The mass of the Greenland ice sheet is declining as mass gain from snow accumulation is exceeded by mass loss from surface meltwater runoff, marine-terminating glacier calving and submarine melting, and basal melting. Here we use the input–output (IO) method to estimate mass change from 1840 through next week. Surface mass balance (SMB) gains and losses come from a semi-empirical SMB model from 1840 through 1985 and three regional climate models (RCMs; HIRHAM/HARMONIE, Modèle Atmosphérique Régional – MAR, and RACMO – Regional Atmospheric Climate MOdel) from 1986 through next week. Additional non-SMB losses come from a marine-terminating glacier ice discharge product and a basal mass balance model. From these products we provide an annual estimate of Greenland ice sheet mass balance from 1840 through 1985 and a daily estimate at sector and region scale from 1986 through next week. This product updates daily and is the first IO product to include the basal mass balance which is a source of an additional ∼24Gt yr−1 of mass loss. Our results demonstrate an accelerating ice-sheet-scale mass loss and general agreement (coefficient of determination, r2, ranges from 0.62 to 0.94) among six other products, including gravitational, volume, and other IO mass balance estimates. Results from this study are available at (Mankoff et al.2021).

1 Introduction

Over the past several decades, mass loss from the Greenland ice sheet has increased (Khan et al.2015; IMBIE Team2019). Different processes dominate the regional mass loss of the ice sheet, and their relative contribution has fluctuated in time (Mouginot and Rignot2019). For example, in the 1970s nearly all sectors gained mass due to positive surface mass balance (SMB), except the northwestern sector, where discharge losses dominated. More recently, in the 2010s, all sectors lost mass, with some sectors losing mass almost entirely via negative SMB and others primarily due to discharge (Fig. 1).

Figure 1Annual mass balance (black lines), surface mass balance (blue lines), and discharge plus basal mass balance (dashed grey) in Gt yr−1 for each of the seven Mouginot and Rignot (2019) regions. The map shows both the named regions (Mouginot and Rignot2019) and the numbered sectors (Zwally et al.2012). Discharge gates are marked in black. Only recent (post-1986) data are shown because reconstructed data are not separated into regions or sectors. Next week is defined as 1 November 2021 based on the date this document was compiled.

There are three common methods for estimating mass balance – changes in gravity (Barletta et al.2013; Groh et al.2019; IMBIE Team2019; Velicogna et al.2020), changes in volume (Simonsen et al.2021a; Sørensen et al.2011; Zwally and Giovinetto2011; Sasgen et al.2012; Smith et al.2020), and the input–output (IO) method (Colgan et al.2019; Mouginot et al.2019; Rignot et al.2019; King et al.2020). Each provides some estimate of where, when, and how the mass is lost or gained, and each method has some limitations. The gravity mass balance (GMB) estimate has low ∼100km spatial resolution (where), monthly temporal resolution (when), and little information on the processes contributing to changes in mass balance components (how). The volume change (VC) mass balance estimate has ∼1km spatial resolution (where), often provided at annual or multi-year temporal resolution (when), and little information on the driving processes (how).

The IO method has a complex spatial resolution (where). The inputs typically come from regional climate models (RCMs) which can reach a spatial resolution of up to 1 km. However, that spatial resolution is generally reduced in the final output to sector or region scale – typically higher than GMB but now lower resolution than VC. The IO temporal resolution (when) is limited by ice velocity data updates, which for the past several years occur every 12 d year-round after the launch of the Sentinel missions (Solgaard et al.2021). The primary issue with the IO method is unknown ice thickness in some locations (e.g., Mankoff et al.2020b). Finally, the IO method can provide insight into the processes (how) by distinguishing between changes caused by SMB (which may be due to changes in positive and/or negative SMB components) vs. changes in other mass loss terms (e.g., calving). Our IO method is also the first IO product to include the basal mass balance (Karlsson et al.2021) – a term implicitly included in the GMB and VC methods but neglected by all previous IO estimates.

In this work we introduce the new Programme for Monitoring of the Greenland ice sheet (PROMICE) Greenland ice sheet mass balance data set based on the IO method, updating the previous product from Colgan et al. (2019). We use the SMB field from one empirical model from 1840 through 1985 and three RCMs from 1986 onward. The combined SMB field used here is comprised of positive SMB terms (precipitation in the form of snowfall, rainfall, condensation/riming, and snow drift deposition) and negative SMB terms (surface melt, evaporation, sublimation, and snow drift erosion). We also use the basal mass balance and an estimate of dynamic ice discharge. Spatial resolution is effectively per sector (Zwally et al.2012) or region (Mouginot and Rignot2019). Temporal resolution is annual from 1840 through 1985 and effectively daily since 1986 – the RCM fields are updated daily and forecasted through next week, and the discharge at marine-terminating glaciers is updated every 12 d with ∼12d resolution, interpolated to daily, and forecasted using historical and seasonal trends through next week. Thus, this study provides a daily-updating estimate of Greenland mass changes from 1840 through next week.

2 Terminology

We use the following terminology throughout the document.

  • This Study refers to the new results presented in this study.

  • Recent refers to the new 1986 through next week daily temporal resolution data at regional and sector scales.

  • Reconstructed refers to the adjusted Kjeldsen et al. (2015) annual temporal resolution data at ice sheet scale used to extend this product from 1986 back through 1840. The 1986 through 2012 portion of the Kjeldsen et al. (2015) data set is used only to adjust the reconstructed data and is then discarded.

  • ROI (region of interest) refers to one or more of the ice sheet sectors or regions (Fig. 1).

  • Sector refers to one of the Zwally et al. (2012) sectors (Fig. 1), expanded here to cover the RCM ice domains which exist slightly outside these sectors in some locations.

  • Region refers to the Mouginot and Rignot (2019) regions (Fig. 1), expanded here to cover the RCM ice domains.

  • SMB is the surface mass balance from an RCM or the average of multiple RCM SMBs. The use should be clear from the context.

  • D is solid ice discharge. It includes both calving and submarine melting at marine-terminating glaciers.

  • BMB is the basal mass balance. It comes from geothermal flux (BMBGF), frictional heating from ice velocity (BMBfriction), and viscous heat dissipation (BMBVHD).

  • MB is the total mass balance including the BMB term (Eq. 3).

  • MB* is the mass balance not including the BMB term (Eq. 4).

  • HIRHAM/HARMONIE, Modèle Atmosphérique Régional (MAR), and Regional Atmospheric Climate MOdel (RACMO) refer to the RCMs, which only provide SMB and runoff in the case of MAR. However, when referencing the different MB products, we use, for example, “MAR MB” rather than repeatedly explicitly stating “MB derived from MAR SMB minus BMB and D”. The use should be clear from the context.

3 Product description

The output of this work is two NetCDF files and two CSV files containing a time series of mass balance and the components used to calculate mass balance. The only difference between the two NetCDF files is the region of interest (ROI) – one for Zwally et al. (2012) sectors and one for Mouginot and Rignot (2019) regions. Each NetCDF file includes the ice sheet mass balance (MB), MB per ROI (sector or region), MB per ROI per RCM, ice sheet SMB, SMB per ROI, ice sheet discharge (D), D per ROI, ice sheet BMB, and BMB per ROI. The CSV files contain a copy of the ice-sheet-summed data, one daily and one annual.

Figure 2Mass balance and its major components. (a) Annual average surface mass balance (blue line), discharge (gray dashed), and their mass balance sum (black line). Here the discharge and basal mass balance (D + BMB) are shown with sign inverted (e.g., -1×(D+BMB)). (b) Same data at daily resolution and limited to 2019 and 2020.


An example of the output is shown in Fig. 2a, which shows mass balance for the entire Greenland ice sheet in addition to SMB and D at annual resolution. Figure 2b shows an example of 2 years at daily temporal resolution. The ice-sheet-wide product includes data from 1840 through next week, but the sector and region-scale products only include data from 1986 through next week, because the 1840 through 1985 reconstructed only exists at ice sheet scale (Fig. 1).

4 Data sources

This section introduces data products that exist prior to and are external to this work (Table 1). In the following Methods section we introduce both the intermediate products we generate using these data sources and the final product that is the output of This Study.

Kjeldsen et al. (2015)Box et al. (2021)Kjeldsen et al. (2015)Box et al. (2021)Langen et al. (2017)Fettweis et al. (2020)Noël et al. (2019)Mankoff et al. (2020b)Mankoff and Solgaard (2020)Karlsson et al. (2021)Karlsson (2021)Kjeldsen et al. (2015)Box et al. (2021)Fettweis et al. (2020)

Table 1Summary of data products used as inputs to This Study.

Download Print Version | Download XLSX

The inputs to this work are the recent SMB fields from the three RCMs, the recent discharge from Mankoff et al. (2020b) (data: Mankoff and Solgaard2020), and the recent basal mass balance fields, of which BMBGF and BMBfriction are direct outputs from Karlsson et al. (2021) (data: Karlsson2021), but the BMBVHD calculations are redone here (see Methods Sect. 5.3) using the MAR runoff field. The reconstructed data (pre-1986) are surface mass balance and discharge from Kjeldsen et al. (2015) (data: Box et al.2021) but adjusted here using the overlapping period (see Methods Sect. 5.4) and runoff from Kjeldsen et al. (2015) (data: Box et al.2021) as a proxy and scaled for BMBVHD (see Methods Sect. 5.3).

4.1 Surface mass balance

We use one reconstructed SMB from 1840 through 1985 and three recent SMBs from 1986 through last month (HIRHAM/HARMONIE, MAR, and RACMO), two through yesterday (HIRHAM/HARMONIE and MAR) and one through next week (MAR).


The HIRHAM/HARMONIE product from the Danmarks Meteorologiske Institut (Danish Meteorological Institute; DMI) is based on an offline subsurface firn/SMB model (Langen et al.2017), which is forced with surface fluxes of energy (turbulent and downward-radiative) and mass (snow, rain, evaporation, and sublimation). These surface fluxes are derived from the HIRHAM5 regional climate model for the reconstructed part of the simulation and from DMI's operational numerical weather forecast model HARMONIE (Iceland–Greenland domain “B”, which covers Iceland, Greenland, and the adjacent seas) for the real-time part. HIRHAM5 is used until 31 August 2017, after which HARMONIE is used.

The HIRHAM5 regional climate model (Christensen et al.2007) combines the dynamical core of the HIRLAM7 numerical weather forecasting model (Eerola2006) with physics schemes from the ECHAM5 general circulation model (Roeckner et al.2003). In the Greenland setup employed here (Lucas-Picher et al.2012), it has a horizontal resolution of 0.05×0.05 on a rotated pole grid (corresponding to 5.5 km resolution) and 31 atmospheric levels. It is forced at 6 h intervals on the lateral boundaries with horizontal wind vectors, temperature, and specific humidity from the ERA-Interim reanalysis (Dee et al.2011). ERA-Interim sea surface temperatures and sea ice concentration are prescribed in ocean grid points. Surface fluxes from HIRHAM5 are passed to the offline subsurface model.

The offline subsurface model was developed to improve firn details for the HIRHAM5 experiments (Langen et al.2017). The subsurface consists of 32 layers with time-varying fractions of snow, ice, and liquid water. Layer thicknesses increase with depth from 6.5 cm water equivalent (w.e.) at the top to 9.2 m w.e. at the bottom, giving a full model depth of 60 m w.e. The processes governing the firn evolution include snow densification, varying hydraulic conductivity, irreducible water saturation and other effects on snow liquid water percolation, and retention. Runoff is calculated from liquid water in excess of the irreducible saturation with a characteristic local timescale that depends on surface slope (Zuo and Oerlemans1996; Lefebre2003). The offline subsurface model is run on the HIRHAM5 5.5 km grid.

HARMONIE (Bengtsson et al.2017) is a nonhydrostatic model in terrain-following sigma coordinates based on the fully compressible Euler equations (Simmons and Burridge1981; Laprise1992). HARMONIE is run at 2.5 km horizontal resolution and with 65 vertical levels. Compared to previous model versions, upper-air 3D variational data assimilation of satellite wind and radiance data, and radio occultation data, radiosonde, aircraft, and surface observations are incorporated. This greatly improves the number of observations in the model, as in situ observations from ground stations and radiosondes only make up approximately 20 % of observations in Greenland (Wang and Randriamampianina2021; Yang et al.2018). The model is driven at the boundaries with European Centre for Medium-Range Weather Forecasts (ECMWF) high-resolution data at 9 km resolution. The 2.5 km HARMONIE output is regridded to the 5.5 km HIRHAM grid before input to the offline subsurface model. The HIRHAM5 and offline models both employ the Citterio and Ahlstrøm (2013) ice mask interpolated to the 5.5 km grid.

4.1.2 MAR

The MAR RCM has been developed by the University of Liège (Belgium) with a focus on the polar regions (Fettweis et al.2020). The MAR atmosphere module (Gallée and Schayes1994) is fully coupled with the soil–ice–snow energy balance vegetation model SISVAT (Gallée et al.2001) simulating the evolution of the first 30 m of snow or ice over the ice sheet with the help of 30 snow layers (with time-varying thickness) or the first 10 m of soil over the tundra area. At its lateral boundary, MAR is forced at 6 h intervals by ERA5 reanalysis and runs at 20 km resolution. The snowpack was initialized in 1950 from a former MARv3.11-based simulation. Its snow model is based on a former version of the CROCUS snow model (Vionnet et al.2012) dealing with all the snowpack processes, including the meltwater retention, transformation of melting snow and grain size, compaction of snow, formation of ice lenses impacting meltwater penetration, warming of the snowpack from rainfall, and complex snow/bare ice albedo. MAR uses the Greenland Ice Mapping Project (GIMP) ice sheet mask and ice sheet topography (Howat et al.2014).

We use MAR version 3.12. With respect to version 3.9, intensively validated over Greenland (Fettweis et al.2020), or the 20 km-based MARv3.10 setup used in Tedesco and Fettweis (2020), MARv3.12 now uses the common polar stereographic projection EPSG 3413. With respect to MARv3.11, fully described in Amory et al. (2021), MARv3.12 ensures now the full conservation of water mass into both soil and snowpack at each time step, takes into account the geographical projection deformations in its advection scheme, better deals with the snow/rain temperature limit with a continuous temperature threshold between 0 and 2 C, increases the evaporation above snow thanks to a saturated humidity computation in SISVAT adapted to freezing temperatures, disallows melt below the 30 m of the resolved snowpack, and includes small improvements and bug fixes with the aim of improving the evaluation of MAR (with both in situ and satellite products) as presented in Fettweis et al. (2020) in addition to small computer time improvements in the parallelization of its code.

In addition to providing SMB, MAR also provides daily runoff over both permanent ice and tundra area. The ice runoff is used for the daily BMBVHD estimate (Sect. 5.3).

As the recent SMB decrease (successfully evaluated with GRACE-based estimates in Fettweis et al.2020) has been fully driven by the increase in runoff (Sasgen et al.2020), we assume the same degree of accuracy between SMB simulated by MAR (evaluated with the PROMICE SMB database in Fettweis et al.2020) and the runoff simulated by MAR.

Weather-forecasted SMB. To provide a real-time state of the Greenland ice sheet, MAR is forced automatically every day by the run of 00:00 UTC from the Global Forecast System (GFS) model providing weather forecasting initialized by the snowpack behaviors of the MAR run from the previous day. This continuous GFS-forced time series (without any reinitialization of MAR) provides SMB and runoff estimates between the period covered by ERA5 and the next 7 d. At the end of each day, ERA5 is used to update the GFS-forced MAR time series until about 5 d before the current date and to provide a homogeneous ERA5-forced MAR times series from 1950 to a few days before the current date. We use both the forecasted SMB and forecasted runoff (for BMBVHD) fields.

4.1.3 RACMO

RACMO v2.3p2 has been developed at the Koninklijk Nederlands Meteorologisch Instituut (Royal Netherlands Meteorological Institute; KNMI). It incorporates the dynamical core of the High-Resolution Limited Area Model (HIRLAM) and the physics parameterizations of the ECMWF Integrated Forecast System cycle CY33r1. A polar version (p) of RACMO has been developed at the Institute for Marine and Atmospheric research of Utrecht University (UU-IMAU) to assess the surface mass balance of glaciated surfaces. The current version RACMO2.3p2 has been described in detail in Noël et al. (2018), and here we repeat the main characteristics.

The ice sheet has an extensive dry interior snow zone, a relatively narrow runoff zone along the low-lying margins, and a percolation zone of varying width in between. To capture these processes on first order, the original single-layer snow model in RACMO has been replaced by a 40-layer snow scheme that includes expressions for dry snow densification and a simple tipping bucket scheme to simulate meltwater percolation, retention, refreezing, and runoff (Ettema et al.2010). The snow layers were initialized in September 1957 using temperature and density from a previous run with the offline IMAU Firn Densification Model (Ligtenberg et al.2018). To simulate drifting snow transport and sublimation, Lenaerts et al. (2012) implemented a drifting snow scheme. Snow albedo depends on snow grain size, cloud optical thickness, solar zenith angle, and impurity content (van Angelen et al.2012). Bare ice albedo is assumed constant and estimated as the 5th percentile value of albedo time series (2000–2015) from the 500 m-resolution MODIS 16 d albedo product (MCD43A3). Minimum/maximum values of 0.30/0.55 are applied to the bare ice albedo, representing ice with high-/low-impurity content (cryoconite, algae).

To simulate as accurately as possible the contemporary climate and surface mass balance of the ice sheet, the following boundary conditions have been applied. The glacier ice mask and surface topography have been downsampled from the 90 m-resolution Greenland Ice Mapping Project (GIMP) digital elevation model (DEM) (Howat et al.2014). At the lateral boundaries, model temperature, specific humidity, pressure, and horizontal wind components at the 40 vertical model levels are relaxed towards 6-hourly ECMWF reanalysis (ERA) data. For this we use ERA-40 between 1958 and 1978 (Uppala et al.2005), ERA-Interim between 1979 and 1989 (Dee et al.2011), and ERA-5 between 1990 and 2020 (Hersbach et al.2020). The relaxation zone is 24 grid cells (∼130km) wide to ensure a smooth transition to the domain interior. This run has active upper-atmosphere relaxation (van de Berg and Medley2016). Over glaciated grid points, surface aerodynamic roughness is assumed constant for snow (1 mm) and ice (5 mm). In this run, RACMO2.3p2 has 5.5 km horizontal resolution over Greenland and the adjacent oceans and land masses, but it was found previously that this is insufficient to resolve the many narrow outlet glaciers. The 5.5 km product is therefore statistically downscaled onto a 1 km grid sampled from the GIMP DEM (Noël et al.2019), employing corrections for biases in elevation and bare ice albedo using a MODIS albedo product at 1 km resolution (Noël et al.2016).

4.1.4 Reconstruction

The Kjeldsen et al. (2015) 173-year (1840 through 2012) mass balance reconstruction is based on the Box (2013) 171-year (1840 through 2010) statistical reconstruction. Kjeldsen et al. (2015) add a more sophisticated meltwater retention scheme (Pfeffer et al.1991), weighting of in situ records in their contribution to the estimated value, and dispersal of annual accumulation to monthly and extend the reconstruction in time through 2012.

The Box (2013) 171-year (1840–2010) reconstruction is developed from linear regression parameters that describe the least squares regression between (a) spatially discontinuous in situ monthly air temperature records (Cappelen et al.2011, 2006; Cappelen2001; Vinther et al.2006) or firn/ice cores (Box et al.2013) and (b) spatially continuous outputs from the regional climate model RACMO version 2.1 (Ettema et al.2010). A 43-year overlap period (1960 through 2012) with the RACMO data is used to determine regression parameters (slope, intercept) on a 5 km grid cell basis. Temperature data define melting degree days, which have a different coefficient for bare ice than snow cover, determined from hydrological-year cumulative SMB. A fundamental assumption is that the calibration factors, regression slope, and offset for the calibration period 1960 through 2012 are stationary over time, for which there is some evidence in Fettweis et al. (2017). Box et al. (2013) describe the methods in more detail.

Figure 3Comparison between This Study and the reconstruction (Kjeldsen et al.2015). All axis units are Gt yr−1. Plotted numbers represent the last two digits of the years for the unadjusted data sets. The matching colored squares show the adjusted data. MB* shown here does not include BMB for either the reconstructed or This Study data. Arrows show statistical properties before and after the adjustment. No adjustment is made to MB*, but it is computed from Eq. (4) both before (numbered) and after (squares) the surface and discharge adjustments.


The reconstructed surface mass balance is adjusted as described in the Methods Sect. 5.4 (Fig. 3).

4.2 Discharge

The recent discharge data are from Mankoff et al. (2020b) (data: Mankoff and Solgaard2020). This product covers all fast-flowing (>100m yr−1) marine-terminating glaciers. The discharge in Mankoff et al. (2020b) is computed at flux gates ∼5km upstream from glacier termini (Mankoff2020) using a wide range of velocity products and ice thickness from BedMachine v4. Discharge across flux gates is derived with a 200 m spatial resolution grid but then summed and provided at glacier resolution. Temporal coverage begins in 1986 with a few velocity estimates and is updated each time a new velocity product is released, which is every ∼12d with a ∼30d lag (Solgaard et al.2021; data: Solgaard and Kusk2021).

Some changes have been implemented since the last publication describing the discharge product (i.e., Mankoff et al.2020b). These are minor and include updating the Khan et al. (2016) (data: Khan2017) surface elevation change product from 2015 through 2019, updating various MEaSUREs velocity products to their latest version, updating the PROMICE Sentinel ice velocity product from Edition 1 ( to Edition 2 (Solgaard et al.2021; Solgaard and Kusk2021), and updating from BedMachine v3 (supplemented in the SE with Millan et al.2018) to use only BedMachine v4 (Morlighem et al.2021).

The reconstructed discharge data (Kjeldsen et al.2015) are estimated via a linear fit between unsmoothed annual discharge spanning 2000 to 2012 (Enderlin et al.2014) and runoff data from Kjeldsen et al. (2015) using a 6-year trailing average. The method for scaling discharge from runoff was introduced by Rignot et al. (2008), who scaled the SMB anomaly with discharge. Sensitivity analyses conducted by Box and Colgan (2013) showed runoff to be the more effective discharge predictor and include a discussion of the physical basis. Although the fitting period of the present data set includes an anomalous period of discharge (2000 through 2005; e.g., Boers and Rypdal2021), the discharge data used by Rignot et al. (2008) and Box and Colgan (2013) also include years 1958 and 1964 that lie near the regression line (see Box and Colgan2013, Fig. 4, and the related Sect. 4, Physical basis). Further, while 2000 through 2005 cover a changing period in Greenlandic discharge (Mankoff et al.2020b; King et al.2020), there were likely other anomalous periods in the past, when glaciers in Greenland experienced considerable increases in discharge, as inferred by geological and geodetic investigations (Andresen et al.2012; Bjørk et al.2012; Khan et al.2015, 2020).

The reconstructed discharge is adjusted as described in the Methods Sect. 5.4.

4.3 Basal mass balance

The BMB (Karlsson et al.2021) comes from mass lost at the bed from BMBGF, BMBfriction from the basal shear velocity, and BMBVHD from surface runoff routed to the bed (i.e., the volume of the subglacial conduits formed from surface runoff; Mankoff and Tulaczyk2017). These fields (data: Karlsson2021) are provided as steady-state annual estimates. We use the BMBGF and BMBfriction products and apply 1/365th to each day, each year. Because BMBVHD is proportional to runoff, an annual estimate is not appropriate for this work with daily resolution. We therefore re-calculate the BMBVHD-induced basal melt as described in the Methods Sect. 5.3.

4.3.1 Geothermal flux

Due to a lack of direct observations, the geothermal flux is poorly constrained under most of the Greenland ice sheet. Different approaches have been employed to infer the value of the BMBGF, often with diverging results (see, e.g., Rogozhina et al.2012; Rezvanbehbahani et al.2019). Lacking substantial validation that favors one BMBGF map over the others, Karlsson et al. (2021) instead use the average of three widely used BMBGF estimates: Fox Maule et al. (2009), Shapiro and Ritzwoller (2004), and Martos et al. (2018). The BMBGF melt rate is calculated as

(1) b ˙ m = E GF ρ i - 1 L - 1 ,

where EGF is available energy at the bed, here the geothermal flux in unit W m−2, ρi is the density of ice (917 kg m−3), and L is the latent heat of fusion (335 kJ kg−1; Cuffey and Paterson2010). BMBGF melting is only calculated where the bed is not frozen. We use the MacGregor et al. (2016) estimate of temperate bed extent and scale Eq. (1) by 0, 0.5, or 1 where the bed is frozen (∼25 % of the ice sheet area), uncertain (∼33 %), or thawed (∼42 %), respectively.

4.3.2 Friction

This heat term stems from the friction produced as ice slides over the bedrock. The term has only been measured in a handful of places (e.g., Ryser et al.2014; Maier et al.2019), and it is unclear how representative those measurements are at ice sheet scales. Karlsson et al. (2021) therefore estimate the frictional heating using the full Stokes Elmer/Ice model that resolves all stresses while relating basal sliding and shear stress using a linear friction law (Gillet-Chaulet et al.2012; Maier et al.2021). The model is tuned to match a multi-decadal surface velocity map (Joughin et al.2018) covering 1995–2015, and it returns an estimated basal friction heat that is used to calculate the basal melt due to friction, similarly to Eq. (1):

(2) b ˙ m = E f ρ i - 1 L - 1 ,

where Ef is energy due to friction. We also apply the 0, 0.5, and 1 scale as used for the BMBGF term (MacGregor et al.2016) in order to mask out areas that are likely frozen.

4.4 Other

ROI regions come from Mouginot and Rignot (2019) and ROI sectors come from Zwally et al. (2012).

4.5 Products used for validation

We validate This Study against five other data products (see Table 2 and Sect. 6). These products are the most recent IO product (Mouginot et al.2019), the previous PROMICE mass balance product (Colgan et al.2019; data: Colgan2021), the two mostly independent methods of estimating ice sheet mass change, GMB (Barletta et al.2013; data: Barletta et al.2020) and VC (Simonsen et al.2021a; data: Simonsen et al.2021b), and the IMBIE2 data (IMBIE Team2019). In addition to this, we evaluate the reconstructed Kjeldsen et al. (2015) (data: Box et al.2021) and This Study data during the overlapping period 1986 through 2012.

Mouginot et al. (2019)Colgan et al. (2019)Kjeldsen et al. (2015)

Table 2Summary of correlation, bias, and RMSE between different products during their overlap periods with This Study. Basal mass balance not included in This Study when comparing against Mouginot and Rignot (2019), Colgan et al. (2019), or Kjeldsen et al. (2015). Peripheral ice masses never included in This Study.

Download Print Version | Download XLSX

5 Methods

The total mass balance for all of Greenland and all the different ROIs involves summing each field (SMB, D, BMB) by each ROI and then subtracting the D and BMB from the SMB fields, or

(3) MB = SMB - D - BMB .

Products that do not include the BMB term (i.e., Mouginot et al.2019, Colgan et al.2019, and Kjeldsen et al.2015) have total mass balance defined as

(4) MB * = SMB - D ,

and when comparing This Study to those products, we compare like terms, never comparing our MB to a different product MB*, except in Fig. 4, where all products are shown together.

Figure 4Comparison between This Study and other mass balance time series. Note that various products do or do not include basal mass balance or peripheral ice masses (see Table 2). This Study annual-resolution data prior to 1986 are the Kjeldsen et al. (2015) data adjusted as described in Sect. 5.4. Sea level rise calculated as -Gt/361.8. Inset highlights changes since 2010. Data product version 74 from 25 October 2021 used to generate this graphic.


Prior to calculating the mass balance, we perform the following steps.

5.1 Surface mass balance

In This Study we generate an output based on each of the three RCMs (HIRHAM/HARMONIE, MAR, and RACMO); however, in addition to these we generate a final and fourth SMB field defined as a combination of (1) the adjusted reconstructed SMB from 1840 through 1985 (Sect. 5.4) and (2) the average of HIRHAM/HARMONIE, MAR, and RACMO from 1986 through a few months ago, the average of HIRHAM/HARMONIE and MAR from a few months ago through yesterday, and MAR from yesterday through next week. See Appendix A for differences between This Study MB and MB derived using each of the RCM SMBs. There is no obvious change or step function at the 1985 to 1986 reconstructed-to-recent change nor as the RACMO and then HIRHAM/HARMONIE RCMs become unavailable a few months ago and yesterday, respectively.

5.2 Projected discharge

We project the discharge from the last observed point from Mankoff et al. (2020b) (generally between 2 weeks and 1 month old) to 7 d into the future at each glacier. We define the long-term trend as the linear least squares fit to the last 3 years of data. The residual is the data minus the long-term trend. We define the seasonal signal as the daily average from each year of the last 3 years of the residual during the temporal window of interest that spans from the most recently available observation through next week. We shift the seasonal signal so that it is 0 on the first projected day. We then assign the value of the last observation, plus the long-term trend, plus the seasonal signal to the recent past-projected and future-forecasted D.

Discharge does not change sign and changes magnitude by approximately 6 % annually over the entire ice sheet (King et al.2018), but surface mass balance changes sign and has both larger and higher frequency variability. From this, the statistical forecast for discharge described above does not impact results as much as the physically based model forecast for surface mass balance.

5.3 Basal mass balance

Because Karlsson et al. (2021) provide a steady-state annual-average estimate of the BMB fields, we divide the BMBGF and BMBfriction fields by 365 to estimate daily average. This is a reasonable treatment of the BMBGF field, which does not have an annual cycle. The BMBfriction field does have a small annual cycle that matches the annual velocity cycle. However, when averaged over all of Greenland, this is only a ∼6 % variation (King et al.2018), and Karlsson et al. (2021) found that basal melt rates are 5 % higher during the summer. Thus, the intra-annual changes are less than the uncertainty. The BMBVHD field varies significantly throughout the year, because it is proportional to surface runoff. We therefore generate our own BMBVHD for This Study.

To estimate recent BMBVHD, we use daily MAR runoff (see Mankoff et al.2020a) and BedMachine v4 (Morlighem et al.2017, 2021) to derive subglacial routing pathways, similarly to Mankoff and Tulaczyk (2017). We assume that all runoff travels to the bed within the grid cell where it is generated, the bed is pressurized by the load of the overhead ice, and the runoff discharges on the day it is generated. We calculate subglacial routing from the gradient of the subglacial pressure head surface, h, defined as

(5) h = z b + k ρ i ρ w ( z s - z b ) ,

with zb the basal topography, k the flotation fraction (1), ρi the density of ice (917 kg m−3), ρw the density of water (1000 kg m−3), and zs the ice surface. Equation (5) comes from Shreve (1972), where the hydropotential has units of Pascals (Pa), but here it is divided by gravitational acceleration g times the density of water ρw to convert the units from Pascals to meters (Pa to m).

We compute h and from h streams and outlets and both the pressure and elevation difference between the source and outlet. The energy available for basal melting is the elevation difference (gravitational potential energy) and two-thirds of the pressure difference, with the remaining one-third consumed to warm the water to match the changing phase transition temperature (Liestøl1956; Mankoff and Tulaczyk2017). We assume all energy, EVHD (in Joules), is used to melt ice with

(6) b m = E VHD ρ i - 1 L - 1 .

Because results are presented per ROI, and to reduce the computational load of this daily estimate, we only calculate the integrated energy released between the RCM runoff source cell and the outlet cell and then assign that to the ROI containing the runoff source cell.

To estimate reconstructed basal mass balance, we treat BMBGF and BMBfriction as steady state as described at the start of this section. For BMBVHD we use the fact that VHD comes from runoff by definition, and from this, reconstructed BMBVHD is calculated using scaled runoff as a proxy. VHD theory suggests that a unit volume of runoff that experiences a 1000 m elevation drop will release enough heat to melt an additional 3 % (Liestøl1956). To estimate the scale factor, we use the 1986 through 2012 overlap between the Kjeldsen et al. (2015) runoff and This Study recent BMBVHD from MAR runoff described above. The correlation between the two has an r2 value of 0.75, a slope of 0.03, and an intercept of 3 Gt yr−1 (Appendix D). From this, we scale the Kjeldsen et al. (2015) reconstructed runoff by 3 % (from the 0.03 slope, unrelated to the theoretical 1000 m drop described earlier) to estimate reconstructed BMBVHD.

5.4 Reconstructed adjustment

We use the reconstructed and recent SMB and D overlap from 1986 through 2012 to adjust the reconstructed data. This Study vs. reconstructed SMB has a slope of 0.6 and an intercept of 166 Gt yr−1 (Fig. 3 SMB), and This Study vs. reconstructed D has a slope of 1.1 and an intercept of 17 Gt yr−1 (Fig. 3 D). The unadjusted reconstructed data slightly underestimate years with high SMB and overestimate years with low SMB (see 1986, 2010, 2011, and 2012 in Fig. 3 SMB). The unadjusted reconstructed data slightly overestimate years with low D and overestimate years with high D.

We adjust the reconstructed data until the reconstructed vs. recent slope is 1 and the intercept is 0 Gt yr−1 for each of the surface mass balance and discharge comparisons (Fig. 3). We then derive the BMBVHD term for reconstructed basal mass balance (Sect. 5.3 and Appendix D), bring in the other BMB terms (Sect. 5.3), and use Eq. (3) to compute the adjusted reconstructed mass balance.

For reconstructed SMB and D, the mean of the recent uncertainty is added to the reconstructed uncertainty during the adjustment. Reconstructed MB uncertainty is then re-calculated as the square root of the sum of the squares of the reconstructed SMB and D uncertainty.

For surface mass balance, the adjustment is effectively a rotation around the mean values, with years with low SMB decreasing and years with high SMB increasing after the adjustment. For discharge, years with low D are slightly reduced, and years with high D have a higher reduction to better match the overlapping estimates.

The adjustment described above treats all biases in the reconstructed data. The primary assumption of our adjustment is that the bias contributions do not change in proportion to each other over time. We attribute the disagreement and need for the adjustment to the demonstrated too-high biases in accumulation and ablation estimates in the 1840–2012 reconstructed SMB field (Fettweis et al.2020), an offset resulting from differences in ice masks (Kjeldsen et al.2015), the inclusion of peripheral glaciers (Kjeldsen et al.2015), other accumulation rate inaccuracies (Lewis et al.2017, 2019), and other unknowns.

5.5 Domains, boundaries, and regions of interest

Few of the ice masks used here are spatially aligned. The Zwally et al. (2012) sectors and the Mouginot and Rignot (2019) regions are often smaller than the RCM ice domains. For example, the RACMO ice domain is 1 718 959 km2, of which 1 696 419 km2 (99 %) are covered by the Mouginot and Rignot (2019) regions and 22 540 km2 (1 %) are not, or 1 678 864 km2 (98 %) are covered by Zwally et al. (2012) and 40 095 km2 (2 %) are not.

Cropping the RCM domain edges would remove the edge cells where the largest SMB losses occur. This effect is minor when SMB is high (years with low runoff, assuming SMB magnitude is dominated by the runoff term). This effect is large when SMB is low (years with high runoff). As an example of the 2010 decade, RACMO SMB has a mean of 251 Gt yr−1 for the decade, with a low of 45 Gt in 2019 and a high of 420 Gt in 2018. For these same extreme years RACMO cropped to Mouginot and Rignot (2019) has a low of 76 Gt (68 % high) and a high of 429 Gt (2 % high). RACMO cropped to Zwally et al. (2012) has a low of 84 Gt (85 % high) and a high of 429 Gt (2 % high).

We therefore grow the ROIs to cover the RCM domains. ROIs are grown by expanding them outward, assigning the new cells the value (ROI classification, that is, sector number or region name; see Fig. 1) of the nearest non-null cell, and then clipping to the RCM ice domain. This is done for each ROI and RCM. Appendix E provides a graphical display of the HIRHAM RCM domain, the Mouginot and Rignot (2019) domain, and our expanded Mouginot and Rignot (2019) domain.

BMBVHD comes from the MAR ice domain runoff but is generated on the BedMachine ice thickness grid, which is smaller than the ice domain in some places. Therefore, the largest runoff volumes per unit area (from the low-elevation edge of the ice sheet) are discarded in these locations.

6 Product evaluation and assessment

We compare to six related data sets (see Table 2 and Sect. 4.5): the most similar and recent IO product (Mouginot et al.2019), the previous PROMICE assessment (Colgan et al.2019), the two mostly independent methods (GMB, Barletta et al.2013, and VC, Simonsen et al.2021a), IMBIE2 (IMBIE Team2019), and the unadjusted reconstructed/recent overlap (Kjeldsen et al.2015).

Our initial comparison (Fig. 4) shows all seven products overlaid in a time series accumulating at the product resolution (daily to annual) from the beginning of the first overlap (1972, Mouginot et al.2019) until 7 d from now (now defined as 25 October 2021 based on the date this document is compiled). Each data set is manually aligned vertically so that the last timestamps appear to overlap, allowing disagreements to grow back in time. We also assume errors are smallest at present and allow errors to grow back in time. The errors for this product are described in the Uncertainty section.

In the sections below, we compare This Study to each of the validation data in more detail. The Mouginot et al. (2019) and Colgan et al. (2019) products allow term-level (SMB, D, and MB*) comparison and the GMB, VC, and IMBIE2 only MB-level comparison. The MB or MB* comparison for each product is summarized in Table 2. All have different masks. Bias (Gt yr−1) is defined as 1ni=1n(xi-yi). Root mean square error (RMSE) (Gt yr−1) is defined as 1ni=1n(xi-yi)2. Sums are computed using ice-sheet-wide annual values, where x is This Study, y is the other product, and a positive bias means that This Study has a larger value.

6.1 Mouginot (2019)

The Mouginot et al. (2019) product spans the 1972 through 2018 period. We only use 1986 and onward because This Study has annual resolution prior to 1986 and Mouginot et al. (2019) data are provided on a non-calendar-year period. The SMB comes from RACMO v2.3p2 downscaled at 1 km and agrees very well with SMB from This Study (r2 0.94, bias 11, RMSE 38, slope 1.1). The minor SMB differences are likely due to mask differences or our use of a three-RCM average SMB estimate.

Mouginot et al. (2019) discharge and our D from Mankoff et al. (2020b) have a 33 Gt yr−1 bias. This difference can mainly be attributed to different discharge estimates in the southeastern and central eastern sectors (Appendix: Mouginot regions). When we include BMB in This Study (diamonds in middle panel of Fig. 5 shifting values to the right), it adds ∼25Gt yr−1 to This Study.

Because MB* is a linear combination of SMB and D terms (Eq. 4), the MB* differences between this product and Mouginot et al. (2019) are dominated by the D term, although it is not apparent because interannual variability is dominated by SMB.

6.2 Colgan (2019)

The Colgan et al. (2019) product spans 1995 through 2015. The SMB term is broadly similar to the RCM-averaged SMB term in This Study, although Colgan et al. (2019) use only an older version of MAR (Fig. 6 top panel). The Colgan et al. (2019) SMB is spatially interpolated over the PROMICE ice sheet ice mask (Citterio and Ahlstrøm2013), which contains more detail on the ice sheet periphery and therefore a larger ablation area than the native coarser MAR ice mask. This Study does not interpolate the SMB field and instead works on the SMB ice domain.

Figure 5Comparison of This Study vs. Mouginot et al. (2019). All axis units are Gt yr−1. Plotted numbers represent the last two digits of the year. Matching colored diamonds show the data when BMB is added to This Study. Printed numbers (r2, bias, RMSE, slope) compare values without BMB.


Figure 6Comparison of This Study vs. Colgan et al. (2019). All axis units are Gt yr−1. Plotted numbers represent the last two digits of the year. Matching colored diamonds show the data when BMB is added to This Study. Printed numbers (r2, bias, RMSE, slope) compare values without BMB.


The largest difference between This Study and Colgan et al. (2019) is that the latter estimate grounding line ice discharge based on corrections to ice volume flow rate measured across the ∼1700m elevation contour. This is far inland relative to the grounding line flux gates used in This Study (from Mankoff2020). This introduces uncertainty into the Colgan et al. (2019) D term from SMB corrections between the 1700 m elevation contour and the terminus (see the large disagreement in Fig. 6 middle panel). This disagreement increases when BMB is included in the results of This Study (shown by the annual values shifting to the right).

The D disagreement is represented differently across sectors (Appendix: Colgan 2019), where sectors 1, 2, 5, and 6 all have correlation coefficients less than ∼0.1, while the remaining sectors 3, 4, 7, and 8 all have correlation coefficients greater than 0.5.

This Study assesses greater D bias (43 Gt yr−1) than Colgan et al. (2019). While Colgan et al. (2019) did not assess BMB, the majority of this discrepancy likely results from Colgan et al. (2019) aliasing the aforementioned downstream correction terms. For example, while This Study shows very little interannual variability in ice discharge in the predominantly land-terminating SW region, Colgan et al. (2019) infer large interannual variability in ice discharge based on large interannual variability in SMB and changes in ablation area ice volume in their Sector 6. The discrepancy between This Study and the Colgan et al. (2019) D[+BMB] is largest during the earliest part of the record (i.e., 1995–2000), decreasing towards the present day, which may suggest that Colgan et al. (2019) particularly overestimated the response in ice discharge to 1990s climate variability.

Similarly to the comparison with Mouginot et al. (2019), the disagreement between This Study and Colgan et al. (2019) is dominated by D disagreement, although it is again not apparent because interannual variability is dominated by SMB.

6.3 GMB

Unlike This Study, the GMB method includes mass losses and gains on peripheral ice masses which should introduce a bias of ∼10 % to 15 % (Colgan et al.2015; Bolch et al.2013). The inclusion of peripheral ice in the GMB product is because the spatial resolution is so low that it cannot distinguish between them and the main ice sheet. There is also signal leakage from other glaciated areas, e.g., the Canadian Arctic. This can have an effect on the estimated signal, especially in sectors 1 and 8 or regions NW and NO. There is also leakage between basins, which becomes a larger issue for smaller basins or where major outlet glaciers are near basin boundaries. GMB may also have an amplified seasonal signal due to changing snow loading in the surrounding land areas that may be mapped as ice sheet mass change variability. This would enhance the seasonal amplitude but not have an impact on the interannual mass change rates. Additionally, different glacial isostatic adjustment (GIA) corrections applied to the gravimetric signal may also lead to differences in GMB estimates on an ice sheet scale but also on a sector scale (e.g., Sutterley et al.2014; Khan et al.2016).

GMB and the IO method (This Study) both report changes in ice sheet mass, but they are measuring two fundamentally different things. The IO method tracks volume flow rate across the ice sheet boundaries. Typically this is meltwater across the ice sheet surface and solid ice across flux gates near the calving edge of the ice sheet, and in This Study also meltwater across the ice sheet basal boundary. That volume is then converted to mass. We consider that mass is “lost” as soon as it crosses the boundary (i.e., the ice melts or ice crosses the flux gate). The GMB method tracks the regional mass changes. Melting ice has no impact on this until the meltwater enters the ocean and a similar mass leaves the far-field GMB footprint. From these differences, the GMB method may be a better estimate of sea level rise, while the IO method may be a better representation of the state of the Greenland ice sheet.

6.4 VC

When deriving surface elevation change from satellite altimetry, data from multiple years are needed to give a stable ice-sheet-wide prediction. Hence, the altimetric mass balance estimates are often reported as averages of single satellite missions.

Although This Study has a small (−𝟷𝟷Gt yr−1) bias in comparison to Simonsen et al. (2021a) VC, there is a relatively high RMSE of 86Gt yr−1 and a mid-range correlation (r2=0.62). This suggests that while both This Study and VC agree on the total mass loss of the ice sheet, they disagree on the precise temporal distribution of this mass loss. It is possible that the outlying 1992 and 2019 years are influenced by the edge of the time series record if not fully sampled, but other outliers exist – the 1992 extreme low melt year and the 2019 extreme melt year as well as the 1995 through 1998 period stand out as years with poor agreement.

Figure 7This Study total mass balance (MB) vs. the gravimetric method (GMB), volume change method (VC), and IMBIE2 estimates of MB. All three include BMB. All axis units are Gt yr−1. Plotted numbers represent the last two digits of the year. GRACE and IMBIE2 include peripheral ice masses.


We suggest that this is due to climate influences on the effective radar horizon across the ice sheet during these years. Weather-driven change in the effective scatter horizon, mapped by the Ku band in the upper snow layer of ice sheets, hampers the conversion of radar-derived elevation change into mass change (Nilsson et al.2015). Simonsen et al. (2021a) used a machine learning approach to derive a temporal calibration field for converting the radar elevation change estimates into mass change. This approach relied on precise mass balance estimates from ICESat to train the model and thereby was able to remove the effects of the changing scattering horizon in the radar data. This VC mass balance is given for monthly time steps (Simonsen et al.2021a); however, the running mean applied to derive radar elevation change will dampen the interannual variability of the mass balance estimate from VC. This is especially true prior to 2010, after which the novel radar altimeter onboard CryoSat-2 allowed for a shortening of the data windowing from 5 to 3 years. This smoothing of the interannual variability is also seen in the intercomparison between This Study and the VC MB, where in addition to the two end-members of the time series (1992 and 2019) the years 1995, 1996, and 1998 seem to be outliers (Fig. 7). These years are notable for high MB, which seems to be captured less precisely by the older radar altimeters due to the longer temporal averaging.


The most widely cited estimate of Greenland mass balance today is the Ice-Sheet Mass Balance Inter-Comparison Exercise 2 (IMBIE2, IMBIE Team2019). IMBIE2 seeks to provide a consensus estimate of monthly Greenland mass balance between 1992 and 2018 that is derived from altimetry, gravimetry, and input–output ensemble members. There are two critical methodological differences between This Study and IMBIE2. Firstly, the gravimetry members of IMBIE2 assess the mass balance of all Greenlandic land ice, including peripheral ice masses, while This Study only assesses the mass balance of the ice sheet proper. Secondly, the input–output members of IMBIE2 do not assess BMB, while This Study does.

The IMBIE2 composite record of ice sheet mass balance equally weights three methods of assessing ice sheet mass balance: input–output, altimetry, and gravimetry. Prior to ca. 2003, however, IMBIE2 is derived solely from IO studies that explicitly exclude BMB (MB is actually MB*). After ca. 2003, by comparison, IMBIE2 includes both satellite altimetry and gravimetry records implicitly sampling BMB. The representation of BMB in the composite IMBIE2 mass balance record therefore shifts before and after ca. 2003.

In comparison to mass balance assessed by IMBIE2, This Study has a small bias of ∼-7Gt yr−1 over the 26-calendar-year comparison period. This apparent agreement may be attributed to the compensating effects of IMBIE2 effectively sampling peripheral ice masses and ignoring BMB, while This Study does the opposite and ignores peripheral ice masses but samples BMB, equal to ∼25Gt yr−1. Over the entire 26-year comparison period, the RMSE with IMBIE2 is 44Gt yr−1 and the correlation is 0.89. This relatively high correlation highlights good agreement in interannual variability between studies, and the RMSE suggests that formal stated uncertainties of each study (ca. ±30 to ±63Gt yr−1 for IMBIE2 and a mean of 86 Gt yr−1 for This Study) are indeed good estimates of the true uncertainty, as assessed by inter-study discrepancies.

7 Uncertainty

We treat the three inputs to the total mass balance (surface mass balance, discharge, and basal mass balance, or SMB, D, and BMB) as independent when calculating the total error. This is a simplification – the RCM SMB and the BMBVHD from RCM runoff are related and D ice thickness and BMBVHD pressure gradients are related, and other terms may have dependencies. However, the two dominant IO terms, SMB inputs and D outputs, are independent on annual timescales, and for simplification we treat all terms as independent. We use Eq. (3) and standard error propagation for SMB, D, and BMB terms (i.e., the square root of the sum of the squares of the SMB plus D plus BMB error terms). For D, extra work is done to calculate uncertainty between the last Mankoff et al. (2020b) D data (up to 30 d old, with an error of ∼9 % or ∼45Gt yr−1) and the forecasted now-plus-7 d D (see Sect. 7.1). Table 3 provides a summary of the uncertainty for each input.

Langen et al. (2017)Fettweis et al. (2020)Noël et al. (2019)Kjeldsen et al. (2015)Mankoff et al. (2020b)(Mankoff2021)Kjeldsen et al. (2015)Karlsson et al. (2021)Karlsson et al. (2021)

Table 3Summary of uncertainty estimates for products used in This Study. This is an approximate and simplified representation – RCM uncertainties are calculated separately for gain and loss terms, because SMB near 0 does not mean uncertainty is near 0. This is also why the final This Study uncertainty is presented with units Gt yr−1.

Download Print Version | Download XLSX

The final This Study MB uncertainty value shown in Table 3 comes from the mean of the annual sum of the MB error term.

7.1 Discharge

The D uncertainty is discussed in detail in Mankoff et al. (2020b), but the main uncertainties come from unknown ice thickness, the assumption of no vertical shear at fast-flowing marine-terminating outlet glaciers, and ice density of 917 kg m−3. Regional ice density can be significantly reduced by crevasses. For example, Mankoff et al. (2020c) identified a snow-covered crevasse field with 20 % crevasse density, meaning at that location regional firn density should be reduced by 20 %.

Temporally, D at daily resolution comes from ∼12d observations upsampled to daily, and those ∼12d resolution observations come from longer time period observations (Solgaard et al.2021). Because the velocity method uses feature tracking, it is correct on average but misses variability within each sample period (e.g., Greene et al.2020).

Spatially, discharge is estimated ∼5km upstream from the grounding lines for ice velocities as low as 100 m yr−1. That ice accelerates toward the margin, but even ice flowing steadily at 1 km yr−1 would take 5 years before that mass is lost. However, at any given point in time, ice that had previously crossed the flux gate is calving or melting into the fjord. The discrepancy here between the flux gate estimated mass loss and the actual mass lost at the downstream terminus is only significant for glaciers that have had large velocity changes at some point in the recent past, large changes in ice thickness, or large changes in the location (retreat or advance) of the terminus. We do not consider SMB changes downstream of the flux gate, because the gates are temporally near the terminus for most of the ice that is fast-flowing, and the largest SMB uncertainty is at the ice sheet margin, where there are both mask issues and high topographic variability.

The forecasted D uncertainty is the average historical uncertainty plus a 1 % increase per day for the past projected and forecasted period.

7.2 ROIs

We work on the three different domains of the three RCMs and expand the ROIs to match the RCMs (see Appendix E). However, some alignment issues cannot be solved. For example, we use BedMachine ice thickness to estimate BMBVHD. Often, the largest BMBVHD occurs near the ice margin under ice with the steepest surface slopes. This is also where the largest runoff often occurs, because the ice margin, at the lowest elevations, is exposed to the warmest air. If these RCM ice grid cells with high runoff are anywhere inside the BedMachine ice domain, that runoff is still included in our BMBVHD estimates because it flows outward and passes through the BedMachine near-ice-edge grid cells with the large pressure gradients. However, any RCM ice runoff outside the BedMachine ice domain (ice thickness is 0) is ignored.

The MAR ice domain is 1 825 600 km2, of which 1 708 400 km2 (94 %) are covered by the BedMachine ice mask and 26 400 km2 (6 %) are not. This 6 % area contributes ∼18 % of runoff on average (range of 16 % to 21 % from 2010 through 2019). This 18 % of runoff is excluded from the VHD calculations and likely contributes more than 18 % to the VHD term, because the border region of the ice sheet has the steepest gradients and the largest volume of subglacial flow.

We encourage RCM developers, BedMachine, and others to use a common and up-to-date mask (see Kjeldsen et al.2020).

7.3 Accumulating uncertainties

When accumulating errors as in Fig. 4, we use only the D and BMBGF uncertainty. The D uncertainty is primarily due to unknown ice thickness and is invariant in time, and the geothermal heat flux is steady state. SMB uncertainty is assumed to have errors randomly distributed in time (for the purposes of Fig. 4). There may be time-invariant biases in the BMBfriction and SMB fields, but treating all uncertainties as biases is incorrect – evidence for that comes from the six other MB estimates. This distinction between bias and random uncertainty is only done for Fig. 4, where errors accumulate in time. The provided data product contains one uncertainty field and does not distinguish between systematic and random uncertainty. We caution others in treating SMB uncertainty as random in time for analyses that go beyond the graphical display used here.

The shaded region in Fig. 4 representing the uncertainty for This Study is computed as 365 d rolling smoothly from 1840 through 1999 of the above-described uncertainty, 1/365th of the annual error at now +7d, and a linear blend, from 2000 to now +7d, between the smoothed reconstructed uncertainty and the present and future more variable uncertainty.

The Mouginot et al. (2019), Colgan et al. (2019), and Kjeldsen et al. (2015) products all provide an error estimate but do not distinguish between temporally fixed errors (biases; should accumulate in time) vs. temporally random errors.

We treat the Mouginot et al. (2019) data the same as This Study. Discharge uncertainty is treated as a bias and accumulates, and surface mass balance uncertainty is treated as random and does not accumulate.

The Colgan et al. (2019) vs. This Study bias and RMSE are 32 and 59 Gt yr−1, respectively. This suggests that in any given year, there could be up to -32±59 or +27/-91Gt yr−1 departure from This Study. From this, we assign a 32 Gt yr−1 bias (35 %; accumulates in time) and a 59 Gt yr−1 RMSE (65 %; random in time).

The adjusted Kjeldsen et al. (2015) data have 0 surface mass balance and discharge bias by definition (Sect. 5.4), but Fig. 4 displays the unadjusted data, and we apply a 36 Gt yr−1 accumulating uncertainty from the unadjusted D bias (Fig. 3).

7.4 Peripheral ice masses

Greenland's peripheral glaciers and ice caps are not included in this product. Nonetheless, we briefly summarize recent mass balance estimates of these areas. Greenlandic peripheral ice contributes more runoff per unit area than the main ice sheet – it is <5 % of the total ice area but contributes ∼15 % to 20 % of the whole island mass loss (Bolch et al.2013). From 2003 to 2009 and using the VC method (altimetry), Gardner et al. (2013) estimate -38±7Gt yr−1 peripheral mass balance. From 2006 to 2016 and using the VC method (DEM differencing), Zemp et al. (2020) estimate -51±17Gt yr−1 peripheral mass balance using Rastner et al. (2012) delineations.

8 Results

From the 181 complete years of data (excluding partial 2021), the mean mass balance is -77±125Gt yr−1, with a minimum of -428±110Gt in 2012 (SMB of 87±8Gt, D of 485±46Gt, BMB of 29±6Gt) and a maximum of 142±83Gt yr−1 in 1996 (SMB of 584±53Gt, D of 420±39Gt, BMB of 21±5Gt).

At the decadal average, the following trends are apparent. Surface mass balance has decreased from a high of ∼450Gt yr−1 in the 1860s to a low of ∼260Gt yr−1 in the 2010s. SMB variability has also increased during this time. Discharge has increased slightly from a low of ∼375Gt yr−1 in the 1860s to a high of ∼490Gt yr−1 in the 2010s. Basal mass balance, from runoff as a proxy, had a high of 26±16Gt yr−1 in the 1930s and a low of 22±5Gt yr−1 in the 1990s but, as with runoff, has been increasing in recent decades.

The total mass balance decadal trend from the 1840s through the 2010s is one of general mass decrease and increased intra-decadal variability. The record begins in the 1840s with -10Gt yr−1, has only 1 (of 19) decades with a mass gain (∼50Gt yr−1 in the 1860s), and has a record low of -250Gt yr−1 in the 2010s.

9 Data availability

The RCM surface mass balance and the VHD basal mass balance components are updated daily, the discharge approximately every 12 d, and all are used to produce the final daily-updating product. The data area is available at (Mankoff et al.2021), with all historical (daily updated) versions archived.

As part of our commitment to making continual and improving updates to the data product, we introduce a GitHub database (, last access: 26 October 2021, Mankoff2021) where users can track progress, make suggestions, and discuss, report, and respond to issues that arise during use of this product.

10 Conclusions

This study is the first to provide a data set containing more than a century and real-time estimates detailing the state of Greenland ice sheet mass balance, with regional or sector spatial and daily temporal resolution products of surface mass balance, discharge, basal mass balance, and the total mass balance.

IMBIE2 highlights that during the GRACE satellite gravimetry era (2003 through 2017), there are usually more than 20 independent estimates of annual Greenland ice sheet mass balance. Just two independent estimates, however, are available prior to 2003. This study will therefore provide additional insight into ice sheet mass balance during the late 1980s and 1990s. IMBIE2 also highlights how the availability of mass balance estimates declines in the year prior to IMBIE2 publication. This reflects a lag period during which mass balance assessments from non-operational products are undergoing peer review. The operational nature of this product supports the timely inclusion of annual MB estimates in community consensus reports such as those from IMBIE and the IPCC.

As such, the data products provided in This Study present the first operational monitoring of the Greenland ice sheet total mass balance and its components. One property of the input–output approach used in This Study is the explanatory capabilities of the data products, allowing scrutiny of the physical origins of recorded mass changes. By excluding peripheral ice masses, This Study allows and invites anyone to keep an eye on the current evolution of the Greenland ice sheet proper. However, as the spatial resolutions of RCMs increase and estimates of peripheral ice thickness become available, our setup allows inclusion of these ice masses to generate a full Greenland-wide product. Moreover, as the determination of each of the individual components of the ice sheet mass balance is expected to improve over time through international research efforts, the total mass balance product presented will also be able to improve, as it is sustained by the Danish–Greenlandic governmental long-term monitoring effort – PROMICE.

Appendix A: RCM differences

Figure A1Comparison of This Study combined RCM product and the HIRHAM/HARMONIE, MAR, and RACMO RCMs. Results shown here are MB, not SMB, but the same D and BMB have been subtracted from each SMB product. (a) Annual MB for the entire time series. (b) Example 2 years (2019 and 2020) at daily resolution. (c) Difference between the three RCM MB products and This Study RCM-averaged product for the same data shown in panel (b).


Appendix B: Mouginot 2019 by region

Figure B1Comparison between This Study (excluding BMB) and Mouginot et al. (2019). Same data and display as Fig. 5 except here displayed by the Mouginot and Rignot (2019) region. Numbers in each graph show r2, bias, and RMSE from top to bottom, respectively. All axis units are Gt yr−1.


Appendix C: Colgan 2019 by sector

Figure C1Comparison between This Study (excluding BMB) and Colgan et al. (2019). Same data and display as Fig. 6 except here displayed by the Zwally et al. (2012) sector. Numbers in each graph show r2, bias, and RMSE from top to bottom, respectively. All axis units are Gt yr−1.


Appendix D: Reconstructed runoff

Figure D1Comparison between MAR runoff and basal viscous heat dissipation derived from that runoff. The slope is used to estimate the reconstructed BMBVHD from reconstructed runoff (see Sect. 5.3). Axis units are Gt yr−1. Plotted numbers represent the last two digits of the year.


Appendix E: RCM coverage

Figure E1HIRHAM RCM coverage by Mouginot and Rignot (2019). Coverage of HIRHAM by Zwally et al. (2012), and MAR and RACMO by Mouginot and Rignot (2019) and Zwally et al. (2012) are similar to the graphic shown here (see Sect. 5.5 for a discussion of RACMO coverage issues). HIRHAM latitude and longitude cover the Equator because we work on the native HIRHAM rotated pole coordinate system.

Appendix F: Software

This work was performed using only open-source software, primarily GRASS GIS (Neteler et al.2012), CDO (Schulzweida2019), NCO (Zender2008), GDAL (GDAL/OGR contributors2020), and Python (Van Rossum and Drake Jr1995), in particular the Jupyter (Kluyver et al.2016), dask (Dask Development Team2016; Rocklin2015), pandas (McKinney2010), numpy (Oliphant2006), x-array (Hoyer and Hamman2017), and Matplotlib (Hunter2007) packages. The entire work was performed in Emacs (Stallman1981) using Org Mode (Schulte et al.2012) on GNU/Linux and using many GNU utilities. The parallel (Tange2011) tool was used to speed up processing.

Appendix G: CRediT

Figure G1Author contributions following the CRediT system (Allen et al.2014; Brand et al.2015; Allen et al.2019).


Author contributions

Author contribution is captured following the CRediT system (Allen et al.2014; Brand et al.2015; Allen et al.2019) and shown graphically in Fig. G1. The following authors contributed in the following ways. Conceptualization from KDM, APA, and RSF. Curation from KDM, XF, PLL, MS, KKK, NBK, BN, MRvdB, AS, and JEB. Implementation from KDM, XF, PLL, KKK, and MDK. Funding from AS, APA, SBA, and RSF. SMB methods from XF, PLL, BN, and MvdB. D methods from KDM, WC, AS, MDK, APA, and RSF. BMB methods from NBK and KDM. General validation by KDM. GRACE validation from WC. VC validation by WC and SBS. Reconstruction methods from KKK, JEB, and KDM. Project administration by KDM, APA, SBA, and RSF. Resources from KDM, XF, PLL, MS, KKK, NBK, BN, MRvdB, AS, and SBA. Software written by KDM, XF, PLL, AS, and MDK. Visualization by KDM. Writing by KDM, XF, PLL, MS, KKK, NBK, BN, MRvdB, WC, JEB, SBS, APA, and RSF. Editing by KDM, XF, PLL, MS, KKK, NBK, BN, MRvdB, WC, JEB, SBS, MDK, and RSF.

Competing interests

The authors declare that they have no conflict of interest.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issues “Extreme environment datasets for the three poles” and “Towards an integrated Arctic observation system to fill gaps of observing system across the atmosphere, ocean, cryosphere, geosphere, and terrestrial ecosystems”. It is not associated with a conference.


The editor and anonymous reviewers (Anonymous2021a, b) provided feedback and helped improve this paper.

Financial support

This research has been supported by the Programme for Monitoring of the Greenland ice sheet (PROMICE). Parts of this work were funded by the INTAROS project under the European Union's Horizon 2020 research and innovation program under grant agreement no. 727890. Brice Noël was funded by NWO VENI grant VI.Veni.192.019. Michiel R. van den Broeke was funded by the Netherlands Earth System Science Centre (NESSC).

Review statement

This paper was edited by David Carlson and reviewed by two anonymous referees.


Allen, L., Scott, J., Brand, A., Hlava, M., and Altman, M.: Publishing: Credit where credit is due, Nature, 508, 312–313,, 2014. a, b

Allen, L., O'Connell, A., and Kiermer, V.: How can we ensure visibility and diversity in research contributions? How the Contributor Role Taxonomy (CRediT) is helping the shift from authorship to contributorship, Learn. Publ., 32, 71–74,, 2019. a, b

Amory, C., Kittel, C., Le Toumelin, L., Agosta, C., Delhasse, A., Favier, V., and Fettweis, X.: Performance of MAR (v3.11) in simulating the drifting-snow climate and surface mass balance of Adélie Land, East Antarctica, Geosci. Model Dev., 14, 3487–3510,, 2021. a

Andresen, C. S., Straneo, F., Ribergaard, M. H., Bjørk, A. A., Andersen, T. J., Kuijpers, A., Nørgaard-Pedersen, N., Kjær, K. H., Schjøth, F., Weckström, K., and Ahlstrøm, A. P.: Rapid response of Helheim Glacier in Greenland to climate variability over the past century, Nat. Geosci., 5, 37–41,, 2012. a

Anonymous: RC1: Comment on essd-2021-131,, 2021a. a

Anonymous: RC2: Comment on essd-2021-131,, 2021b. a

Barletta, V. R., Sørensen, L. S., and Forsberg, R.: Scatter of mass changes estimates at basin scale for Greenland and Antarctica, The Cryosphere, 7, 1411–1432,, 2013. a, b, c

Barletta, V. R., Sørensen, L. S., Simonsen, S. B., and Forsberg, R.: Gravitational mass balance of Greenland 2003 to present v2.1, DTU Data [data set],, 2020. a

Bengtsson, L., Andrae, U., Aspelien, T., Batrak, Y., Calvo, J., de Rooy, W., Gleeson, E., Hansen-Sass, B., Homleid, M., Hortal, M., Ivarsson, K.-I., Lenderink, G., Niemelä, S., Nielsen, K. P., Onvlee, J., Rontu, L., Samuelsson, P., Muñoz, D. S., Subias, A., Tijm, S., Toll, V., Yang, X., and Køltzow, M. Ø.: The HARMONIE–AROME Model Configuration in the ALADIN–HIRLAM NWP System, Mon. Weather Rev., 145, 1919–1935,, 2017. a

Bjørk, A. A., Kjær, K. H., Korsgaard, N. J., Khan, S. A., Kjeldsen, K. K., Andresen, C. S., Box, J. E., Larsen, N. K., and Funder, S.: An aerial view of 80 years of climate-related glacier fluctuations in southeast Greenland, Nat. Geosci., 5, 427–432,, 2012. a

Boers, N. and Rypdal, M.: Critical slowing down suggests that the western Greenland Ice Sheet is close to a tipping point, P. Natl. Acad. Sci. USA, 118, e2024192118,, 2021. a

Bolch, T., Sørensen, L. S., Simonsen, S. B., Mölg, N., Machguth, H., Rastner, P., and Paul, F.: Mass loss of Greenland's glaciers and ice caps 2003–2008 revealed from ICESat data, Geophys. Res. Lett., 40, 875–881,, 2013. a, b

Box, J., Colgan, W. T., and Kjeldsen, K.: Greenland SMB, D and TMB annual time series 1840–2012, GEUS Dataverse [data set],, 2021. a, b, c, d, e, f

Box, J. E.: Greenland ice sheet mass balance reconstruction. Part II: surface mass balance (1840–2010), J. Climate, 26, 6974–6989,, 2013. a, b

Box, J. E. and Colgan, W.: Greenland ice sheet mass balance reconstruction. Part III: marine ice loss and total mass balance (1840–2010), J. Climate, 26, 6990–7002,, 2013. a, b, c

Box, J. E., Cressie, N., Bromwich, D. H., Jung, J.-H., van den Broeke, M. R., van Angelen, J. H., Forster, R. R., Miège, C., Mosley-Thompson, E., Vinther, B., and McConnell, J. R.: Greenland ice sheet mass balance reconstruction. Part I: Net snow accumulation (1600–2009), J. Climate, 26, 3919–3934,, 2013. a, b

Brand, A., Allen, L., Altman, M., Hlava, M., and Scott, J.: Beyond authorship: attribution, contribution, collaboration, and credit, Learn. Publ., 28, 151–155,, 2015. a, b

Cappelen, J.: The observed climate of Greenland, 1958–99-with climatological standard normals, 1961–90, Danish Meteorological Institute, Copenhagen, Denmark, 2001. a

Cappelen, J., Laursen, E. V., Jørgensen, P. V., and Kern-Hansen, C.: DMI Monthly Climate Data Collection 1768–2005: Denmark, The Faroe Islands and Greenland, DMI, Copenhagen, Denmark, 2006. a

Cappelen, J., Laursen, E., Jørgensen, P., and Kern-Hansen, C.: DMI monthly climate data collection 1768–2010, Denmark, the Faroe Islands and Greenland, Tech. rep., DMI Technical Report 11-05, Copenhagen, 2011. a

Christensen, O. B., Drews, M., Hesselbjerg Christensen, J., Dethloff, K., Ketelsen, K., Hebestadt, I., and Rinke, A.: The HIRHAM Regional Climate Model. Version 5 (beta), no. 06-17 in Denmark, Danish Meteorological Institute, Technical Report, Danish Climate Centre, Danish Meteorological Institute, Copenhagen, Denmark, 2007. a

Citterio, M. and Ahlstrøm, A. P.: Brief communication “The aerophotogrammetric map of Greenland ice masses”, The Cryosphere, 7, 445–449,, 2013. a, b

Colgan, W.: Greenland ice sheet mass balance assessment (1995–2019), GEUS Dataverse [data set],, 2021. a

Colgan, W., Abdalati, W., Citterio, M., Csatho, B., Fettweis, X., Luthcke, S., Moholdt, G., Simonsen, S. B., and Stober, M.: Hybrid glacier Inventory, Gravimetry and Altimetry (HIGA) mass balance product for Greenland and the Canadian Arctic, Remote Sens. Environ., 168, 24–39,, 2015. a

Colgan, W., Mankoff, K. D., Kjeldsen, K. K., Bjørk, A. A., Box, J. E., Simonsen, S. B., Sørensen, L. S., Khan, S. A., Solgaard, A. M., Forsberg, R., Skourup, H., Stenseng, L., Kristensen, S. S., Hvidegaard, S. M., Citterio, M., Karlsson, N., Fettweis, X., Ahlstrøm, A. P., Andersen, S. B., van As, D., and Fausto, R. S.: Greenland ice sheet mass balance assessed by PROMICE (1995–2015), Geol. Surv. Den. Greenl., 43, e2019430201,, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x

Cuffey, K. M. and Paterson, W. S. B.: The physics of glaciers, 4th edn., Academic Press, Burlington, MA USA, 2010. a

Dask Development Team: Dask: Library for dynamic task scheduling, available at: (last access: 26 October 2021), 2016. a

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, I., Biblot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Greer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Holm, E. V., Isaksen, L., Kallberg, P., Kohler, M., Matricardi, M., McNally, A. P., Mong-Sanz, B. M., Morcette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thepaut, J. N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597,, 2011. a, b

Eerola, K.: About the performance of HIRLAM version 7.0, Hirlam Newsletter, 51, 93–102, 2006. a

Enderlin, E. M., Howat, I. M., Jeong, S., Noh, M.-J., van Angelen, J. H., and van den Broeke, M. R.: An improved mass budget for the Greenland Ice Sheet, Geophys. Res. Lett., 41, 866–872,, 2014. a

Ettema, J., van den Broeke, M. R., van Meijgaard, E., van de Berg, W. J., Box, J. E., and Steffen, K.: Climate of the Greenland ice sheet using a high-resolution climate model – Part 1: Evaluation, The Cryosphere, 4, 511–527,, 2010. a, b

Fettweis, X., Box, J. E., Agosta, C., Amory, C., Kittel, C., Lang, C., van As, D., Machguth, H., and Gallée, H.: Reconstructions of the 1900–2015 Greenland ice sheet surface mass balance using the regional climate MAR model, The Cryosphere, 11, 1015–1033,, 2017. a

Fettweis, X., Hofer, S., Krebs-Kanzow, U., Amory, C., Aoki, T., Berends, C. J., Born, A., Box, J. E., Delhasse, A., Fujita, K., Gierz, P., Goelzer, H., Hanna, E., Hashimoto, A., Huybrechts, P., Kapsch, M.-L., King, M. D., Kittel, C., Lang, C., Langen, P. L., Lenaerts, J. T. M., Liston, G. E., Lohmann, G., Mernild, S. H., Mikolajewicz, U., Modali, K., Mottram, R. H., Niwano, M., Noël, B., Ryan, J. C., Smith, A., Streffing, J., Tedesco, M., van de Berg, W. J., van den Broeke, M., van de Wal, R. S. W., van Kampenhout, L., Wilton, D., Wouters, B., Ziemen, F., and Zolles, T.: GrSMBMIP: intercomparison of the modelled 1980–2012 surface mass balance over the Greenland Ice Sheet, The Cryosphere, 14, 3935–3958,, 2020. a, b, c, d, e, f, g, h, i

Fox Maule, C., Purucker, M. E., and Olsen, N.: Inferring magnetic crustal thickness and geothermal heat flux from crustal magnetic field models, Danish Climate Centre Report, Copenhagen, Denmark, 2009. a

Gallée, H. and Schayes, G.: Development of a Three-Dimensional Meso-γ Primitive Equation Model: Katabatic Winds Simulation in the Area of Terra Nova Bay, Antarctica, Mon. Weather Rev., 122, 671–685,<0671:doatdm>;2, 1994. a

Gallée, H., Guyomarc'h, G., and Brun, E.: Impact Of Snow Drift On The Antarctic Ice Sheet Surface Mass Balance: Possible Sensitivity To Snow-Surface Properties, Bound.-Lay. Meteorol., 99, 1–19,, 2001. a

Gardner, A. S., Moholdt, G., Cogley, J. G., Wouters, B., Arendt, A. A., Wahr, J., Berthier, E., Hock, R., Pfeffer, W. T., Kaser, G., Ligtenberg, S. R. M., Bolch, T., Sharp, M. J., Hagen, J. O., Van Den Broeke, M. R., and Paul, F.: A Reconciled Estimate of Glacier Contributions to Sea Level Rise: 2003 to 2009, Science, 340, 852–857,, 2013. a

GDAL/OGR contributors: GDAL/OGR Geospatial Data Abstraction software Library, Open Source Geospatial Foundation, available at: (last access: 26 October 2021), 2020. a

Gillet-Chaulet, F., Gagliardini, O., Seddik, H., Nodet, M., Durand, G., Ritz, C., Zwinger, T., Greve, R., and Vaughan, D. G.: Greenland ice sheet contribution to sea-level rise from a new-generation ice-sheet model, The Cryosphere, 6, 1561–1576,, 2012. a

Greene, C. A., Gardner, A. S., and Andrews, L. C.: Detecting seasonal ice dynamics in satellite images, The Cryosphere, 14, 4365–4378,, 2020. a

Groh, A., Horwath, M., Horvath, A., Meister, R., Sørensen, L. S., Barletta, V. R., Forsberg, R., Wouters, B., Ditmar, P., Ran, J., Klees, R., Su, X., Shang, K., Guo, J., Shum, C. K., Schrama, E., and Shepherd, A.: Evaluating GRACE Mass Change Time Series for the Antarctic and Greenland Ice Sheet – Methods and Results, Geosciences, 9, 415,, 2019. a

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., Chiara, G. D., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049,, 2020. a

Howat, I. M., Negrete, A., and Smith, B. E.: The Greenland Ice Mapping Project (GIMP) land classification and surface elevation data sets, The Cryosphere, 8, 1509–1518,, 2014. a, b

Hoyer, S. and Hamman, J. J.: xarray: N-D labeled Arrays and Datasets in Python, Journal of Open Research Software, 5,, 2017. a

Hunter, J. D.: Matplotlib: A 2D graphics environment, Comput. Sci. Eng., 9, 90–95, 2007. a

Joughin, I., Smith, B. E., and Howat, I.: Greenland Ice Mapping Project: ice flow velocity variation at sub-monthly to decadal timescales, The Cryosphere, 12, 2211–2227,, 2018. a

Karlsson, N. B.: Greenland Ice Sheet Basal Melt, GEUS Dataverse [data set],, 2021. a, b, c

Karlsson, N. B., Solgaard, A. M., Mankoff, K. D., Gillet-Chaulet, F., MacGregor, J. A., Box, J. E., Citterio, M., Colgan, W. T., Larsen, S. H., Kjeldsen, K. K., Korsgaard, N. J., Benn, D. I., Hewitt, I., and Fausto, R. S.: A first constraint on basal melt-water production of the Greenland ice sheet, Nat. Commun., 12, 3461,, 2021. a, b, c, d, e, f, g, h, i, j

Khan, S. A.: Greenland Ice Sheet Surface Elevation Change, GEUS Data Center [data set],, 2017. a

Khan, S. A., Aschwanden, A., Bjørk, A. A., Wahr, J., Kjeldsen, K. K., and Kjær, K. H.: Greenland ice sheet mass balance: a review, Rep. Prog. Phys., 78, 046801,, 2015. a, b

Khan, S. A., Sasgen, I., Bevis, M., van Dam, T., Bamber, J. L., Wahr, J., Willis, M., Kjær, K. H., Wouters, B., Helm, V., Csatho, B., Fleming, K., Bjørk, A. A., Aschwanden, A., Knudsen, P., and Munneke, P. K.: Geodetic measurements reveal similarities between post-Last Glacial Maximum and present-day mass loss from the Greenland ice sheet, Science Advances, 2, e1600931,, 2016. a, b

Khan, S. A., Bjørk, A. A., Bamber, J. L., Morlighem, M., Bevis, M., Kjær, K. H., Mouginot, J., Løkkegaard, A., Holland, D. M., Aschwanden, A., Zhang, B., Helm, V., Korsgaard, N. J., Colgan, W., Larsen, N. K., Lui, L., Hansen, K., Barletta, V., Dahl-Jensen, T. S., Søndergaard, A. S., Csatho, B. M., Sasgen, I., Box, J., and Schenk, T.: Centennial response of Greenland's three largest outlet glaciers, Nat. Commun., 11, 5718,, 2020. a

King, M. D., Howat, I. M., Jeong, S., Noh, M. J., Wouters, B., Noël, B., and van den Broeke, M. R.: Seasonal to decadal variability in ice discharge from the Greenland Ice Sheet, The Cryosphere, 12, 3813–3825,, 2018. a, b

King, M. D., Howat, I. M., Candela, S. G., Noh, M. J., Jeong, S., Noël, B. P. Y., van den Broeke, M. R., Wouters, B., and Negrete, A.: Dynamic ice loss from the Greenland Ice Sheet driven by sustained glacier retreat, Communications Earth and Environment, 1, 1,, 2020. a, b

Kjeldsen, K. K., Korsgaard, N. J., Bjørk, A. A., Khan, S. A., Box, J. E., Funder, S., Larsen, N. K., Bamber, J. L., Colgan, W., Broeke, M. v. d., Siggaard-Andersen, M.-L., Nuth, C., Schomacker, A., Andresen, C. S., Willerslev, E., and Kjær, K. H.: Spatial and temporal distribution of mass loss from the Greenland Ice Sheet since AD 1900, Nature, 528, 396–400,, 2015. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z

Kjeldsen, K. K., Khan, S. A., Colgan, W. T., MacGregor, J. A., and Fausto, R. S.: Time-Varying Ice Sheet Mask: Implications on Ice-Sheet Mass Balance and Crustal Uplift, J. Geophys. Res.-Earth, 125, e2020JF005775,, 2020. a

Kluyver, T., Ragan-Kelley, B., Pérez, F., Granger, B., Bussonnier, M., Frederic, J., Kelley, K., Hamrick, J., Grout, J., Corlay, S., Ivanov, P., Avila, D., Abdalla, S., and Willing, C.: Jupyter Notebooks – a publishing format for reproducible computational workflows,, 2016.  a

Langen, P. L., Fausto, R. S., Vandecrux, B., Mottram, R. H., and Box, J. E.: Liquid Water Flow and Retention on the Greenland Ice Sheet in the Regional Climate Model HIRHAM5: Local and Large-Scale Impacts, Front. Earth Sci., 4, 110,, 2017. a, b, c, d

Laprise, R.: The Euler Equations of Motion with Hydrostatic Pressure as an Independent Variable, Mon. Weather Rev., 120, 197–207,<0197:teeomw>;2, 1992. a

Lefebre, F.: Modeling of snow and ice melt at ETH Camp (West Greenland): A study of surface albedo, J. Geophys. Res., 108, 4231,, 2003. a

Lenaerts, J. T. M., van den Broeke, M. R., van Angelen, J. H., van Meijgaard, E., and Déry, S. J.: Drifting snow climate of the Greenland ice sheet: a study with a regional climate model, The Cryosphere, 6, 891–899,, 2012. a

Lewis, G., Osterberg, E., Hawley, R., Whitmore, B., Marshall, H. P., and Box, J.: Regional Greenland accumulation variability from Operation IceBridge airborne accumulation radar, The Cryosphere, 11, 773–788,, 2017. a

Lewis, G., Osterberg, E., Hawley, R., Marshall, H. P., Meehan, T., Graeter, K., McCarthy, F., Overly, T., Thundercloud, Z., and Ferris, D.: Recent precipitation decrease across the western Greenland ice sheet percolation zone, The Cryosphere, 13, 2797–2815,, 2019. a

Liestøl, O.: Glacier Dammed Lakes in Norway, vol. 81, Fabritius and Sønners Forlag, Oslo, 1956. a, b

Ligtenberg, S. R. M., Kuipers Munneke, P., Noël, B. P. Y., and van den Broeke, M. R.: Brief communication: Improved simulation of the present-day Greenland firn layer (1960–2016), The Cryosphere, 12, 1643–1649,, 2018. a

Lucas-Picher, P., Wulff-Nielsen, M., Christensen, J. H., Aðalgeirsdóttir, G., Mottram, R., and Simonsen, S. B.: Very high resolution regional climate model simulations over Greenland: Identifying added value, J. Geophys. Res.-Atmos., 117, D02108,, 2012. a

MacGregor, J. A., Fahnestock, M. A., Catania, G. A., Aschwanden, A., Clow, G. D., Colgan, W. T., Gogineni, P. S., Morlighem, M., Nowicki, S. M. J., Paden, J. D., Price, S. F., and Seroussi, H.: A synthesis of the basal thermal state of the Greenland Ice Sheet, J. Geophys. Res.-Earth, 121, 1328–1350,, 2016. a, b

Maier, N., Humphrey, N., Harper, J., and Meierbachtol, T.: Sliding dominates slow-flowing margin regions, Greenland Ice Sheet, Science Advances, 5, eaaw5406,, 2019. a

Maier, N., Gimbert, F., Gillet-Chaulet, F., and Gilbert, A.: Basal traction mainly dictated by hard-bed physics over grounded regions of Greenland, The Cryosphere, 15, 1435–1451,, 2021. a

Mankoff, K.: Greenland Ice Sheet solid ice discharge from 1986 through last month: Gates, GEUS Dataverse [data set],, 2020. a, b

Mankoff, K. and Solgaard, A.: Ice Discharge, GEUS Dataverse [data set],, 2020. a, b, c

Mankoff, K., Fettweis, X., Solgaard, A., Langen, P., Stendel, M., Noël, B., van den Broeke, M. R., Karlsson, N., Box, J. E., and Kjeldsen, K.: Greenland ice sheet mass balance from from 1840 through next week, GEUS Dataverse [data set],, 2021. a, b

Mankoff, K. D.: Total Mass Balance, GitHub [data set], available at:, last access: 4 October 2021. a, b

Mankoff, K. D. and Tulaczyk, S. M.: The past, present, and future viscous heat dissipation available for Greenland subglacial conduit formation, The Cryosphere, 11, 303–317,, 2017. a, b, c

Mankoff, K. D., Noël, B., Fettweis, X., Ahlstrøm, A. P., Colgan, W., Kondo, K., Langley, K., Sugiyama, S., van As, D., and Fausto, R. S.: Greenland liquid water discharge from 1958 through 2019, Earth Syst. Sci. Data, 12, 2811–2841,, 2020a. a

Mankoff, K. D., Solgaard, A., Colgan, W., Ahlstrøm, A. P., Khan, S. A., and Fausto, R. S.: Greenland Ice Sheet solid ice discharge from 1986 through March 2020, Earth Syst. Sci. Data, 12, 1367–1383,, 2020b. a, b, c, d, e, f, g, h, i, j, k, l

Mankoff, K. D., van As, D., Lines, A., Bording, T., Elliott, J., Kraghede, R., Cantalloube, H., Oriot, H., Dubois-Fernandez, P., Ruault du Plessis, O., Vest Christiansen, A., Auken, E., Hansen, K., Colgan, W., and Karlsson, N. B.: Search and recovery of aircraft parts in ice-sheet crevasse fields using airborne and in situ geophysical sensors, J. Glaciol., 66, 496–508,, 2020c. a

Martos, Y. M., Jordan, T. A., Catalán, M., Jordan, T. M., Bamber, J. L., and Vaughan, D. G.: Geothermal Heat Flux Reveals the Iceland Hotspot Track Underneath Greenland, Geophys. Res. Lett., 45, 8214–8222,, 2018. a

McKinney, W.: Data Structures for Statistical Computing in Python, in: Proceedings of the 9th Python in Science Conference, edited by: van der Walt, S. and Millman, J., 28 June–3 July, Austin, Texas, 51–56,, 2010. a

Millan, R., Rignot, E., Mouginot, J., Wood, M., Bjørk, A. A., and Morlighem, M.: Vulnerability of Southeast Greenland glaciers to warm Atlantic Water from Operation IceBridge and Ocean Melting Greenland data, Geophys. Res. Lett., 45, 2688–2696,, 2018. a

Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauché, N., Dowdeswell, J. A., Dorschel, B., Fenty, I., Hogan, K., Howat, I. M., Hubbard, A., Jakobsson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., Mouginot, J., Noël, B. P. Y., Cofaigh, C. Ó., Palmer, S., Rysgaard, S., Seroussi, H., Siegert, M. J., Slabon, P., Straneo, F., van den Broeke, M. R., Weinrebe, W., Wood, M., and Zinglersen, K. B.: BedMachine v3: Complete bed topography and ocean bathymetry mapping of Greenland from multi-beam echo sounding combined with mass conservation, Geophys. Res. Lett., 44, 11051–11061,, 2017. a

Morlighem, M., Williams, C., Rignot, E., An, L., Arndt, J. E., Bamber, J., Catania, G., Chauché, N., Dowdeswell, J. A., Dorschel, B., Fenty, I., Hogan, K., Howat, I., Hubbard, A., Jakobsson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., Mouginot, J., Noël, B., O'Cofaigh, C., Palmer, S. J., Rysgaard, S., Seroussi, H., Siegert, M. J., Slabon, P., Straneo, F., van den Broeke, M. R., Weinrebe, W., Wood, M., and Zinglersen, K.: IceBridge BedMachine Greenland, Version 4, all subsets used, NASA National Snow and Ice Data Center Distributed Active Archive Center [data set], Boulder, Colorado, USA,, 2021. a, b

Mouginot, J. and Rignot, E.: Glacier catchments/basins for the Greenland Ice Sheet, Dryad [data set],, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p

Mouginot, J., Rignot, E., Bjørk, A. A., van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B., and Wood, M.: Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018, P. Natl. Acad. Sci. USA, 116, 9239–9244,, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p

Neteler, M., Bowman, M. H., Landa, M., and Metz, M.: GRASS GIS: a multi-purpose Open Source GIS, Environ. Modell. Softw., 31, 124–130,, 2012. a

Nilsson, J., Vallelonga, P., Simonsen, S. B., Sørensen, L. S., Forsberg, R., Dahl-Jensen, D., Hirabayashi, M., Goto-Azuma, K., Hvidberg, C. S., Kjær, H. A., and Satow, K.: Greenland 2012 melt event effects on CryoSat-2 radar altimetry, Geophys. Res. Lett., 42, 3919–3926, 2015. a

Noël, B., van de Berg, W. J., Machguth, H., Lhermitte, S., Howat, I., Fettweis, X., and van den Broeke, M. R.: A daily, 1 km resolution data set of downscaled Greenland ice sheet surface mass balance (1958–2015), The Cryosphere, 10, 2361–2377,, 2016. a

Noël, B., van de Berg, W. J., van Wessem, J. M., van Meijgaard, E., van As, D., Lenaerts, J. T. M., Lhermitte, S., Kuipers Munneke, P., Smeets, C. J. P. P., van Ulft, L. H., van de Wal, R. S. W., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 1: Greenland (1958–2016), The Cryosphere, 12, 811–831,, 2018. a

Noël, B., van de Berg, W. J., Lhermitte, S., and van den Broeke, M. R.: Rapid ablation zone expansion amplifies north Greenland mass loss, Science Advances, 5, eaaw0123,, 2019. a, b, c

Oliphant, T. E.: A guide to NumPy, vol. 1, Trelgol Publishing USA, United States, 2006. a

Pfeffer, W. T., Meier, M. F., and Illangasekare, T. H.: Retention of Greenland runoff by refreezing: Implications for projected future sea level change, J. Geophys. Res., 96, 22117,, 1991. a

Rastner, P., Bolch, T., Mölg, N., Machguth, H., Le Bris, R., and Paul, F.: The first complete inventory of the local glaciers and ice caps on Greenland, The Cryosphere, 6, 1483–1495,, 2012. a

Rezvanbehbahani, S., Stearns, L. A., van der Veen, C. J., Oswald, G. K. A., and Greve, R.: Constraining the geothermal heat flux in Greenland at regions of radar-detected basal water, J. Glaciol., 65, 1023–1034,, 2019. a

Rignot, E., Mouginot, J., Scheuchl, B., van den Broeke, M., van Wessem, M. J., and Morlighem, M.: Four decades of Antarctic Ice Sheet mass balance from 1979–2017, P. Natl. Acad. Sci. USA, 116, 1095–1103,, 2019. a

Rignot, E. J., Box, J. E., Burgess, E. W., and Hanna, E.: Mass balance of the Greenland ice sheet from 1958 to 2007, Geophys. Res. Lett., 35, L20502,, 2008. a, b

Rocklin, M.: Dask: Parallel Computation with Blocked algorithms and Task Scheduling, in: Proceedings of the 14th Python in Science Conference, edited by: Huff, K. and Bergstra, J., 6–12 July, Austin, TX, 130–136, 2015. a

Roeckner, E., Bäuml, G., Bonaventura, L., Brokopf, R., Esch, M., Giorgetta, M., Hagemann, S., Kirchner, I., Kornblueh, L., Manzini, E., Rhodin, A., Schlese, U., Schulzweida, U., and Tompkins, A.: The atmospheric general circulation model ECHAM 5. Part I: Model description, Report No. 349, Max-Planck-Institut für Meteorologie, 2003.  a

Rogozhina, I., Hagedoorn, J. M., Martinec, Z., Fleming, K., Soucek, O., Greve, R., and Thomas, M.: Effects of uncertainties in the geothermal heat flux distribution on the Greenland Ice Sheet: An assessment of existing heat flow models, J. Geophys. Res.-Earth, 117, F02025,, 2012. a

Ryser, C., Lüthi, M. P., Andrews, L. C., Hoffman, M. J., Catania, G. A., Hawley, R. L., Neumann, T. A., and Kristensen, S. S.: Sustained high basal motion of the Greenland ice sheet revealed by borehole deformation, J. Glaciol., 60, 647–660,, 2014. a

Sasgen, I., Van Den Broeke, M., Bamber, J. L., Rignot, E., Sørensen, L. S., Wouters, B., Martinec, Z., Velicogna, I., and Simonsen, S. B.: Timing and origin of recent regional ice-mass loss in Greenland, Earth Planet. Sc. Lett., 333–334, 293–303,, 2012. a

Sasgen, I., Wouters, B., Gardner, A. S., King, M. D., Tedesco, M., Landerer, F. W., Dahle, C., Save, H., and Fettweis, X.: Return to rapid ice loss in Greenland and record loss in 2019 detected by the GRACE-FO satellites, Communications Earth and Environment, 1, 293–303,, 2020. a

Schulte, E., Davison, D., Dye, T., and Dominik, C.: A multi-language computing environment for literate programming and reproducible research, J. Stat. Softw., 46, 1–24, 2012. a

Schulzweida, U.: CDO User Guide, Zenodo [code],, 2019. a

Shapiro, N. M. and Ritzwoller, M. H.: Inferring Surface Heat Flux Distributions Guided by a Global Seismic Model: Particular Application to Antarctica, Earth Planet. Sc. Lett., 223, 213–224,, 2004. a

Shreve, R. L.: Movement of water in glaciers, J. Glaciol., 11, 205–214, 1972. a

Simmons, A. J. and Burridge, D. M.: An Energy and Angular-Momentum Conserving Vertical Finite-Difference Scheme and Hybrid Vertical Coordinates, Mon. Weather Rev., 109, 758–766,<0758:aeaamc>;2, 1981. a

Simonsen, S. B., Barletta, V. R., Colgan, W., and Sørensen, L. S.: Greenland Ice Sheet mass balance (1992-2020) from calibrated radar altimetry, Geophys. Res. Lett., 48, e2020GL091216,, 2021a. a, b, c, d, e, f

Simonsen, S. B., Barletta, V. R., Colgan, W., and Sørensen, L. S.: Greenland Ice Sheet mass balance (1992–2020) from calibrated radar altimetry, DTU [data set],, 2021b. a

Smith, B., Fricker, H. A., Gardner, A. S., Medley, B., Nilsson, J., Paolo, F. S., Holschuh, N., Adusumilli, S., Brunt, K., Csatho, B., Harbeck, K., Markus, T., Neumann, T., Siegfried, M. R., and Zwally, J. H.: Pervasive ice sheet mass loss reflects competing ocean and atmosphere processes, Science, 368, eaaz5845,, 2020. a

Solgaard, A. and Kusk, A.: Greenland Ice Velocity from Sentinel-1 Edition 2, GEUS Dataverse [data set],, 2021. a, b

Solgaard, A., Kusk, A., Merryman Boncori, J. P., Dall, J., Mankoff, K. D., Ahlstrøm, A. P., Andersen, S. B., Citterio, M., Karlsson, N. B., Kjeldsen, K. K., Korsgaard, N. J., Larsen, S. H., and Fausto, R. S.: Greenland ice velocity maps from the PROMICE project, Earth Syst. Sci. Data, 13, 3491–3512,, 2021. a, b, c, d

Stallman, R. M.: EMACS the Extensible, Customizable Self-Documenting Display Editor, in: Proceedings of the ACM SIGPLAN SIGOA Symposium on Text Manipulation, Association for Computing Machinery, New York, NY, USA, 147–156,, 1981. a

Sutterley, T. C., Velicogna, I., Rignot, E. J., Mouginot, J., Flament, T., van den Broeke, M. R., van Wessem, J. M., and Reijmer, C. H.: Mass loss of the Amundsen Sea Embayment of West Antarctica from four independent techniques, Geophys. Res. Lett., 41, 8421–8428,, 2014. a

Sørensen, L. S., Simonsen, S. B., Nielsen, K., Lucas-Picher, P., Spada, G., Adalgeirsdottir, G., Forsberg, R., and Hvidberg, C. S.: Mass balance of the Greenland ice sheet (2003–2008) from ICESat data – the impact of interpolation, sampling and firn density, The Cryosphere, 5, 173–186,, 2011. a

Tange, O.: GNU Parallel – The Command-Line Power Tool, ;login: The USENIX Magazine, 36, 42–47, available at: (last access: 24 October 2021), 2011. a

Tedesco, M. and Fettweis, X.: Unprecedented atmospheric conditions (1948–2019) drive the 2019 exceptional melting season over the Greenland ice sheet, The Cryosphere, 14, 1209–1223,, 2020. a

The IMBIE Team: Mass balance of the Greenland Ice Sheet from 1992 to 2018, Nature, 579, 233–239,, 2019. a, b, c, d, e

Uppala, S. M., Kållberg, P. W., Simmons, A. J., Andrae, U., Da Costa Bechtold, V., Fiorino, M., Gibson, J., Haseler, J., Hernandez, A., Kelly, G. A., Li, X., Onogi, K., Saarinen, S., Sokka, N., Allan, R. P., Anderson, E., Arpe, K., Balmaseda, M. A., Beljaars, A. C. M., Van De Berg, L., Bidlot, J., Bormann, N., Caires, S., Chevallier, F., Dethof, A., Dragosavac, M., Fisher, M., Fuentes, M., Hagemann, S., Hólm, E., Hoskins, B. J., Isaksen, L., Janssen, P. A. E. M., Jenne, R., Mcnally, A. P., Mahfouf, J.-F., Morcrette, J.-J., Rayner, N. A., Saunders, R. W., Simon, P., Sterl, A., Trenbreth, K. E., Untch, A., Vasiljevic, D., Viterbo, P., and Woollen, J.: The ERA-40 re-analysis, Q. J. Roy. Meteor. Soc., 131, 2961–3012,, 2005. a

van Angelen, J. H., Lenaerts, J. T. M., Lhermitte, S., Fettweis, X., Kuipers Munneke, P., van den Broeke, M. R., van Meijgaard, E., and Smeets, C. J. P. P.: Sensitivity of Greenland Ice Sheet surface mass balance to surface albedo parameterization: a study with a regional climate model, The Cryosphere, 6, 1175–1186,, 2012. a

van de Berg, W. J. and Medley, B.: Brief Communication: Upper-air relaxation in RACMO2 significantly improves modelled interannual surface mass balance variability in Antarctica, The Cryosphere, 10, 459–463,, 2016.  a

Van Rossum, G. and Drake Jr., F. L.: Python reference manual, Centrum voor Wiskunde en Informatica Amsterdam, 1995. a

Velicogna, I., Mohajerani, Y., A, G., Landerer, F., Mouginot, J., Noel, B., Rignot, E., Sutterley, T., Broeke, M., Wessem, M., Rignot, E., Sutterley, T., van den Broeke, M., van Wessem, M., and Wieseand, D.: Continuity of Ice Sheet Mass Loss in Greenland and Antarctica From the GRACE and GRACE Follow-On Missions, Geophys. Res. Lett., 47, e2020GL087291,, 2020. a

Vinther, B. M., Andersen, K. K., Jones, P. D., Briffa, K. R., and Cappelen, J.: Extending Greenland temperature records into the late eighteenth century, J. Geophys. Res., 111, D11105,, 2006. a

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791,, 2012. a

Wang, Z. Q. and Randriamampianina, R.: The Impact of Assimilating Satellite Radiance Observations in the Copernicus European Regional Reanalysis (CERRA), Remote Sens.-Basel, 13, 426,, 2021. a

Yang, X., Palmason, B., Sattler, K., Thorsteinsson, S., Amstrup, B., Dahlbom, M., Hansen Sass, B., Pagh Nielsen, K., and Petersen, G. N.: IGB, the Upgrade to the Joint Operational HARMONIE by DMI and IMO in 2018, ALADIN-HIRLAM Newsletter, 93–96, 2018. a

Zemp, M., Huss, M., Eckert, N., Thibert, E., Paul, F., Nussbaumer, S. U., and Gärtner-Roer, I.: Brief communication: Ad hoc estimation of glacier contributions to sea-level rise from the latest glaciological observations, The Cryosphere, 14, 1043–1050,, 2020. a

Zender, C. S.: Analysis of self-describing gridded geoscience data with netCDF Operators (NCO), Environ. Modell. Softw., 23, 1338–1342,, 2008. a

Zuo, Z. and Oerlemans, J.: Modelling albedo and specific balance of the Greenland ice sheet: calculations for the Søndre Strømfjord transect, J. Glaciol., 42, 305–317,, 1996. a

Zwally, H. J. and Giovinetto, M. B.: Overview and Assessment of Antarctic Ice-Sheet Mass Balance Estimates: 1992–2009, Surv. Geophys., 32, 1–26,, 2011. a

Zwally, H. J., Giovinetto, M. B., Beckley, M. A., and Saba, J. L.: Antarctic and Greenland Drainage Systems, available at: (last access: 27 April 2020), GSFC Cryospheric Sciences Laboratory, 2012. a, b, c, d, e, f, g, h, i, j, k

Short summary
We estimate the daily mass balance and its components (surface, marine, and basal mass balance) for the Greenland ice sheet. Our time series begins in 1840 and has annual resolution through 1985 and then daily from 1986 through next week. Results are operational (updated daily) and provided for the entire ice sheet or by commonly used regions or sectors. This is the first input–output mass balance estimate to include the basal mass balance.