ms.org

Table of contents

• About this document
  • Workflow
• Introduction
  • Terminology
• Product description
• Data sources
  • Surface mass balance
  • Discharge
  • Basal mass balance
  • Other
  • Products used for validation
• Methods
  • Surface mass balance
  • Projected discharge
  • Basal mass balance
  • Reconstructed adjustment
  • Domains, boundaries, and regions of interest
• Product evaluation and assessment
  • Mouginot (2019)
  • Colgan (2019)
  • GMB
  • VC
  • IMBIE
• Uncertainty
  • Discharge
  • ROIs
  • Accumulating uncertainties
  • Peripheral ice masses
• Results
• Data availability
• Conclusions
• Appendix
  • RCM differences
  • Mouginot 2019 by region
  • Colgan 2019 by sector
  • Reconstructed runoff
  • RCM coverage
  • Software
  • CRediT
• Misc journal sections
• References

About this document

This document is an Emacs Org Mode plain-text file with code and text embedded. If you are viewing:

• A DOC, Google Doc, or PDF file, then it was generated by exporting from Org. Not all of the Org parts (code, results, comments, etc.) were exported. The Org source file is available upon request, and may be embedded in the PDF. Most non-Apple PDF viewers provide easy access to embedded or attached files.
• A webpage somewhere, then this is a subset of the code and text that the website render has decided to display to you through the browser. You can choose to view the raw source and/or download it and view it locally on your computer.
• A file with a org extension in something other than Emacs, then you are seeing the canonical version and the full source, but without any syntax highlighting, document structure, or the ability to execute the code blocks.
• An Org file within Emacs, then this is the canonical version. You should be able to fully interact and reproduce the contents of this document, although it may require 3rd-party applications (Python, etc.) a similar Emacs configuration, and the data files. This is available upon request.

Workflow

To recreate this work

• Open this file in Emacs Org Mode.
• check that you have the necessary software dependencies installed. See section: Code.
• Download and set up the necessary data files as per the Data section
• Tangle the embedded code blocks.
  • Execute C-c C-v C-t to run the (org-babel-tangle) function
• Run make
  • This should probably be run in an external terminal because it takes hours to days…
• Update Babel result blocks throughout the document by
  • Cleaning all result blocks with C-u C-c C-v k or (org-babel-remove-result-one-or-many t), then
  • Executing all blocks (without :eval no) using C-c C-v C-b or (org-babel-execute-buffer)

Introduction

\introduction

Over the past several decades, mass loss from the Greenland ice sheet has increased citep:khan_2015,the-imbie-team_2019. Different processes dominate the regional mass loss of the ice sheet, and their relative contribution has fluctuated in time citep:mouginot_2019_data. 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. \ref{fig:overview}).

There are three common methods for estimating mass balance – changes in gravity citep:barletta_2013,groh_2019,the-imbie-team_2019,velicogna_2020, changes in volume citep:simonsen_2021,sorensen_2011,zwally_2011,sasgen_2012,smith_2020, and the input-output (IO) method citep:colgan_2019,mouginot_2019,rignot_2019,king_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 ~100 km 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 ~1 km 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 days year-round after the launch of the Sentinel missions citep:solgaard_2021. The primary issue with the IO method is unknown ice thickness in some locations (e.g. citet:mankoff_2020_solid). 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 citep:karlsson_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 the Monitoring of the Greenland ice sheet (PROMICE) Greenland ice sheet mass balance data set based on the IO method, updating the previous product from citet:colgan_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 citep:zwally_2012_data or region citep:mouginot_2019_data. 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 days with ~12 day 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.

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 citet:kjeldsen_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 citet:kjeldsen_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. \ref{fig:overview}).
• Sector refers to one of the citet:zwally_2012_data sectors (Fig. \ref{fig:overview}), expanded here to cover the RCM ice domains which exist slightly outside these sectors in some locations.
• Region refers to the citet:mouginot_2019_data regions (Fig. \ref{fig:overview}), 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 (BMB_GF), frictional heating from ice velocity (BMB_friction), and viscous heat dissipation (BMB_VHD).
• MB is the total mass balance including the BMB term (Eq. \ref{eq:MB}).
• MB* is the mass balance not including the BMB term (Eq. \ref{eq:MBstar}).
• 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.

Product description

The output of this work is two NetCDF files and two CSV files containing 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 citet:zwally_2012_data sectors and one for citet:mouginot_2019_data 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 contains a copy of the ice-sheet-summed data, one daily and one annual.

An example of the output is shown in Fig. \ref{fig:mb_ts}a, which shows mass balance for the entire Greenland ice sheet in addition to SMB and D at annual resolution. Figure \ref{fig:mb_ts}b 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. \ref{fig:overview}).

Data sources

This section introduces data products that exist prior to and are external to this work (Table \ref{tab:data}). 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.

The inputs to this work are the recent SMB fields from the three RCMs, the recent discharge from citet:mankoff_2020_solid (data: citet:mankoff_D_dataverse), and the recent basal mass balance fields, of which BMB_GF and BMB_friction are direct outputs from citet:karlsson_2021 (data: citet:karlsson_2021_data), but the BMB_VHD calculations are redone here (see Methods Sect. \ref{sec:methods:VHD}) using the MAR runoff field. The reconstructed data (pre-1986) are surface mass balance and discharge from citet:kjeldsen_2015 (data: citet:kjeldsen_2015_data), but adjusted here using the overlapping period (see Methods Sect. \ref{sec:methods:reconstructed_adjustment}), and runoff from citet:kjeldsen_2015 (data: citet:kjeldsen_2015_data) as a proxy and scaled for BMB_VHD (see Methods Sect. \ref{sec:methods:VHD}).

Product	Period	Reference	Data/notes
Reconstructed SMB	1840 through 1985	citet:kjeldsen_2015	citet:kjeldsen_2015_data
Reconstructed D	1840 through 1985	citet:kjeldsen_2015	citet:kjeldsen_2015_data
HIRHAM/HARMONIE SMB	1986 through yesterday	citet:langen_2017
MAR SMB	1986 through next week	citet:fettweis_2020
RACMO SMB	1986 through last month	citet:noel_2019
D	1986 through last month	citet:mankoff_2020_solid	citet:mankoff_D_dataverse
BMBGF; BMBfriction	1840 through next week	citet:karlsson_2021	citet:karlsson_2021_data
BMBVHD	1840 through 1985	citet:kjeldsen_2015	citet:kjeldsen_2015_data reconstructed runoff
BMBVHD	1986 through next week	citet:fettweis_2020	MAR runoff

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).

HIRHAM/HARMONIE

The HIRHAM/HARMONIE product from the Danmarks Meteorologiske Institut (Danish Meteorological Institute; DMI) is based on an offline subsurface firn/SMB model citep:langen_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 2017-08-31 after which HARMONIE is used.

The HIRHAM5 regional climate model citep:christensen_2007 combines the dynamical core of the HIRLAM7 numerical weather forecasting model citep:eerola_2006 with physics schemes from the ECHAM5 general circulation model citep:roeckner_2003. In the Greenland setup employed here citep:lucas-picher_2012, it has a horizontal resolution of 0.05 ° x 0.05 ° on a rotated pole grid (corresponding to 5.5 km resolution), and 31 atmospheric levels. It is forced at 6 hour intervals on the lateral boundaries with horizontal wind vectors, temperature, and specific humidity from the ERA-Interim reanalysis citep:dee_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 citep:langen_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 citep:zuo_1996,lefebre_2003. The offline subsurface model is run on the HIRHAM5 5.5 km grid.

HARMONIE citep:bengtsson_2017 is a nonhydrostatic model in terrain-following sigma coordinates based on the fully compressible Euler equations citep:simmons_1981,laprise_1992. 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 citep:wang_2021,yang_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 citet:citterio_2013 ice mask interpolated to the 5.5 km grid.

MAR

The MAR RCM has been developed by the University of Liège (Belgium) with a focus on the polar regions citep:fettweis_2020. The MAR atmosphere module citep:gallee_1994 is fully coupled with the soil-ice-snow energy balance vegetation model SISVAT citep:gallee_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 hour intervals by ERA5 reanalysis and runs at 20 km resolution. The snowpack was initialised in 1950 from a former MARv3.11-based simulation. Its snow model is based on a former version of the CROCUS snow model citep:vionnet_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 citep:howat_2014.

We use MAR version 3.12. With respect to version 3.9, intensively validated over Greenland citep:fettweis_2020 or the 20 km-based MARv3.10 setup used in citet:tedesco_2020, MARv3.12 now uses the common polar stereographic projection EPSG 3413. With respect to MARv3.11 fully described in citet:amory_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 citet:fettweis_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 BMB_VHD estimate (Section \ref{sec:methods:VHD}).

As the recent SMB decrease (successfully evaluated with GRACE-based estimates in citet:fettweis_2020) has been fully driven by the increase in runoff citep:sasgen_2020, we assume the same degree of accuracy between SMB simulated by MAR (evaluated with the PROMICE SMB database citep:fettweis_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 h 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 days. At the end of each day, ERA5 is used to update the GFS-forced MAR time series until about 5 days 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 BMB_VHD) fields.

RACMO

The 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 citet:noel_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 citep:ettema_2010. The snow layers were initialized in September 1957 using temperature and density from a previous run with the offline IMAU Firn Densification Model citep:ligtenberg_2018. To simulate drifting snow transport and sublimation, citet:lenaerts_2012 implemented a drifting snow scheme. Snow albedo depends on snow grain size, cloud optical thickness, solar zenith angle, and impurity content citep:van-angelen_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-day 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) citep:howat_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 citep:uppala_2005, ERA-Interim between 1979 and 1989 citep:dee_2011, and ERA-5 between 1990 and 2020 citep:hersbach_2020. The relaxation zone is 24 grid cells (~130 km) wide to ensure a smooth transition to the domain interior. This run has active upper-atmosphere relaxation citep:berg_2016. 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 citep:noel_2019, employing corrections for biases in elevation and bare ice albedo using a MODIS albedo product at 1 km resolution citep:noel_2016.

Reconstruction

The citet:kjeldsen_2015 173-year (1840 through 2012) mass balance reconstruction is based on the citet:box_2013_II 171-year (1840 through 2010) statistical reconstruction. citet:kjeldsen_2015 add a more sophisticated meltwater retention scheme citep:pfeffer_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 citet:box_2013_II 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 citep:cappelen_2011,cappelen_2001,cappelen_2006,vinther_2006 or firn/ice cores citep:box_2013_I and b) spatially continuous outputs from the regional climate model RACMO version 2.1 citep:ettema_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 citet:fettweis_2017. citet:box_2013_I describe the methods in more detail.

The reconstructed surface mass balance is adjusted as described in the Methods Sect. \ref{sec:methods:reconstructed_adjustment} (Fig. \ref{fig:reconstructed_adj}).

Discharge

The recent discharge data are from citet:mankoff_2020_solid (data: citet:mankoff_D_dataverse). This product covers all fast-flowing (> 100 m yr^-1) marine-terminating glaciers. The discharge in citet:mankoff_2020_solid is computed at flux gates ~5 km upstream from glacier termini citep:mankoff_D_gates, 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 ~12 days with a ~30 day lag (citet:solgaard_2021; data: citet:solgaard_2021_data).

Some changes have been implemented since the last publication describing the discharge product (i.e., citet:mankoff_2020_solid). These are minor and include updating the citet:khan_2016 (data: citet:khan_2016_data) 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 (doi:10.22008/promice/data/sentinel1icevelocity/greenlandicesheet/v1.0.0) to Edition 2 (citet:solgaard_2021,solgaard_2021_data), and updating from BedMachine v3 (supplemented in the SE with citet:millan_2018) to use only BedMachine v4 citep:NSIDC_BedMachine_GL.

The reconstructed discharge data citep:kjeldsen_2015 are estimated via a linear fit between unsmoothed annual discharge spanning 2000 to 2012 citep:enderlin_2014_improved and runoff data from citet:kjeldsen_2015 using a 6-year trailing average. The method for scaling discharge from runoff was introduced by citet:rignot_2008_mass, who scaled the SMB anomaly with discharge. Sensitivity analyses conducted by citet:box_2013_III 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., citet:boers_2021), the discharge data used by citet:rignot_2008_mass and citet:box_2013_III also include years 1958 and 1964 that lie near the regression line (see citet:box_2013_III Fig. 4, and the related Sect. 4, Physical basis). Further, while 2000 through 2005 cover a changing period in Greenlandic discharge citep:mankoff_2020_solid,king_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 citep:andresen_2012,bjork_2012,khan_2015,khan_2020.

The reconstructed discharge is adjusted as described in the Methods Sect. \ref{sec:methods:reconstructed_adjustment}.

Basal mass balance

The BMB citep:karlsson_2021 comes from mass lost at the bed from BMB_GF, BMB_friction from the basal shear velocity, and BMB_VHD from surface runoff routed to the bed (i.e. the volume of the subglacial conduits formed from surface runoff; citet:mankoff_2017_VHD).

These fields (data: citet:karlsson_2021_data) are provided as steady-state annual estimates. We use the BMB_GF and BMB_friction products and apply 1/365th to each day, each year. Because BMB_VHD is proportional to runoff, an annual estimate is not appropriate for this work with daily resolution. We therefore re-calculate the BMB_VHD-induced basal melt as described in the Methods Sect. \ref{sec:methods:VHD}.

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 BMB_GF often with diverging results (see e.g., citet:rogozhina_2012,rezvanbehbahani_2019). Lacking substantial validation that favors one BMB_GF map over the others, citet:karlsson_2021 instead use the average of three widely used BMB_GF estimates: citet:foxmaule_2009,shapiro_2004, and citet:martos_2018. The BMB_GF melt rate is calculated as

\begin{equation} \dot{b}_m = E_GF \, ρ_i^-1 \, L^-1, \end{equation}

where \(E_GF\) 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; citet:cuffey_2010). BMB_GF melting is only calculated where the bed is not frozen. We use the citet:macgregor_2016 estimate of temperate bed extent and scale Eq. \ref{eq:GF} by 0, 0.5, or 1 where the bed is frozen (~25 % of the ice sheet area), uncertain (~33 %), or thawed (~42 %), respectively.

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., citet:ryser_2014_sustained,maier_2019) and it is unclear how representative those measurements are at ice sheet scales. citet:karlsson_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 citep:gillet-chaulet_2012,maier_2021. The model is tuned to match a multi-decadal surface velocity map citep:joughin_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. \ref{eq:GF}:

\begin{equation} \dot{b}_m = E_f \, ρ_i^-1 \, L^-1, \end{equation}

where \(E_f\) is energy due to friction. We also apply the 0, 0.5, and 1 scale as used for the BMB_GF term citep:macgregor_2016 in order to mask out areas that are likely frozen.

Other

ROI regions come from citet:mouginot_2019_data and ROI sectors come from citet:zwally_2012_data.

Products used for validation

We validate This Study against five other data products (see Table \ref{tab_summary} and Sect. \ref{sec:validation}). These products are the most recent IO product citep:mouginot_2019, the previous PROMICE mass balance product (citet:colgan_2019; data: citet:colgan_2019_data), the two mostly independent methods of estimating ice sheet mass change, GMB (citet:barletta_2013; data: citet:barletta_2013_data) and VC (citet:simonsen_2021; data: citet:simonsen_2021_data), and the IMBIE2 data citep:the-imbie-team_2019. In addition to this, we evaluate the reconstructed citet:kjeldsen_2015 (data: citet:kjeldsen_2015_data) and This Study data during the overlapping period 1986 through 2012.

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

\begin{equation} MB = SMB - D - BMB. \end{equation}

Products that do not include the BMB term (i.e., citet:mouginot_2019,colgan_2019, and citet:kjeldsen_2015) have total mass balance defined as

\begin{equation} MB^* = SMB - D, \end{equation}

and when comparing This Study to those products, we compare like terms, never comparing our MB to a different product MB^*, except in Fig. \ref{fig:mb_cumsum}, where all products are shown together.

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

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. \ref{sec:methods:reconstructed_adjustment}), 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 \ref{apdx:RCM_diff} 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.

Projected discharge

We project the discharge from the last observed point from citet:mankoff_2020_solid (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 citep:king_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.

Basal mass balance

Because citet:karlsson_2021 provide a steady-state annual-average estimate of the BMB fields, we divide the BMB_GF and BMB_friction fields by 365 to estimate daily average. This is a reasonable treatment of the BMB_GF field, which does not have an annual cycle. The BMB_friction 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 citep:king_2018, and citet:karlsson_2021 found that basal melt rates are 5 % higher during the summer. Thus, the intra-annual changes are less than the uncertainty. The BMB_VHD field varies significantly throughout the year, because it is proportional to surface runoff. We therefore generate our own BMB_VHD for This Study.

To estimate recent BMB_VHD we use daily MAR runoff (see citet:mankoff_2020_liquid) and BedMachine v4 citep:morlighem_2017,NSIDC_BedMachine_GL to derive subglacial routing pathways, similarly to citet:mankoff_2017_VHD. 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

\begin{equation} h = z_b + k \frac{ρ_i}{ρ_w} (z_s - z_b), \end{equation}

with \(z_b\) 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 \(z_s\) the ice surface. Eq. eq:head comes from citet: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 citep:liestol_1956,mankoff_2017_VHD. We assume all energy, \(E_\mathrm{VHD}\) (in Joules), is used to melt ice with

\begin{equation} b_m = E_\mathrm{VHD} \, ρ_i^-1 \, L^-1. \end{equation}

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 BMB_GF and BMB_friction as steady state as described at the start of this section. For BMB_VHD we use the fact that VHD comes from runoff by definition, and from this, reconstructed BMB_VHD 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 % citep:liestol_1956. To estimate the scale factor we use the 1986 through 2012 overlap between the citet:kjeldsen_2015 runoff and This Study recent BMB_VHD from MAR runoff described above. The correlation between the two has an r^2 value of 0.75, a slope of 0.03, and an intercept of -3 Gt yr^-1 (Appendix \ref{apdx:reconstructed_runoff}). From this, we scale the citet:kjeldsen_2015 reconstructed runoff by 3 % (from the 0.03 slope, unrelated to the theoretical 1000 m drop described earlier) to estimate reconstructed BMB_VHD.

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. \ref{fig:reconstructed_adj} SMB), and This Study vs reconstructed D has a slope of 1.1 and an intercept of -17 Gt yr^-1 (Fig. \ref{fig:reconstructed_adj} 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. \ref{fig:reconstructed_adj} 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. \ref{fig:reconstructed_adj}). We then derive the BMB_VHD term for reconstructed basal mass balance (Sect. \ref{sec:methods:VHD} and Appendix \ref{apdx:reconstructed_runoff}), bring in the other BMB terms (Sect. \ref{sec:methods:VHD}), and use Eq. \ref{eq:MB} 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 citep:fettweis_2020, an offset resulting from differences in ice masks citep:kjeldsen_2015, the inclusion of peripheral glaciers citep:kjeldsen_2015, other accumulation rate inaccuracies citep:lewis_2017,lewis_2019, and other unknowns.

Domains, boundaries, and regions of interest

Few of the ice masks used here are spatially aligned. The citet:zwally_2012_data sectors and the citet:mouginot_2019_data regions are often smaller than the RCM ice domains. For example, the RACMO ice domain is 1,718,959 km², of which 1,696,419 km² (99 %) are covered by the citet:mouginot_2019_data regions and 22,540 km² (1 %) are not, or 1,678,864 km² (98 %) are covered by citet:zwally_2012_data and 40,095 km² (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 citet:mouginot_2019_data has a low of 76 Gt (68 % high) and a high of 429 Gt (2 % high). RACMO cropped to citet:zwally_2012_data 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. \ref{fig:overview}) of the nearest non-null cell, and then clipping to the RCM ice domain. This is done for each ROI and RCM. Appendix \ref{apdx:RCM_coverage} provides a graphical display of the HIRHAM RCM domain, the citet:mouginot_2019_data domain, and our expanded citet:mouginot_2019_data domain.

BMB_VHD 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.

Product evaluation and assessment

We compare to six related data sets (see Table \ref{tab_summary} and Sect. \ref{sec:validation_products}): the most similar and recent IO product citep:mouginot_2019, the previous PROMICE assessment citep:colgan_2019, the two mostly independent methods (GMB citep:barletta_2013 and VC citep:simonsen_2021), IMBIE2 citep:the-imbie-team_2019, and the unadjusted reconstructed/recent overlap citep:kjeldsen_2015.

Our initial comparison (Fig. \ref{fig:mb_cumsum}) 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, citet:mouginot_2019) until 7 d from now (now defined as src_bash{date –iso} 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 citet:mouginot_2019 and citet:colgan_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 \ref{tab_summary}. All have different masks. Bias [Gt yr^-1] is defined as \(\frac{1}{n} ∑_i=1^n(x_i - y_i)\). Root mean square error (RMSE) (Gt yr^-1) is defined as \(\sqrt{\frac{1}{n} ∑_i=1^n(x_i - y_i)^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.

Other product	r^2	Bias	RMSE	Fig.	Overlap	Notes
citet:mouginot_2019	0.94	12	37	\ref{fig:M2019}	1986 – 2018	No basal mass balance
citet:colgan_2019	0.87	-32	59	\ref{fig:C2019}	1995 – 2015	No basal mass balance
GMB	0.86	32	63	\ref{fig:GMB_VC_IMBIE}	2002 – 2020	Includes peripheral masses
VC	0.62	-11	86	\ref{fig:GMB_VC_IMBIE}	1992 – 2019	Multi-year smooth
IMBIE2	0.89	-7	44	\ref{fig:GMB_VC_IMBIE}	1992 – 2018	No BMB when using IO; BMB when using GMB or VC
citet:kjeldsen_2015	0.80	6	60	\ref{fig:reconstructed_adj}	1986 – 2012	No basal mass balance; includes peripheral masses

Mouginot (2019)

The citet:mouginot_2019 product spans the 1972 through 2018 period. We only use 1986 and onward because This Study has annual resolution prior to 1986 and citet:mouginot_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 (r² call_tso(row=”mouginot”, col=”r”) {{{results(0.94)}}}, bias call_tso(row=”mouginot”, col=”bias”) {{{results(12)}}}, RMSE call_tso(row=”mouginot”, col=”RMSE”) {{{results(37)}}} , slope 1.1). The minor SMB differences are likely due to mask differences or our use of a three-RCM average SMB estimate.

citet:mouginot_2019 discharge and our D from citet:mankoff_2020_solid 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. \ref{fig:M2019} shifting values to the right), it adds ~25 Gt yr^-1 to This Study.

Because MB is a linear combination of SMB and D terms (Eq \ref{eq:MBstar}), the MB differences between this product and citet:mouginot_2019 are dominated by the D term, although it is not apparent because interannual variability is dominated by SMB.

Colgan (2019)

The citet:colgan_2019 product spans 1995 through 2015. The SMB term is broadly similar to the RCM-averaged SMB term in This Study, although citet:colgan_2019 use only an older version of MAR (Fig. \ref{fig:C2019} top panel). The citet:colgan_2019 SMB is spatially interpolated over the PROMICE ice sheet ice mask citep:citterio_2013, 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.

The largest difference between This Study and citet:colgan_2019 is that the latter estimate grounding line ice discharge based on corrections to ice volume flow rate measured across the ~1700 m elevation contour. This is far inland relative to the grounding line flux gates used in This Study (from citet:mankoff_D_gates). This introduces uncertainty into the citet:colgan_2019 D term from SMB corrections between the 1700 m elevation contour and the terminus (see the large disagreement in Fig. \ref{fig:C2019} 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 citet:colgan_2019. While citet:colgan_2019 did not assess BMB, the majority of this discrepancy likely results from citet:colgan_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, citet:colgan_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 citet:colgan_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 citet:colgan_2019 particularly overestimated the response in ice discharge to 1990s climate variability.

Similarly to the comparison with citet:mouginot_2019, the disagreement between This Study and citet:colgan_2019 is dominated by D disagreement, although it is again not apparent because interannual variability is dominated by SMB.

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% citep:colgan_2015_hybrid,bolch_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., citet:sutterley_2014,khan_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.

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 (call_tso(row=”VC”, col=”bias”) {{{results(-11)}}} Gt yr$^-1$) bias in comparison to citet:simonsen_2021 VC, there is a relatively high RMSE of call_tso(row=”VC”, col=”RMSE”) {{{results(86)}}} Gt yr$^-1$ and a mid-range correlation (r^2 = call_tso(row=”VC”, col=”r”) {{{results(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.

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 citep:nilsson_2015. citet:simonsen_2021 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 citep:simonsen_2021; 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. \ref{fig:GMB_VC_IMBIE}). 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.

IMBIE

The most widely cited estimate of Greenland mass balance today is the Ice-Sheet Mass Balance Inter-Comparison Exercise 2 (IMBIE2, citet:the-imbie-team_2019). 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 ~ call_tso(row=”IMBIE”, col=”bias”) {{{results(-7)}}} Gt 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 ~25 Gt yr^-1. Over the entire 26-year comparison period, the RMSE with IMBIE2 is call_tso(row=”IMBIE”, col=”RMSE”) {{{results(44)}}} Gt yr$^-1$ and the correlation is call_tso(row=”IMBIE”, col=”r”) {{{results(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 (c. ±30 to ±63 Gt 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.

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 BMB_VHD from RCM runoff are related, and D ice thickness and BMB_VHD 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 \ref{eq:MB} 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 citet:mankoff_2020_solid D data (up to 30 days old, with an error of ~9 % or ~45 Gt yr$^-1$) and the forecasted now-plus-7-day D (see Sect. \ref{sec:uncertainty:D}). Table \ref{tab:uncertainty} provides a summary of the uncertainty for each input.

The final This Study MB uncertainty value shown in Table \ref{tab:uncertainty} comes from the mean of the annual sum of the MB error term.

\begin{table}[ht] \centering \begin{tabular}{p{2.7cm} | p{2.2cm} | p{11cm}} Term & Uncertainty [±] & Notes \ \hline HIRHAM / \newline HARMONIE SMB & 15 \% & \citet{langen_2017}. The mean accumulation bias (-5\%) and ablation bias (-7\%) tend to cancel out, but this cannot be expected to be the case on single-basin, short-term scales where uncertainty is estimated to be larger.
MAR SMB & 15 \% & \citet{fettweis_2020}. The mean bias between the model and the measurements was 15 \% with a maximum of 1000 mmWE yr$^-1$. GrSMBMIP uses integrated values over several months of SMB, suggesting larger uncertainty of modeled runoff at the daily timescale. \ RACMO SMB & 15 \% & \citet{noel_2019}. Average 5\% runoff bias compared to annual cumulative discharge from the Watson River. Increases to a maximum of 20 \% for extreme runoff years. \ This Study SMB & 9 \% & Average of 15 \% SMB uncertainties above, assuming uncorrelated. \ Reconstructed SMB & $∼20$ \% & From \citet{kjeldsen_2015} Table 1. \ Recent D & $∼45$ Gt yr$^-1$ & \(∼9\) \%. \citet{mankoff_2020_solid} updated \citep{mankoff_2020_ice_github}. \ Reconstructed D & \(∼10\) \% & From \citet{kjeldsen_2015} Table 1. \ BMB_\mathrm{GF} & 50 \% & 5.3 +4/-1.4 Gt yr$^-1$ from \citet{karlsson_2021} Table 1, using the average of the three available methods. \ BMB_{\mathrm{friction}} & 20 \% & 11.8 $±$3.4 Gt yr$^-1$ from \citet{karlsson_2021} Table 1.\ BMB_\mathrm{VHD} & 15 \% & MAR runoff uncertainty.\ This Study MB & $∼86$ Gt yr$^-1$ & Eq \ref{eq:MB}, assuming all uncertainty is uncorrelated. \end{tabular} \caption{Summary 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].} \label{tab:uncertainty} \end{table}

Discharge

The D uncertainty is discussed in detail in citet:mankoff_2020_solid, 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, citet:mankoff_2020_A380 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 ~12 day observations upsampled to daily, and those ~12 day resolution observations come from longer time period observations citep:solgaard_2021. Because the velocity method uses feature tracking, it is correct on average but misses variability within each sample period (e.g., citet:greene_2020).

Spatially, discharge is estimated ~5 km 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.

ROIs

We work on the three different domains of the three RCMs and expand the ROIs to match the RCMs (see Appendix \ref{apdx:RCM_coverage}). However, some alignment issues cannot be solved. For example, we use BedMachine ice thickness to estimate BMB_VHD. Often, the largest BMB_VHD 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 BMB_VHD 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 km², of which 1,708,400 km² (src_octave{round(1708400/1825600*100)} {{{results(94)}}} %) are covered by the BedMachine ice mask and 26,400 km² (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 citet:kjeldsen_2020).

Accumulating uncertainties

When accumulating errors as in Fig. \ref{fig:mb_cumsum}, we use only the D and BMB_GF 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. \ref{fig:mb_cumsum}). There may be time-invariant biases in the BMB_friction 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. \ref{fig:mb_cumsum} 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. \ref{fig:mb_cumsum} representing the uncertainty for This Study is computed as 365 day rolling smoothly from 1840 through 1999 of the above-described uncertainty, 1/365th of the annual error at now + 7 days, and a linear blend, from 2000 to now + 7 days, between the smoothed reconstructed uncertainty and the present and future more variable uncertainty.

The citet:mouginot_2019, citet:colgan_2019, and citet:kjeldsen_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 citet:mouginot_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 citet:colgan_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/-91 Gt 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 citet:kjeldsen_2015 data have 0 surface mass balance and discharge bias by definition (Sect. \ref{sec:methods:reconstructed_adjustment}), but Fig. \ref{fig:mb_cumsum} displays the unadjusted data, and we apply a 36 Gt yr^-1 accumulating uncertainty from the unadjusted D bias (Fig. \ref{fig:reconstructed_adj}).

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 citep:bolch_2013. From 2003 to 2009 and using the VC method (altimetry), citet:gardner_2013 estimate -38 ±7 Gt yr^-1 peripheral mass balance. From 2006 to 2016 and using the VC method (DEM differencing), citet:zemp_2020 estimate -51 ±17 Gt yr^-1 peripheral mass balance using citet:rastner_2012 delineations.

Results

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

At the decadal average, the following trends are apparent. Surface mass balance has decreased from a high of ~450 Gt yr^-1 in the 1860s to a low of ~260 Gt yr^-1 in the 2010s. SMB variability has also increased during this time. Discharge has increased slightly from a low of ~375 Gt yr^-1 in the 1860s to a high of ~490 Gt yr^-1 in the 2010s. Basal mass balance, from runoff as a proxy, had a high of 26 ±16 Gt yr^-1 in the 1930s and a low of 22 ±5 Gt 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 ~-10 Gt yr^-1, has only 1 (of 19) decades with a mass gain (~50 Gt yr^-1 in the 1860s), and has a record low of ~-250 Gt yr^-1 in the 2010s.

Data availability

The RCM surface mass balance and the VHD basal mass balance components are updated daily, the discharge approximately every 12 days, and all are used to produce the final daily-updating product. The data are available at https://doi.org/10.22008/FK2/OHI23Z citep:this_study_data, 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 (https://github.com/GEUS-Glaciology-and-Climate/mass_balance/; last visited \today; citet:this_study_github) where users can track progress, make suggestions, and discuss, report, and respond to issues that arise during use of this product.

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

RCM differences

Mouginot 2019 by region

Colgan 2019 by sector

Reconstructed runoff

RCM coverage

Software

This work was performed using only open-source software, primarily GRASS GIS citep:neteler_2012, CDO citep:CDO, NCO citep:NCO, GDAL citep:GDAL, and Python citep:python, in particular the Jupyter citep:kluyver_2016, dask citep:dask_sw,dask_paper, pandas citep:pandas, geopandas citep:geopandas, numpy citep:numpy, x-array citep:xarray, and Matplotlib citep:matplotlib packages. The entire work was performed in Emacs citep:emacs using Org Mode citep:schulte_2012 on GNU/Linux and using many GNU utilities. The parallel citep:parallel tool was used to speed up processing.

CRediT

\begin{figure}[!h] \centering \includegraphics[width=0.5\textwidth]{./fig/credit.png} \caption{\label{fig:credit}Author contributions following the CRediT system \citep{allen_2014,brand_2015,allen_2019}} \end{figure}

Misc journal sections

\authorcontribution{Author contribution is captured following the CRediT system \citep{allen_2014,brand_2015,allen_2019} and shown graphically in Figure \ref{fig:credit}. 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.}

\competinginterests{The authors declare that they have no conflict of interest.}

References