Greenland ice sheet mass balance from 1840 through next week

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, MAR, and RACMO) from 1986 through next week. Additional non-SMB losses come from a marine 5 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 ~24 Gt 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 10 other IO mass balance estimates. Results from this study are available at https://doi.org/10.22008/FK2/OHI23Z (Mankoff et al., 2021).


Introduction
Over the past several decades, mass loss from the Greenland ice sheet has increased The IMBIE Team, 2019). Different processes dominate the regional mass loss of the ice sheet, and their relative contribution has fluctuated in 15 time . For example, in the 1970s nearly all sectors gained mass due to positive SMB, except the northwest 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).
There are three common methods for estimating mass balance -changes in gravity (Barletta et al., 2013;Groh et al., 2019; The IMBIE Team, 2019; Velicogna et al., 2020), changes in volume (Simonsen et al., 2021a;Sørensen et al., 2011;Zwally and 20 Giovinetto, 2011;Sasgen et al., 2012;Smith et al., 2020), and the input/output (IO) method (Colgan 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~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 on 25 annual or multi-year temporal resolution (when), and with 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 30 Sentinel missions . 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 ) -a term implicitly included in the GMB and VC methods but neglected by all previous IO estimates. 35 In this work we introduce the new PROMICE Greenland ice sheet mass balance dataset 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, and condensation/riming) and negative SMB terms (surface melt, evaporation, sublimation, and snow drift erosion/deposition). We also use the basal mass balance, and an estimate of dynamic ice discharge. Spatial 40 resolution is effectively per sector (Zwally et al., 2012) or region . 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 marine mass balance (i.e., 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 dailyupdating estimate of Greenland mass changes from 1840 through next week. 45

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 region and sector scale -Reconstructed refers to the adjusted Kjeldsen et al. (2015) annual temporal resolution data at ice sheet scale used to 50 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, then discarded. for each of the seven  regions. The map shows both the named regions  and the numbered sectors (Zwally et al., 2012). Marine mass balance 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 2021-08-20 based on the date this document was compiled.
The inputs to this work are the recent SMB fields from the three RCMs, the recent marine mass balance data from Mankoff et al. (2020b) (data: Mankoff and Solgaard (2020)), and the recent basal mass balance fields, of which BMB GF  are direct outputs from Karlsson et al. (2021) (data: ), but the BMB VHD calculations are redone here (see Methods Sect. 5.3) using the MAR runoff field. The reconstructed data (pre 1986) are surface and marine mass balance from 85 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 BMB VHD (see Methods Sect. 5.3).

Surface mass balance
We use one reconstructed SMB from 1840 through 1985, and three recent SMB 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 (Langen et al., 2017), which is forced with surface fluxes of energy (turbu- The HIRHAM5 regional climate model (Christensen et al., 2007) combines the dynamical core of the HIRLAM7 numerical weather forecasting model (Eerola, 2006) 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°x 0.05°on a 100 rotated pole grid (corresponding to 5.5 km resolution), and 31 atmospheric levels in the vertical. It is forced at 6 hr 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).

105
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 Oerlemans, 1996;Lefebre, 110 2003). The offline subsurface model is run on the HIRHAM5 5.5 km grid.
For the real-time data we use DMI's operational numerical weather forecast model HARMONIE (Bengtsson et al., 2017), a nonhydrostatic model in terrain-following sigma coordinates based on the fully compressible Euler equations (Simmons and Burridge, 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, radio occultation data, 115 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 Randriamampianina, 2021;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 the offline model both 120 employ the Citterio and Ahlstrøm (2013) ice mask interpolated to the 5.5 km grid.

MAR
The Modèle Atmosphérique Régional (MAR) RCM has been developed by the University of Liège (Belgium) with a focus on the polar regions . The MAR atmosphere module (Gallée and Schayes, 1994) is fully coupled with the soil-ice-snow energy balance vegetation model SISVAT (Gallée et al., 2001) simulating the evolution of the 30 first meters of 125 snow/ice over the ice sheet with the help of 30 snow layers (with time varying thickness) or the 10 first meter of soil over the tundra area. At its lateral boundary, MAR is 6 hourly forced by the reanalysis ERA5 and runs at a resolution of 20 km. The snow pack has been initialised 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,

130
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).
The new version, MAR 3.12, is used here. With respect to version 3.9 intensively validated over Greenland  or the 20 km based MARv3.10 set-up 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 assures 135 now the full conservation of water mass into both soil and snowpack at each time step, takes into account of 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, 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 In addition to providing SMB, MAR also provides daily runoff over both permanent ice and tundra area -this is used for the daily BMB VHD estimate (Section 5.3).
As the recent SMB decrease (successfully evaluated with GRACE based estimates in ) has been fully driven by the increase in runoff , we assume the same degree of accuracy between SMB simulated by 145 MAR (evaluated with the PROMICE SMB database ) and the runoff simulated by MAR. 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 , and here we repeat the main characteristics.

160
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 in 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 are initialized in September 1957 using temperature and density from a previous run with the offline IMAU Firn Densification

170
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 down-sampled from the 90 m resolution Greenland Ice Mapping Project (GIMP) digital elevation model (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 1958and 1978(Uppala et al., 2005,

175
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 (~130 km) wide to ensure a smooth transition to the domain interior. This run has active upper atmosphere relaxation (van de Berg and Medley, 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.

180
The 5.5 km product is therefore statistically downscaled onto a 1 km grid sampled from the GIMP DEM , employing corrections for biases in elevation and bare ice albedo using a MODIS albedo product at 1 km resolution (Noël et al., 2016).

Reconstructed
The Kjeldsen et al. (2015) 173-year (1840 through 2012) mass balance reconstruction is based on the Box (2013) 171-year 185 (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; 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;Cappelen, 190 2001; Cappelen et al., 2006;Vinther et al., 2006)) or firn/ice cores  and b) spatially continuous outputs from regional climate model RACMO version 2.1 (Ettema et al., 2010). A 43-year overlap period (1960 through 2012) with the RACMO data are 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 hydrologicalyear cumulative SMB. A fundamental assumption is that the calibration factors, regression slope, and offset for the calibration The reconstructed surface mass balance is adjusted as described in the Methods Sect. 5.4 (Fig. 3).

Marine mass balance
The recent marine mass balance data are the discharge ( The reconstructed marine mass balance data  are estimated via a linear fit between unsmoothed annual marine mass balance spanning 2000 to 2012 (Enderlin et al., 2014) and runoff data from  using a 6-year trailing average. The method for scaling discharge from runoff was introduced by (Rignot et al., 2008), who scaled the SMB discharge predictor, and include a discussion of the physical basis. Although the fitting period of the present dataset includes an anomalous period of discharge (2000 through 2005; e.g., Boers and Rypdal (2021)), the discharge data used by Rignot et al. (2008) and Box and Colgan (2013)  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 Bjørk et al., 2012;Khan et al., 2015Khan et al., , 2020.
The reconstructed marine mass balance is adjusted as described in the Methods Sect. 5.4.

Basal mass balance 225
The basal mass balance (BMB; Karlsson et al. (2021)) comes from mass lost at the bed from geothermal flux (BMB GF ), frictional heating (BMB friction ) from the basal shear velocity, and viscous heat dissipation (BMB VHD ) from surface runoff routed to the bed (i.e. the volume of the subglacial conduits formed from surface runoff; Mankoff and Tulaczyk (2017)).
These fields (data: Karlsson (2021)) 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 230 appropriate for this work with daily resolution. We therefore re-calculate the BMB VHD -induced basal melt as described in Methods Sect. 5.3.

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., Rogozhina et al. (2012); where E GF is available energy at the bed, here the geothermal flux in units W m -2 , ρ i is the density of ice (917 kg m -3 ), and 240 L is the latent heat of fusion (335 kJ kg -1 ; Cuffey and Paterson (2010)). BMB GF 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.

Friction
This heat term stems from the friction produced as ice slides over the bedrock. The term has only been measured in a handful 245 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, similar to Eq. 1: 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 (MacGregor et al., 2016) in order to mask out areas that are likely frozen.

Other
ROI regions come from  and ROI sectors come from Zwally et al. (2012).

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 , the previous PROMICE mass balance product (Colgan et al. (2019); data: Colgan (2021)), 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 (The IMBIE Team, 2019). In addition to 260 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.

Methods
The total mass balance for all of Greenland and all the different ROIs involves summing each field (SMB, MMB, BMB) by each ROI, then subtracting the MMB and BMB from the SMB fields, or, have total mass balance defined as and when comparing This Study to those products, we compare like terms, never comparing our MB to a different product 270 MB * , except in Fig. 4 where all products are shown together.
Prior to calculating the mass balance, we perform the following steps.

Projected marine mass balance
We project the marine mass balance from the last observed point from Mankoff et al. (2020b) (generally between 2 weeks and 1 month old) to seven days into the future at each glacier. We define the long-term trend as the linear least squares fit to the last three years of observed data. We define the seasonal signal as the mean of the last three years of observed data during 285 the temporal window of interest that spans from the most recently available observation through to present plus one week. We then assign the long term trend plus the seasonal signal value to next week's MMB (i.e., now + 7 days), and finally linearly interpolate between the last observed value and the forecasted next week's value, to provide the recent past projected and future forecasted MMB.
Marine mass balance does not change sign and changes magnitude by approximately 6 % annually over the entire ice sheet 290 (King et al., 2018), but surface mass balance changes sign and has much larger and short-term variability. From this, the statistical forecast for marine mass balance described above does not impact results as much as the physically-based model forecast for surface mass balance.

Basal mass balance
Because Karlsson et al. (2021) provide a steady-state annual-average estimate of the BMB fields, we divide the BMB GF and 295 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 (King et al., 2018), and Karlsson et al. (2021) found that basal melt rates are 5 % higher for summer maps. 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 300 for this study.
To estimate recent BMB VHD we use daily MAR runoff (see Mankoff et al. (2020a)) and BedMachine v4 (Morlighem et al., 2017(Morlighem et al., , 2021 to derive subglacial routing pathways, similar 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, 305 h, defined as 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. 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). 310 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øl, 1956;Mankoff and Tulaczyk, 2017). We assume all energy, E VHD (in Joules), is used to melt ice with

315
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 320 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øl, 1956). To estimate the scale factor we use the 1986 through 2012 overlap between Kjeldsen et al. (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, 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 325 described earlier) to estimate reconstructed BMB VHD .

Reconstructed adjustment
We use the reconstructed and recent surface (SMB) and marine (MMB) mass balance 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 165 Gt yr -1 (Fig. 3 SMB), and This Study vs reconstructed MMB has a slope of 1.1 and an intercept of -6 Gt yr -1 (Fig. 3 MMB). The unadjusted reconstructed We adjust the reconstructed data until the reconstructed vs. recent slope is 1 and intercept is 0 Gt yr -1 for each of the surface and marine mass balance comparisons (Fig. 3). We then derive the BMB VHD term for reconstructed basal mass balance (Sect. For reconstructed SMB and MMB, 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 MMB uncertainty.

340
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 marine mass balance, years with low MMB are slightly reduced, and years with high MMB 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 345 adjustment to the demonstrated too-high biases in accumulation and ablation estimates in the 1840-2012 reconstructed SMB field , an offset resulting from differences in ice masks , the inclusion of peripheral glaciers , other accumulation rate inaccuracies (Lewis et al., 2017(Lewis et al., , 2019, and other unknowns. 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 365 sheet) are discarded in these locations.

Domains, boundaries, and regions of interest
We compare to six related data sets (see Table 2 and Sect. 4.5): The most similar and recent IO product , 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 (The IMBIE Team, 2019), and the unadjusted reconstructed/recent overlap 370 .
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 seven days from now (now defined as 2021-08-13). 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 375 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.  Because MB * is a linear combination of SMB and MMB terms (Eq 4), the MB * differences between this product and  are dominated by the MMB term, although it is not apparent because interannual variability is dominated 390 by SMB.

Colgan (2019)
The dominated by MMB disagreement, although it is again not apparent because interannual variability is dominated by SMB.

Gravimetric Mass Balance (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% 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 420 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 425 applied to the gravimetric signal may also lead to differences in GMB estimates on ice sheet scale, but also on 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 430 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. We suggest that this is due to climate influences on the effective radar horizon across the ice sheet during these years.

445
Weather-driven changes in the effective scatter horizon, mapped by 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 450 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 455 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, 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 460 methodological differences between This Study and IMBIE2. Firstly, the gravimetry members of IMBIE2 assess mass balance of all Greenlandic land ice, including peripheral ice masses, while This Study only assesses 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: 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, marine, and basal mass balance, or SMB, MMB, 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 MMB ice thickness and BMB VHD pressure gradients are related, and other terms may have dependencies. However, the two dominant IO terms, SMB inputs and MMB outputs, are independent on annual time scales, and for simplification we 480 treat all terms as independent. We use Eq 3 and standard error propagation for SMB, MMB, and BMB terms (i.e., the square root of the sum of the squares of the SMB plus MMB plus BMB error terms). For the MMB, extra work is done to calculate uncertainty between the last Mankoff 2020 MMB data (up to 30 days old, with error of~9 % or~45 Gt yr −1 ) and the forecasted now-plus-7-day MMB (see Sect. 7.1). Table 3 provides a summary of the uncertainty for each input.
The final This Study MB uncertainty value shown in Table 3 comes from mean of the annual sum of the MB error term.  Langen et al. (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. Eq 3, assuming all uncertainty is uncorrelated. Table 3. 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 ].

Marine mass balance
The MMB 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 %.

490
Temporally, MMB at daily resolution comes from~12 day observations up-sampled to daily, and those~12 day resolution observations come from longer time period observations . Because the velocity method using feature tracking, it is correct on average but misses variability within each sample period (e.g., Greene et al. (2020)).
Spatially, MMB 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 steady at 1 km yr -1 would take 5 years before that mass is lost. However, 495 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 500 sheet margin where there are both mask issues and high topographic variability.
The forecasted MMB uncertainty is estimated from the same forecasted calendar days of the last three years of observed MMB at each glacier. MMB uncertainty for the first forecasted day is the baseline MMB uncertainty plus the largest daily change observed during the forecast period from the last three years. The MMB uncertainty for the second forecasted day is the baseline uncertainty plus the uncertainty from the first forecasted day plus the second largest daily change observed 505 during the forecast period from the last three years. We repeat this for the approximately 30 days of forecasted MMB. This implementation takes into account the larger variability (uncertainty) during the seasonal transition between the lower winter and higher summer discharge, or the smaller variability during the winter period.

Regions of interest (ROI)
We work on the three different domains of the three RCMs, and expand the ROIs to match the RCMs (see Appendix E).

510
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, if 515 these RCM ice grid cells with high runoff are outside the BedMachine ice domain (ice thickness is 0), there is no reasonable way to include that runoff in our BMB VHD budget, and these grid cells are ignored.
The MAR ice domain is 1,828,800 km 2 of which 1,711,200 km 2 (94 %) are covered by the BedMachine ice mask, and 26,000 km 2 (6 %) are not. This 6 % area contributes~15 % of the runoff (excluded from the VHD calculations) and likely a higher percent of the VHD, because the border region of the ice sheet has the steepest gradients and the largest volume of 520 subglacial flow.
We encourage RCM developers, BedMachine, and others to use a common and up-to-date mask (see Kjeldsen et al. (2020)).

Accumulating uncertainties
When accumulating errors as in Fig. 4, we use only the MMB and BMB GF uncertainty. The MMB 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 525 have errors randomly distributed in time (for the purposes of Fig. 4). 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. 4 where errors accumulate in time. The provided data product contains one uncertainty field and does not distinguish between systematic and random uncertainty.

20
The shaded region in Fig. 4 representing the uncertainty for This Study is computed as a 365 day rolling smooth from 1840 530 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 , 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  data the same as This Study. Marine mass balance uncertainty is treated as a bias and 535 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 -35 and 61 Gt yr -1 respectively. This suggests that in any given year, there could be up to -35 ± 61 or +26/-96 Gt yr -1 departure from This Study. From this, we assign a 35 Gt yr -1 bias (36 %; accumulates in time) and a 61 Gt yr -1 RMSE (64 %; random in time).
The adjusted Kjeldsen et al. (2015) data have 0 surface and marine mass balance bias by definition (Sect. 5.4), but Fig. 4 540 displays the unadjusted data, and we apply a 36 Gt yr -1 accumulating uncertainty from the unadjusted MMB bias (Fig. 3).

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 sheetthey are < 5 % of the total ice area but contribute~15 to 20 % of the whole island mass loss .
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 low of~260 Gt yr -1 in the 2010s. SMB variability has also increased during this time. Marine mass balance 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 555 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, is increasing in recent decades.
The total mass balance decadal trend from the 1840s through the 2010s is one of general mass decrease and increased intradecadal variability. The record begins in the 1840s with~-10 Gt yr -1 , has only one (of 19) decades with a mass gain (~50 Gt yr -1 in the 1860s), and a record low of~-250 Gt yr -1 in the 2010s. The RCM surface mass balance, and the VHD basal mass balance components are updated daily, the marine mass balance approximately every 12 days, and all are used to produce the final daily-updating product. The data area available at https: //doi.org/10.22008/FK2/YG3IWC , with all historical (daily updated) versions archived.
As part of our commitment to make continual and improving updates to the data product, we introduce a GitHub database 565 (https://github.com/GEUS-Glaciology-and-Climate/mass_balance/; last visited August 13, 2021) where users can track progress, make suggestions, discuss, report and respond to issues that arise during use of this product.

Conclusions
This study is the first to provide a dataset 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, marine 570 mass balance, 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 twenty 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 on 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 575 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 580 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 resolution 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 585 presented will also be able to improve, as it is sustained by the Danish-Greenlandic governmental long-term monitoring effort -the Programme for Monitoring of the Greenland ice sheet (PROMICE).  . All axes units are Gt yr -1 . Plotted numbers represent the last two digits of the years for the unadjusted data sets. The matching colored diamonds 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 (diamonds) the surface and marine mass balance adjustments.  Table 2). This Study annual-resolution data prior to 1986 is 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.    33 Appendix E: RCM coverage Figure E1. HIRHAM RCM coverage by . Coverage of HIRHAM by Zwally et al. (2012), and MAR and RACMO by  and Zwally et al. (2012) is similar to graphic shown here (See section 5.5 for discussion of RACMO coverage issues). HIRHAM latitude and longitude covers the equator because we work on the native HIRHAM rotated pole coordinate system.