Greenland Ice Sheet solid ice discharge from 1986 through March 2020

Abstract

We present a 1986 through March 2020 estimate of Greenland Ice Sheet ice discharge. Our data include all discharging ice that flows faster than 100 m yr${}^{-\mathrm{1}}$ and are generated through an automatic and adaptable method, as opposed to conventional handpicked gates. We position gates near the present-year termini and estimate problematic bed topography (ice thickness) values where necessary. In addition to using annual time-varying ice thickness, our time series uses velocity maps that begin with sparse spatial and temporal coverage and end with near-complete spatial coverage and 12 d updates to velocity. The 2010 through 2019 average ice discharge through the flux gates is $\sim \mathrm{487}\pm \mathrm{49}$ Gt yr${}^{-\mathrm{1}}$. The 10 % uncertainty stems primarily from uncertain ice bed location (ice thickness). We attribute the $\sim \mathrm{50}$ Gt yr${}^{-\mathrm{1}}$ differences among our results and previous studies to our use of updated bed topography from BedMachine v3. Discharge is approximately steady from 1986 to 2000, increases sharply from 2000 to 2005, and then is approximately steady again. However, regional and glacier variability is more pronounced, with recent decreases at most major glaciers and in all but one region offset by increases in the northwest region through 2017 and in the southeast from 2017 through March 2020. As part of the journal's living archive option and our goal to make an operational product, all input data, code, and results from this study will be updated as needed (when new input data are available, as new features are added, or to fix bugs) and made available at https://doi.org/10.22008/promice/data/ice_discharge and at https://github.com/mankoff/ice_discharge (last access: 6 June 2020, ).

1 What is new in this edition

The data have been updated repeatedly between the first edition and this edition. The data will continue to be updated (currently every 12 d), although reference papers will only come out when major changes occur in the processing algorithms or input data. Therefore users are encouraged to regularly check for data updates at 10.22008/promice/data/ice_discharge when using the data.

A post-peer-review website is available at https://github.com/mankoff/ice_discharge (last access: 6 June 2020, ), where we document changes to the paper and use the GitHub Issues feature to collect suggested improvements to the paper, document those improvements as they are implemented, document problems that made it through review, and mention similar papers that have been published since this was accepted. The git commit for this edition is 8acd36f.

In this edition the NSIDC 0478 ice velocity data have been updated from v2 to v2.1. These data are used for the baseline velocity and gate selection. Therefore the gate locations, number of gates, and number of gate pixels have all changed. Overall there are $\sim \mathrm{3}$ % fewer pixels, but only eight fewer gates, and some of the gate reduction is due to combining two gates on the same glacier that had a 1-pixel gap in the previous edition. The effect of these changes is $<\mathrm{2.5}$ % of the estimated discharge. In this edition there are now 5829 pixels and 267 gates. Many gates remain the same pixel for pixel, but some do not meet the same inclusion criteria applied to these new baseline data, while some glaciers are included in this edition that were not in the previous edition. Many gates are also a few pixels narrower.

The NSIDC 0646 ice velocity data have been updated from v2 to v2.1. This update increases coverage and discharge in the 1980s by $\sim \mathrm{25}$ to 40 Gt yr${}^{-\mathrm{1}}$ ($\sim \mathrm{6}$ % to 10 %) due to higher velocity estimates than the previous product that covered that time period with annual averages . This change highlights that ice-sheet-wide differences between velocity products can be nontrivial (see ). The time series has also been extended through both the updated NSIDC 0646 data and 61 Sentinel-1 velocity maps from 2018 through the present (29 March 2020). We have also added 48 additional MEaSUREs (; hereafter NSIDC 0731) monthly average velocity maps from 1 December 2014 through 30 November 2018.

We have updated the time series graphics (Figs. , , and ) in the following manner: any observation (gate-, region-, or ice-sheet-wide) where coverage is $<\mathrm{50}$ % is discarded from the graphic (low-coverage data are still included in the downloadable data), and annual average is only computed if there are three or more samples in a year.

Finally, the Supplement includes significantly more meta-data about the input data used in this work to aid in both reproducibility by third parties and in tracking the impact of additional and updated input data on future versions of this work.

Results show a continued steady total discharge. The contributions from the central west (CW) region continue to decrease, while the central east (CE) region continues to increase, and CE and CW are now approximately tied for the third-largest discharging region. The top three individual contributing glaciers remain dynamic – Sermeq Kujalleq (English: Jakobshavn Glacier; Danish: Jakobshavn Isbræ) continued its rapid discharge decline in 2017 and 2018, returning to approximately its discharge from year 2000, but increased discharge slightly in 2019. For some time in 2018 and all data points so far in 2020 (i.e., though March) Helheim was the top Greenlandic glacier contributing to sea-level rise, although not with statistical significance (error bars still overlap).

2 Introduction

The mass of the Greenland Ice Sheet is decreasing (e.g., ). Most ice sheet mass loss – as iceberg discharge, submarine melting, and meltwater runoff – enters the fjords and coastal seas, and therefore ice sheet mass loss directly contributes to sea-level rise . Greenland's total ice loss can be estimated through a variety of independent methods, for example direct mass change estimates from GRACE , or by using satellite altimetry to estimate surface elevation change, which is then converted into mass change (using a firn model, e.g., ). However, partitioning the mass loss between ice discharge ($D$) and surface mass balance (SMB) remains challenging (see ). Correctly assessing mass loss, as well as the attribution of this loss (SMB or $D$), is critical to understanding the process-level response of the Greenland Ice Sheet to climate change and thus improving models of future ice sheet changes and associated sea-level rise .

The total mass of an ice sheet, or a drainage basin, changes if the mass gain (SMB inputs, primarily snowfall) is not balanced by the mass loss ($D$ and SMB outputs, the latter generally meltwater runoff). This change is typically termed ice sheet mass balance (MB), and the formal expression for this rate of change in mass is (e.g., )

1$${\displaystyle \frac{\mathrm{d}M}{\mathrm{d}t}}=\mathit{\rho }\underset{A}{\int }b\mathrm{d}A-\underset{g}{\int }Q\mathrm{d}g,$$ where $\mathit{\rho }$ is the average density of ice, $b$ is an area mass balance, and $Q$ is the discharge flux. The left-hand side of the equation is the rate of change of mass, the first term on the right-hand side is the area $A$ integrated surface mass balance (SMB), and the second term is the discharge $D$ mass flow rate that drains through gate $g$. Equation () is often simplified to 2$$\text{MB}=\text{SMB}-D,$$ where MB is the mass balance, and referred to as the “input–output” method (e.g., ). Virtually all studies agree on the trend of Greenland mass balance, but large discrepancies persist in both the magnitude and attribution. Magnitude discrepancies include, for example, reporting a mass imbalance of $-\mathrm{250}\pm \mathrm{21}$ Gt yr${}^{-\mathrm{1}}$ during 2003 to 2010, reporting $-\mathrm{181}\pm \mathrm{28}$ Gt yr${}^{-\mathrm{1}}$ during 2003 to 2008, and reporting a mass imbalance of $-\mathrm{265}\pm \mathrm{19}$ Gt yr${}^{-\mathrm{1}}$ during 2004 to 2008. Some of these differences may be due to different ice sheet area masks used in the studies. Attribution discrepancies include, for example, attributing the majority (64 %) of mass loss to changes in SMB during the 2005 to 2009 period but attributing the majority (85 %) of mass loss to changes in $D$ during the 2004 to 2008 period.

Discharge may be calculated through several methods, including mass flow rate through gates (e.g., ), or solving as a residual from independent mass balance terms (e.g., ). The gate method that we use in this study incorporates ice thickness and an estimated vertical profile from the observed surface velocity to calculate the discharge. A typical formulation of discharge across a gate $D_{\mathrm{g}}$ is 3$$D_{\mathrm{g}}=\mathit{\rho }VHw,$$ where $\mathit{\rho }$ is the average density of ice, $V$ is depth-average gate-perpendicular velocity, $H$ is the ice thickness, and $w$ is the gate width. Uncertainties in $V$ and $H$ naturally influence the estimated discharge. At fast-flowing outlet glaciers, $V$ is typically assumed to be equal at all ice depths, and observed surface velocities can be directly translated into depth-averaged velocities (as in ). To minimize uncertainty from SMB or basal mass balance corrections downstream of a flux gate, the gate should be at the grounding line of the outlet glacier. Unfortunately, uncertainty in bed elevation (translating to ice thickness uncertainty) increases toward the grounding line.

Conventional methods of gate selection involve handpicking gate locations, generally as linear features (e.g., ), or visually approximating ice-orthogonal gates at one point in time (e.g., ). Manual gate definition is suboptimal. For example, the largest discharging glaciers draw from an upstream radially diffusing region that may not easily be represented by a single linear gate. Approximately flow-orthogonal curved gates may not be flow orthogonal on the multidecade timescale due to changing flow directions. Manual gate selection makes it difficult to update gate locations, corresponding with glacier termini retreat or advance, in a systematic and reproducible fashion. We therefore adopt an algorithmic approach to generate gates based on a range of criteria.

Here, we present a discharge dataset based on gates selected in a reproducible fashion by a new algorithm. Relative to previous studies, we employ ice velocity observation over a longer period with higher temporal frequency and denser spatial coverage. We use ice velocity from 1986 through 2019, including 12 d velocities for the last $\sim \mathrm{500}$ d of the time series, and discharge at 200 m pixel resolution capturing all ice flowing faster than 100 m yr${}^{-\mathrm{1}}$ that crosses glacier termini into fjords.

3 Input data

Historically, discharge gates were selected along well-constrained flight lines of airborne radar data . Recent advances in ice thickness estimates through NASA Operation IceBridge , NASA Oceans Melting Greenland (OMG; ), fjord bathymetry , and methods to estimate thickness from surface properties (e.g., ) have been combined into digital bed elevation models such as BedMachine v3 or released as independent datasets . From these advances, digital bed elevation models have become more robust at tidewater glacier termini and grounding lines. The incorporation of flight-line ice thickness data into higher-level products that include additional methods and data means gates are no longer limited to flight lines (e.g., ).

Ice velocity data are available with increasing spatial and temporal resolution (e.g., ). Until recently, ice velocity mosaics were limited to once per year during winter , and they are still temporally limited, often to annual resolution, prior to 2000 (e.g., ). Focusing on recent times, ice-sheet-wide velocity mosaics from the Sentinel-1A and 1B are now available every 12 d (http://PROMICE.org, last access: 6 June 2020). The increased availability of satellite data has improved ice velocity maps both spatially and temporally, thereby decreasing the need to rely on spatial and temporal interpolation of velocities from annual/winter mosaics .

The discharge gates in this study are generated using only surface speed and an ice mask. We use the MEaSUREs Greenland Ice Sheet Velocity Map from InSAR Data, Version 2 , hereafter termed “MEaSUREs 0478” due to the National Snow and Ice Data Center (NSIDC) dataset ID number. We use the BedMachine v3 ice mask.

For ice thickness estimates, we use surface elevation from GIMP (; NSIDC dataset ID 0715), adjusted through time with surface elevation change from and bed elevations from BedMachine v3 replaced by where available. Ice sector and region delineation is from . Ice velocity data are obtained from a variety of products including Sentinel-1A and 1B derived by PROMICE (see Appendix ), MEaSUREs 0478, MEaSUREs 0646 , , and . Official glacier names come from . Other glacier names come from . See Table  for an overview of datasets used in this work.

Table 1

Summary of data sources used in this work.

Property	Name used in this paper	Reference
Basal topography	BedMachine
Basal topography for southeast
Surface elevation	GIMP 0715
Surface elevation change	Surface elevation change	,
Baseline velocity	MEaSUREs 0478
Velocity	Sentinel	Appendix
Velocity	MEaSUREs 0646
Velocity	MEaSUREs 0731	,
Velocity	pre-2000
Sectors and regions	Sectors and regions
Names		,

This work uses 367 different velocity maps, biased toward post-2015 when 12 d ice velocities become available from the Sentinel-1 satellites. The temporal distribution is $\sim \mathrm{10}$ maps per year from 1986 to 2013, 14 in 2014, 25 in 2015, 36 in 2016, 69 in 2017, 42 in 2018, and 24 in 2019.

4 Methods

4.1 Terminology

We use the following terminology, most displayed in Fig. :

• “Pixels” are individual $\mathrm{200}\mathrm{m}\times \mathrm{200}\mathrm{m}$ raster discharge grid cells. We use the nearest neighbor when combining datasets that have different grid properties.

• “Gates” are contiguous (including diagonal) clusters of pixels.

• “Sectors” are spatial areas that have $\mathrm{0}$, $\mathrm{1}$, or $>\mathrm{1}$ gate(s) plus any upstream source of ice that flows through the gate(s) and come from .

• “Regions” are groups of sectors, also from , and are labeled by approximate geographic region.

• The “baseline” period is the average 2015, 2016, and 2017 winter velocity from MEaSUREs 0478.

• “Coverage” is the percentage of total, region, sector, or gate discharge observed at any given time. By definition coverage is 100 % during the baseline period. From the baseline data, the contribution to total discharge of each pixel is calculated, and coverage is reported for all other maps that have missing observations (Fig. ). Total estimated discharge is always reported because missing pixels are gap filled (see “Missing or invalid data” section below).

• “Fast-flowing ice” is defined as ice that flows more than 100 m yr${}^{-\mathrm{1}}$.

• Names are reported using the official Greenlandic names from ; if an alternate name exists e.g., fromor an English version, then this is shown in parentheses.

4.2 Gate location

Gates are algorithmically generated for fast-flowing ice (greater than 100 m yr${}^{-\mathrm{1}}$) close to the ice sheet terminus determined by the baseline-period data. We apply a 2D inclusive mask to the baseline data for all ice flowing faster than 100 m yr${}^{-\mathrm{1}}$. We then select the mask edge where it is near the BedMachine ice mask (not including ice shelves), which effectively provides grounding line termini. We buffer the termini 5000 m in all directions, creating ovals around the termini and once again down-select to fast-flowing ice pixels. This procedure results in gates 5000 m upstream from the baseline terminus that bisect the baseline fast-flowing ice. We manually mask some land- or lake-terminating glaciers which are initially selected by the algorithm due to fast flow and mask issues.

We select a 100 m yr${}^{-\mathrm{1}}$ speed cutoff because slower ice, taking longer to reach the terminus, is more influenced by SMB. For the influence of this threshold on our results see the Discussion section and Fig. .

Figure 2

Heatmap and table showing ice sheet discharge as a function of gate buffer distance and ice speed cutoff. The colors of the numbers change for readability.

[Figure omitted. See PDF]

We select gates at 5000 m upstream from the baseline termini, which means that gates are likely $>\mathrm{5000}$ m from the termini further back in the historical record . The choice of a 5000 m buffer follows from the fact that it is near terminus and thus avoids the need for (minor) SMB corrections downstream, yet it is not too close to the terminus where discharge results are sensitive to the choice of distance-to-terminus value (Fig. ), which may be indicative of bed (ice thickness) errors.

4.3 Thickness

We derive thickness from surface and bed elevation. We use GIMP 0715 surface elevations in all locations, and the BedMachine bed elevations in most locations, except southeast Greenland where we use the bed. The GIMP 0715 surface elevations are all time stamped per pixel. We adjust the surface through time by linearly interpolating elevation changes from , which covers the period from 1995 to 2016. We use the average of the first and last 3 years for earlier and later times, respectively. Finally, from the fixed bed and temporally varying surface, we calculate the time-dependent ice thickness at each gate pixel.

4.4 Missing or invalid data

The baseline data provide velocity at all gate locations by definition, but individual nonbaseline velocity maps often have missing or invalid data. Also, thickness provided by BedMachine is clearly incorrect in some places (e.g., fast-flowing ice that is 10 m thick, Fig. ). We define invalid data and fill in missing data as described below.

4.4.1 Invalid velocity

We flag invalid (outlier) velocities by treating each pixel as an individual time series, applying a 30-point rolling window, flagging values more than 2 standard deviations outside the mean, and repeating this filter three times. We also drop the 1972 to 1985 years from because there is low coverage and extremely high variability when using our algorithm.

This outlier detection method appears to correctly flag outliers (see , for unfiltered time series graphs) but likely also flags some true short-term velocity increases. The effect of this filter is a $\sim \mathrm{1}$ % reduction in discharge most years but more in years with high discharge – a reduction of 3.2 % in 2013, 4.3 % in 2003, and more in the 1980s when the data are noisy. Any analysis using these data and focusing on individual glaciers or short-term changes (or lack thereof) should reevaluate the upstream data sources.

4.4.2 Missing velocity

We generate an ice speed time series by assigning the PROMICE, MEaSUREs 0478, MEaSUREs 0646, and pre-2000 products to their respective reported time stamps (even though these are time-span products) or to the middle of their time span when they cover a long period such as the annual maps from . We ignore that any individual velocity map or pixel has a time span and not a time stamp. Velocities are sampled only where there are gate pixels. Missing pixel velocities are linearly interpolated in time, except for missing data at the beginning of the time series which are back- and forward filled with the temporally nearest value for that pixel (Fig. ). We do not spatially interpolate missing velocities because the spatial changes around a missing data point are most likely larger than the temporal changes. We visually represent the discharge contribution of directly observed pixels, termed coverage (Fig. ), as time series graphs and opacity of dots and error bars in the figures. The figures only display data where coverage is $\ge \mathrm{50}$ %, but the provided data files include coverage from 0 % to 100 %. Therefore, the gap-filled discharge contribution at any given time is equal to 100 minus the coverage. Discharge is always reported as estimated total discharge even when coverage is less than 100 %.

4.4.3 Invalid thickness

The thickness data appear to be incorrect in some locations. For example, many locations have fast-flowing ice but report ice thickness as 10 m or less (Fig. a). We accept all ice thickness greater than 20 m and construct from this a thickness vs. ${\log }_{\mathrm{10}}$-speed relationship. For all ice thickness less than or equal to 20 m thick we adjust thickness based on this relationship (Fig. b). We selected the 20 m thickness cutoff after visually inspecting the velocity distribution (Fig. a). This thickness adjustment adds 20 Gt yr${}^{-\mathrm{1}}$ to our baseline-period discharge estimate with no adjustment. In the Appendix  and Table  we discuss the discharge contribution of these adjusted pixels and a comparison among this and other thickness adjustments.

4.5 Discharge

We calculate discharge per pixel using density (917 kg m${}^{-\mathrm{3}}$), filtered and filled ice speed, projection-corrected pixel width, and adjusted ice thickness derived from time-varying surface elevation and a fixed bed elevation (Eq. ). We assume that any change in surface elevation corresponds to a change in ice thickness and thereby neglect basal uplift, erosion, and melt, which combined are orders of magnitude less than surface melting (e.g., ). We also assume depth-averaged ice velocity is equal to the surface velocity.

We calculate discharge using the gate-orthogonal velocity at each pixel and at each timestamp – all velocity estimates are gate orthogonal at all times, regardless of gate position, orientation, or changing glacier velocity direction over time.

Annual averages are calculated by taking the sparse daily data with gaps, linearly interpolating at each pixel between observations to daily resolution, and then averaging to annual. The difference between this method and averaging only the observed samples is $\sim \mathrm{3}$ % median (5 % average, and a maximum of 10 % when examining the entire ice sheet and all years in our data). It is occasionally larger at individual glaciers when a year has few widely spaced samples of highly variable velocity.

Discharge uncertainty

A longer discussion related to our and others' treatments of errors and uncertainty is in the Appendix, but here we describe how we estimate the uncertainty related to the ice discharge following a simplistic approach. This yields an uncertainty of the total ice discharge of approximately 10 % throughout the time series.

At each pixel we estimate the maximum discharge, $D_{max}$, from

4

Bottom panel: time series of ice discharge from the Greenland Ice Sheet. Dots represent when observations occurred (limited to coverage $>\mathrm{50}$ %). Orange stepped line is annual average (limited to three or more observations in a year). Coverage (percentage of total discharge observed at any given time) is shown in the top panel and also by the opacity of the dots' interior and error bars on lower panel. When coverage is $<\mathrm{100}$ %, total discharge is estimated and shown.

[Figure omitted. See PDF]

At the region scale, the SE glaciers (see Fig.  for regions) are responsible for 136 to 164 ($\pm \mathrm{12}$ %) Gt yr${}^{-\mathrm{1}}$ of discharge (approximately one-third of ice-sheet-wide discharge) over the 1986 through 2019 period. By comparison, the predominantly land-terminating NO, NE, and NW together were also responsible for about one-third of total ice sheet discharge during this time (Fig. ). The discharge from most regions has been approximately steady or declining for the past decade. The NW is the only region exhibiting a persistent long-term increase in discharge – from $\sim \mathrm{90}$ to $\mathrm{115}$ Gt yr${}^{-\mathrm{1}}$ (22 % increase) over the 1999 through 2017 period ($+\sim \mathrm{1.4}$ Gt yr${}^{-\mathrm{1}}$ or $+\sim \mathrm{1.2}$ % yr${}^{-\mathrm{1}}$). This 1999 through 2017 annual average increase in NW discharge offsets declining discharge from other regions, but the NW increase stopped in 2018 and discharge in the NW dropped by 5 Gt yr${}^{-\mathrm{1}}$ (4 %) in 2019. This NW decline is then offset by a SE region increase. The largest contributing region, SE, contributed a high of $\mathrm{164}\pm \mathrm{19}$ Gt in 2004 but dropped to $\sim \mathrm{150}\pm \mathrm{18}$ Gt yr${}^{-\mathrm{1}}$ for the past decade.

Figure 5

Bottom panel: time series of ice discharge by region. Same graphical properties as Fig. .

[Figure omitted. See PDF]

Focusing on eight major contributors at the individual sector or glacier scale (Fig. ), Sermeq Kujalleq (Jakobshavn Isbræ) has slowed down from an annual average high of $\sim \mathrm{51}$ Gt yr${}^{-\mathrm{1}}$ in 2013 to $\sim \mathrm{34}$ Gt yr${}^{-\mathrm{1}}$ in 2018, likely due to ocean cooling . We exclude Ikertivaq from the top eight because that gate spans multiple sectors and outlets, while the other top dischargers are each a single outlet. The 2013 to 2016 slowdown of Sermeq Kujalleq (Fig. ) is compensated for by the many glaciers that make up the NW region (Fig. ). The large 2017 and 2018 reduction in discharge at Sermeq Kujalleq is partially offset by a large increase in the second-largest contributor, Helheim Gletsjer (Helheim Glacier; Fig. ), and a small increase in the third-largest contributor, Kangerlussuaq . Helheim discharged more ice than Sermeq Kujalleq in early 2018 and for all data estimates to date (through March) in 2020, although error bars still overlap.

Figure 6

Bottom panel: time series of ice discharge showing the eight major discharging glaciers from Fig. . Same graphical properties as Fig. .

[Figure omitted. See PDF]

6 Discussion

Different ice discharge estimates among studies likely stem from three categories: (1) changes in true discharge, (2) different input data (ice thickness and velocity), and (3) different assumptions and methods used to analyze data. Improved estimates of true discharge are the goal of this and many other studies, but changes in true discharge (category 1) can happen only when a work extends a time series into the future because historical discharge is fixed. Thus, any interstudy discrepancies in historical discharge must be due to category 2 (different data) or category 3 (different methods). Most studies use both updated data and new or different methods but do not always provide sufficient information to disentangle the two. This is inefficient. To more quantitatively discuss interstudy discrepancies, it is imperative to explicitly consider all three potential causes of discrepancy. Only when results are fully reproducible – meaning all necessary data and code are available (see ) – can new works confidently attribute discrepancies relative to old works. Therefore, in addition to providing new discharge estimates, we attempt to examine discrepancies among our estimates and other recent estimates. Without access to code and data from previous studies, it is challenging to take this examination beyond a qualitative discussion.

The algorithm-generated gates we present offer some advantages over traditional handpicked gates. Our gates are shared publicly, are generated by code that can be audited by others, and are easily adjustable within the algorithmic parameter space. This both allows sensitivity testing of gate location (Fig. ) and allows gate positions to systematically evolve with glacier termini (not done here). The total ice discharge we estimate is $\sim \mathrm{10}$ % less than the total discharge of two previous estimates and similar to that of , who attributes their discrepancy with to the latter using only summer velocities, which have higher annual average values than seasonally comprehensive velocity products. The gate locations also differ among studies, and glaciers with baseline velocity less than 100 m yr${}^{-\mathrm{1}}$ are not included in our study due to our velocity cutoff threshold, but this should not lead to substantially different discharge estimates (Fig. ).

Our gate selection algorithm also does not place gates in northeast Greenland at Storstrømmen, Bredebræ (Bredebrae), or their confluence, because during the baseline period that surge glacier was in a slow phase. We do not manually add gates at these glaciers. The last surge ended in 1984 , prior to the beginning of our time series, and these glaciers are therefore not likely to contribute substantial discharge even in the early period of discharge estimates.

We instead attribute the majority of our discrepancy with to the use of differing bed topography in southeast Greenland. When we compare our top 10 highest discharging glaciers in 2000 with those reported by , we find that the Køge Bugt (also known as Køge Bay) discharge reported by is $\sim \mathrm{31}$ Gt, but our estimate is only $\sim \mathrm{16}$ Gt ($\sim \mathrm{17}$ Gt in , and similar in ). The bed elevation dataset that likely uses the same bed data employed by has a major depression in the central Køge Bugt bed. This region of enhanced ice thicknesses is not present in the BedMachine dataset that we, , and employ (Fig. ). If the Køge Bugt gates of are in this location, then those gates overlie ice thicknesses that are about twice those reported in BedMachine v3. With all other values held constant, this results in roughly twice the discharge. Although we do not know whether BedMachine or is more correct, conservation of mass suggests that a substantial subglacial depression should be evident as either depressed surface elevation or velocity .

We are unable to attribute the remaining discrepancy between our discharge estimates and those by . It is likely a combination of different seasonal velocity sampling , our evolving surface elevation from , or other previously unpublished algorithmic or data differences, of which many possibilities exist.

Our ice discharge estimates agree well with the most recently published discharge estimate (, also used by ), except that our discharge is slightly less. We note that our uncertainty estimates include the estimates, but the opposite does not appear be true. The minor differences are likely due to different methods. use seasonally varying ice thicknesses, derived from seasonally varying surface elevations, and a Monte Carlo method to temporally interpolate missing velocity data to produce discharge estimates. In comparison, we use linear interpolation of both yearly surface elevation estimates and temporal data gaps. It is not clear whether linear or higher-order statistical approaches are best suited for interpolation as annual cycles begin to shift, as is the case with Sermeq Kujalleq (Jakobshavn Isbræ) after 2015. There are benefits and deficiencies with both methods. Linear interpolation may alias large changes if there are no other observations nearby in time. Statistical models of past glacier behavior may not be appropriate when glacier behavior changes.

It is unlikely that discharge estimates using gates that are only approximately flow orthogonal and time invariant have large errors due to this, because it is unlikely that glacier flow direction changes significantly, but our gate-orthogonal treatment may be the cause of some differences among our approach and other works. Discharge calculated using nonorthogonal methodology would overestimate true discharge.

7 Code and data availability

This work in its entirety is available at https://doi.org/10.22008/promice/data/ice_discharge . The glacier-scale, sector, region, and Greenland summed ice sheet discharge dataset is available at https://doi.org/10.22008/promice/data/ice_discharge/d/v02 , where it will be updated as more velocity data become available. The gates can be found at https://doi.org/10.22008/promice/data/ice_discharge/gates/v02 , the code at https://doi.org/10.22008/promice/data/ice_discharge/code/v0.0.1 , and the surface elevation change at https://doi.org/10.22008/promice/data/DTU/surface_elevation_change/v1.0.0 .

8 Conclusions

We have presented a novel dataset of flux gates and a 1986 through 2019 glacier-scale ice discharge estimate for the Greenland Ice Sheet. These data are underpinned by an algorithm that both selects gates for ice flux and then computes ice discharges.

Our results are similar to the most recent discharge estimate but begin in 1986 – although there are fewer samples prior to 2000. From our discharge estimate we show that, over the past $\sim \mathrm{30}$ years, ice sheet discharge was $\sim \mathrm{440}$ Gt yr${}^{-\mathrm{1}}$ prior to 2000, rose to over 500 Gt yr${}^{-\mathrm{1}}$ from 2000 to 2005, and has held roughly steady since 2005 at near 500 Gt yr${}^{-\mathrm{1}}$. However, when viewed at a region or sector scale, the system appears more dynamic with spatial and temporal increases and decreases canceling each other out to produce the more stable ice sheet discharge. We note that there does not appear to be any dynamic connection among the regions, and any increase in one region that was offset by a decrease in another has likely been due to chance. If in coming years when changes occur the signals have matching signs, then ice sheet discharge would decrease or increase, rather than remain fairly steady.

The application of our flux-gate algorithm shows that ice-sheet-wide discharge varies by $\sim \mathrm{30}$ Gt yr${}^{-\mathrm{1}}$ due only to gate position, or $\sim \mathrm{40}$ Gt yr${}^{-\mathrm{1}}$ due to gate position and cutoff velocity (Fig. ). This variance is approximately equal to the uncertainty associated with ice-sheet-wide discharge estimates reported in many studies (e.g., ). We highlight a major discrepancy with the ice discharge data of , and we suspect this discharge discrepancy – most pronounced in southeast Greenland – is associated with the choice of digital bed elevation model, specifically a deep hole in the bed at Køge Bugt.

Transparency in data and methodology are critical to move beyond a focus of estimating discharge quantities towards more operational mass loss products with realistic errors and uncertainty estimates. The convention of devoting a paragraph, or even page, to methods is insufficient given the complexity, pace, and importance of Greenland Ice Sheet research . Therefore the flux gates, discharge data, and the algorithm used to generate the gates, discharge, and all figures from this paper are available. We hope that the flux gates, data, and code we provide here is a step toward helping others both improve their work and discover the errors in ours.

Appendix A Errors and uncertainties

Here we describe our error and uncertainty treatments. We begin with a brief philosophical discussion of common uncertainty treatments, our general approach, and then the influence of various decisions made throughout our analysis, such as gate location and treatments of unknown thicknesses.

Traditional and mathematically valid uncertainty treatments divide errors into two classes: systematic (bias) and random. The primary distinction is that systematic errors do not decrease with more samples, and random errors decrease as the number of samples or measurements increases. The question is then which errors are systematic and which are random. A common treatment is to decide that errors within a region are systematic and among regions are random. This approach has no physical basis – two glaciers a few hundred meters apart but in different regions are assumed to have random errors, but two glaciers thousands of kilometers apart but within the same region are assumed to have systematic errors. It is more likely the case that all glaciers narrower than some width or deeper than some depth have systematic errors even if they are on opposite sides of the ice sheet, if ice thickness is estimated with the same method (i.e., the systematic error is likely caused by the sensor and airplane, not the location of the glacier).

The decision to have $R$ random samples (where $R$ is the number of regions, usually $\sim \mathrm{18}$ based on ) is also arbitrary. Mathematical treatment of random errors means that, even if the error is 50 %, 18 measurements reduce it to only 11.79 %.

This reduction is unlikely to be physically meaningful. Our 173 sectors, 267 gates, and 5829 pixels means that, even if errors were 100 % for each, we could reduce it to 7.5 %, 6.1 %, or 1.3 % respectively. We note that the area error introduced by the common EPSG:3413 map projection is $-\mathrm{5}$ % in the north and $+\mathrm{8}$ % in the south. While this error is mentioned in some other works (e.g., ) it is often not explicitly mentioned.

We do not have a solution for the issues brought up here, except to discuss them explicitly and openly so that those, and our own, error treatments are clearly presented and understood to likely contain errors themselves.

A1 Invalid thickness

We assume ice thicknesses $<\mathrm{20}$ m are incorrect where ice speed is $>\mathrm{100}$ m yr${}^{-\mathrm{1}}$. Of 5829 pixels, 5194 have valid thickness, and 635 (12 %) have invalid thickness. However, the speed at the locations of the invalid thicknesses is generally much less (and therefore the assumed thickness is less), and the influence on discharge is less than an average pixel with valid thickness (Table ).

Table A1

Statistics of pixels with and without valid thickness. Numbers represent speed (m yr${}^{-\mathrm{1}}$) except for the “count” row.

	Good pixels	Bad pixels
Count	5205	624
Mean	857	272
SD	1117	239
Minimum	100	100
25 %	236	130
50 %	506	181
75 %	995	291
Maximum	10 044	1505

When aggregating by gate, there are 267 gates. Of these, 179 (67 %) have no bad pixels and 88 (33 %) have some bad pixels, 64 have $>\mathrm{50}$ % bad pixels, and 62 (23 %) are all bad pixels.

We adjust these thickness using a poor fit (correlation coefficient: 0.3) of the ${\log }_{\mathrm{10}}$ of the ice speed to a thickness where the relationship is known (thickness $>\mathrm{20}$ m). We set errors equal to one half the thickness (i.e., ${\mathit{\sigma }}_H=\pm \mathrm{0.5}H$). We also test the sensitivity of this treatment to simpler treatments and have the following five categories: NoAdj.

No adjustments made. Assume BedMachine thicknesses are all correct.

NoAdj+Millan.

Same as NoAdj but using bed where available.

300.

If a gate has some valid pixel thicknesses, set the invalid thicknesses to the minimum of the valid thicknesses. If a gate has no valid thickness, set the thickness to 300 m.

400.

Set all thicknesses $<\mathrm{50}$ to 400 m

Fit.

Use the thickness–speed relationship described above.

Table  shows the estimated baseline discharge to these four treatments.

Table A2

Effect of different thickness adjustments on baseline discharge.

Treatment	Discharge (Gt)
NoAdj	$\mathrm{472}\pm \mathrm{49}$
NoAdj+Millan	$\mathrm{481}\pm \mathrm{49}$
300	$\mathrm{489}\pm \mathrm{49}$
400	$\mathrm{495}\pm \mathrm{52}$
Fit	$\mathrm{493}\pm \mathrm{51}$

Finally, Fig.  shows the geospatial locations, concentration, and speed of gates with and without bad pixels.

Figure A1

Gate locations and thickness quality. (a) Locations of all gates. Black dots represent gates with 100 % valid thickness pixels, blue with partial, and red with none. (b) Percent of bad pixels in each of the 267 gates, arranged by region. (c) Average speed of gates. Color same as left panel.

[Figure omitted. See PDF]

A2 Missing velocity

We estimate discharge at all pixel locations for any time when there exists any velocity product. Not every velocity product provides velocity estimates at all locations, and we fill in where there are gaps by linearly interpolating velocity at each pixel in time. We calculate coverage, the discharge-weighted percent of observed velocity at any given time (Fig. ), and display coverage as (1) line plots over the time series graphs, (2) opacity of the error bars, and (3) opacity of the infilling of time series dots. Linear interpolation and discharge-weighted coverage is illustrated in Fig. , where pixel A has a velocity value at all three times, but pixel B has a filled gap at time $t_{\mathrm{3}}$. The concentration of valid pixels is 0.5, but the weighted concentration, or coverage, is $\mathrm{9}/\mathrm{11}$ or $\sim \mathrm{0.82}$. When displaying these three discharge values, $t_{\mathrm{1}}$ and $t_{\mathrm{4}}$ would have opacity of 1 (black), and $t_{\mathrm{3}}$ would have opacity of 0.82 (dark gray).

Appendix B

Differences between BedMachine and near Køge Bugt. Panel (a) is baseline ice speed, (b) BedMachine thickness, (c)  thickness, and (d) difference computed as BedMachine $-$ Bamber. The curved line is the gate used in this work.

[Figure omitted. See PDF]

Appendix C Sentinel-1 ice velocity maps

We use ESA Sentinel-1 synthetic aperture radar (SAR) data to derive ice velocity maps covering the Greenland Ice Sheet margin using offset tracking assuming surface parallel flow using the digital elevation model from the Greenland Ice Mapping Project (GIMP DEM, NSIDC 0645) by . The operational interferometric postprocessing (IPP) chain , developed at the Technical University of Denmark (DTU) Space and upgraded with offset tracking for ESA’s Climate Change Initiative (CCI) Greenland project, was employed to derive the surface movement. The Sentinel-1 satellites have a repeat cycle of 12 d, and due to their constellation, each track has a 12 d repeat cycle. We produce a Greenland-wide product that spans two repeat cycles of Sentinel-1A. The product is a mosaic of all the ice velocity maps based on 12 d pairs produced from all the tracks from Sentinel-1A and 1B covering Greenland during those two cycles. The product thus has a total time span of 24 d. Twelve-day pairs are also included in each mosaic from track 90, 112, and 142 covering the ice sheet margin in the south as well as other tracks on an irregular basis in order to increase the spatial resolution. and have exploited the high temporal resolution of the product to investigate dynamics of glaciers. The maps are available from 13 September 2016 and onward, are updated regularly, and are available from http://promice.org (last access: 6 June 2020).

Appendix D Software

This work was performed using only open-source software, primarily GRASS GIS and Python , in particular the Jupyter , pandas , numpy , statsmodel , x-array , and Matplotlib packages. The entire work was performed in Emacs using Org mode . The parallel tool was used to speed up processing. We used proj4 to compute the errors in the EPSG 3413 projection. All code used in this work is available at https://github.com/mankoff/ice_discharge (last access: 6 June 2020).

Author contributions

KDM conceived of the algorithm approach and wrote the code. KDM, WIC, and RSF iterated over the algorithm results and methods. ASO provided the velocity data. SAK supplied the surface elevation change data. All authors contributed to the scientific discussion, writing, and editing of the manuscript.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

We thank the contributors and coauthors of the previous edition , as well as the reviewers and editors for their constructive input that helped improve the paper. Sentinel ice velocity maps were produced from Copernicus Sentinel-1 image data, processed by ESA data as part of PROMICE, and were provided by the Geological Survey of Denmark and Greenland (GEUS) at http://www.promice.org (last access: 6 June 2020).

Financial support

This research has been supported by the Programme for Monitoring of the Greenland Ice Sheet (PROMICE) and the European Union's Horizon 2020 research and innovation program (INTAROS, grant no. 727890).

Review statement

This paper was edited by Reinhard Drews and reviewed by two anonymous referees.