Greenland liquid water discharge from 1958 through 2019

Greenland runoff, from ice mass loss and increasing rainfall, is increasing. That runoff then discharges and impacts the physical, chemical, and biological properties of the adjacent fjords. However, where and when the discharge occurs is not readily available in an open database. Here we provide datasets of high-resolution Greenland hydrologic outlets, basins, and streams, and a daily 1958 through 2019 time series of Greenland liquid water discharge for each outlet. The data include 24507 ice marginal outlets and upstream basins, and 29635 land coast outlets and upstream basins, derived from the 100 m 5 ArcticDEM and 150 m BedMachine. At each outlet there are daily discharge data for 22644 days ice sheet runoff routed subglacially to ice margin outlets, and land runoff routed to coast outlets from two regional climate models (RCMs; MAR and RACMO). Our sensitivity study of how outlet location changes for every inland cell based on subglacial routing assumptions, shows that most inland cells where runoff occurs are not highly sensitive to those routing assumptions, and outflow location does not move far. We compare RCM results with 10 streams instrumented with gauges spanning four orders of magnitude 10 of daily discharge. Results show that for daily discharge at individual basin scale the 95 % prediction interval generally falls within plus-or-minus a factor of five (half an order of magnitude, or +500%/-80%). Results from this study are available at doi:10.22008/promice/freshwater (Mankoff, 2020a) and code is available at http://github.com/mankoff/freshwater (Mankoff, 2020b).


15
Over the past decades, liquid runoff from Greenland has increased (Mernild and Liston, 2012;Bamber et al., 2018;Trusel et al., 2018;Perner et al., 2019) contributing to mass decrease (Sasgen et al., 2020). When that runoff leaves the ice sheet and discharges into fjords and coastal seas, it influences a wide range of physical (Straneo et al., 2011;Mortensen et al., 2013;Bendtsen et al., 2015;Cowton et al., 2015;Mankoff et al., 2016;Fried et al., 2019;Cowton et al., 2019;Beckmann et al., 2019), chemical (Kanna et al., 2018;Balmonte et al., 2019), and biological (Kamenos et al., 2012;Kanna et al., 2018;20 Balmonte et al., 2019) systems (Catania et al., 2019). The scales of the impacts range from instantaneous at the ice-ocean boundary to decadal in the distal ocean (Gillard et al., 2016). The influence of freshwater on multiple domains and disciplines (Catania et al., 2019) is the reason several past studies have estimated runoff and discharge at various temporal and spatial In Eq. 1, M E is melt, RA is rainfall, RT is retention, and RF is refreezing. In RACMO, retention occurs only when firn is present (not with bare ice). MAR does have a delay for bare ice runoff. Neither have a delay for land runoff. Both RCM outputs 55 were provided regridded to the same 1 km grid using an offline statistical down-scaling technique based on local vertical runoff gradient applied to the sub-grid topography (Noël et al., 2016;Fettweis et al., 2020). MAR (v 3.11;Delhasse et al. (2019)) ran with 7.5 km resolution and ERA5 6-hour forcing. RACMO (v 2.3p2;Noël et al. (2018Noël et al. ( , 2019) ran with 5.5 km resolution and ERA-Interim 6-hour forcing. Runoff is assigned an uncertainty of ±15 % (Sect. 4.3.3). 60 We use 10 river discharge daily time series to validate the results of this work. The name, location, time coverage, and relevant data and scientific publications associated with each of these observational data are listed in Table 1. water). When referring to the surface classification in the RCM, we explicitly state "RCM mask".

River discharge observations
-MAR and RACMO refer to the RCMs, but when comparing discharge estimates between them or to observations, we use MAR and RACMO to refer to our discharge product derived from the MAR and RACMO RCM runoff variables, rather than repeatedly explicitly stating "discharged derived from [MAR|RACMO] runoff". The use should be clear from 75 context.

Streams, outlets, and basins
Streams are calculated from the hydraulic head h which is the DEM surface for land surface routing, or the subglacial pressure head elevation for subglacial routing. h is defined as We introduce the graphics here as part of the methods to reduce replication in figure captions -we show 10 nearly identical graphics (Figs. 7 and 9 through 17) for 10 different observation locations, and each graphic uses the same template of six panels.
For each figure (Figs. 7,9 to 17), the top panel (a) shows a satellite basemap with the land portion of the basin of interest (if 135 it exists) outlined in dark green, the streams within that basin in light green, the basin outlet as an orange filled diamond, and the stream gauge location as an orange unfilled diamond. Ice basin(s) that drain to the outlet are outlined in thick dark blue if they exist, and all other ice basins in thin dark blue. Both MAR and RACMO use the same domains. The RCM ice domain is in light blue, and the RCM land domain is not shown, but is outside the light blue ice domain (not including the water). The scale of each map varies, but the basins lines (green and dark blue) are on a 100 m grid, and the RCM grid cells (light blue) are of RCM land are discarded because the basins classify part or all of these cells as ice or ocean. Of the 1,781,816 km 2 of basin 165 ice cells, the RCMs cover 1,760,912 km 2 (~99 %) with their ice cells, and 20,904 km 2 (~1 %) of basin grid cells are filled in (the RCMs have these as land or ocean). 21,793 km 2 of RCM ice are discarded, because the basins classify part or all of these cells as land or ice (Table and

Comparison with previous similar work
Our static products -streams, outlets, and basins -have been previously estimated. Lewis and Smith (2009) identified 293 distinct hydrologic ice basins and provided a data set of ice basins and ice margin outlets. Our work, a decade later, has~2 175 orders of magnitude more basins and outlets because of the higher resolution of the input data, and includes additional data. We provide ice basins, ice margin outlets, ice streams with metadata, land basins, coastal outlets, and land streams with metadata.
Lewis and Smith (2009) generated basins from a 5 km DEM, compared to the 100 m DEM used here. Routing with a 5 km DEM that does not capture small-scale topography is likely to cause some basins and outlets to drain into an incorrect fjordwe find that some land basins delineated with even the 150 m BedMachine land surface may drain into the incorrect fjord, but 180 we did not find similar errors with the 100 m ArcticDEM product used in this work.
Our time-series product -discharge, also has existing similar products. The most recent of these is from Bamber et al. (2018) who provide a data product at lower spatial resolution (5 km), lower temporal resolution (monthly), and only coastal discharge, not coastal basins, ice basins, or ice margin outlets and discharge. However, Bamber et al. (2018) surpasses our product in that spatial coverage includes a larger portion of the Arctic including Iceland, Svalbard, and Arctic Canada. Furthermore, by a factor of 0.78 (MAR) and 0.73 (RACMO), and 5 to 95 % quantile of 0.30 to 2.08. The top 2/3rds of observed discharge spans 230~2 orders of magnitude (width of horizontal lines, from~10 1 to~10 3 m 3 s -1 ). The ratio of the RCMs to these observations for the top 2/3rds has a 5 to 95 % quantile range of~1 order of magnitude (distance between horizontal lines, from log 10 0.3 to log 10 2.08 = 0.84). The 5 to 95 % quantile range of the ratio between the RCMs and the observations is therefore half the range of the observations. Put differently, days with high observed discharge may have modeled discharge within ±0.5 order of magnitude, or plus-or-minus a factor of five, or +500/-80 %. The modeled discharge is not likely to move farther than this 235 from the observations, and high discharge remains high.
The bottom third of discharge is where the largest disagreement occurs. The mean model values are near the observed -the ratio of RCM to observed discharge is scaled by 0.67 for MAR (~33 % low) and 1.08 for RACMO (~8 % high), but the 5 to 95 % quantile range of the ratio between RCM and observations is large. Although large uncertainties for low discharge may not seem to matter for some uses (e.g. estimates of total discharge from Greenland, which is dominated by the largest quantities 240 of discharge), it may matter for other uses. The bottom 1/3 quantile of observed discharge spans 3 orders of magnitude (10 -2 to~10 1 ) but the uncertainty of the RCM to observations ratio spans~4 and~2 orders of magnitude for MAR and RACMO respectively (~10 -3 to~2.2x10 1 MAR;~10 -1 to 2.2x10 1 RACMO).

245
The partial (last calendar year) discharge time series shows MAR and RACMO agree well with each other, but have a maximum of 500 m 3 s -1 while observations are up to 4x more (Fig. 7b). Low discharge (both early and late season) is over-estimated and high discharge is under-estimated, approximately equal for both RCMs (Fig. 7c). The low discharge over-estimate ranges from a mean multiple of 1.68 (MAR) and 1.57 (RACMO) to a +95 % quantile ratio of~70 (MAR) and~52 (RACMO). The highdischarge under-estimate has a mean multiple of~0.5 for both MAR and RACMO, and a 5 to 95 quantile range of between 250 0.23 to 1.09.
The Watson River discharge presented here is approximately half of the van As et al. (2018) discharge for high discharge.
The large underestimate for high discharge may be due to either errors in the basin delineation used in this study, errors in the stage-discharge relationship used by , errors in the RCM runoff estimates, or a combination of the above three. All three of these error sources increase with high discharge (and associated melt): Basin delineation becomes 255 less certain with inland distance from the ice sheet margin. The river stage-discharge conversion becomes less certain at high stage levels. Runoff calculations become less certain from a snow surface than an ice surface, because of e.g. snow density, subsurface refreezing, and surface darkening.
The complexity of estimating the area of the Watson River catchment is described by Monteban et al. (2020), who note that previous studies have used values ranging from 6131 km 2 (Mernild et al., 2010b) to 12547 km 2 (van As et al., 2012). Our basin 260 is smaller than the basin used in van  and similar to Mernild et al. (2018) who attributed the difference between their modeled outflow and observations from van  to their decision to use surface rather than subglacial routing, and applied a correction term. We find that our basin does not include a separate basin to the south that is part of the Watson River ice basin in van As et al. (2018) (from Lindbäck et al. (2015) and Lindbäck et al. (2014)). We are able to recreate the van To solve this, we manually select two large ice basins to the south of the Watson River ice basin. Modeled and observed discharge agree after including these two basins (Fig. 8), suggesting basin delineation, not stage-discharge or RCM runoff is the primary cause for this disagreement. Furthermore, it is the additional width at lower elevation from the larger basin, not the increased inland high-elevation area, that likely contributes the runoff needed to match At the Watson River outlet, there is no reason to suspect this product underestimates observed discharge by 50 %. The observations are needed to highlight the disagreement. Once this disagreement is apparent, it is also not clear what to do to reduce the disagreement without the previous efforts by Lindbäck et al. (2015) and Lindbäck et al. (2014). Basin delineation is discussed in more detail in the uncertainty section (Sect. 4.3.2). The other two "problematic" areas highlighted above (Sect. 4.2) can be detected and improved without observational support.

Leverett Glacier
The Leverett Glacier basin area is 1361 km 2 and 100 % ice (Fig 9a). The partial (last calendar year) discharge time series shows MAR and RACMO agree well with each other and with the observations (Fig. 9b), with no seasonal dependence (Fig 9c). The 5 to 95 % prediction interval for MAR is generally within the 1:5 and 5:1 bands, with a larger spread for RACMO (Fig 9c).

280
High model discharge is 3 % higher than observed (MAR) or 25 % higher than observed (RACMO), and the 5 to 95 quantile range of the ratio is between 0.73 and 1.62 (MAR) and 0.83 and 2.02 (RACMO). Low model discharge is also centered near the observations, but as always larger errors exist for low discharge (Fig 9d).
This basin is problematic because the basin feeding the outlet is small (< 5 km 2 ), but even without the observational record satellite imagery shows a large river discharging from the ice sheet here. Meanwhile, a large (100s of km 2 ) ice basin does 285 discharge just a few 100 m away, but not upstream of this gauge location. We therefore adjust the gauge location onto the ice (equivalent to selecting a different outlet) so that our database access software selects what appears to be the correct basin given the size of the stream in the satellite imagery ( Fig. 9).
The plots shown here use the adjusted gauge location and modeled discharge appears to match the observed discharge. When plotting (not shown) the modeled discharge for the outlet just upstream of the true gauge location, results are clearly incorrect.

290
This issue -small basins at the margin and incorrect outlet location -is persistent throughout this product and discussed in more detail in Sect. 4.3.2.
The Leverett Glacier basin is a subset of the Watson River outlet basin (Sect. 4.2.2). The strong agreement here supports our claim that the Watson River disagreement is not from the RCM runoff or the stage-discharge relationship, but more likely due to basin area. The correct Watson River basin should include some basins outside of the Leverett Glacier basin that still drain

Kiattuut Sermiat
The Kiattuut Sermiat discharge basin area is 693 km 2 , of which 391 km 2 (56 %) are land and 302 km 2 (44 %) are ice. The basin area is incorrectly large because the land basin reported and shown includes the entire basin that contains the discharge point, of which some is downstream (Fig 10a). However, only~25 % of runoff comes from the land, and only a small portion of the 300 land basin is downstream of the gauge location, so this is not enough to explain the discharge vs. observation disagreement.
The partial (last calendar year) discharge time series shows MAR and RACMO agree well with each other, but are significantly higher than the observations (Fig. 10b). Both low and high discharge are over-estimated, but the 5 to 95 % quantile range of the ratio are within a factor of five (Fig 10c), with a mean ratio between 1.71 (RACMO bottom 1/3rd of discharge) to 2.44 (MAR high 2/3rds discharge)

305
The Kiattuut Sermiat gauge is in a problematic location in terms of determining the actual (non-theoretical) upstream contributing area. Similar to the Leverett Glacier gauge location, the issues here can be estimated independent of observational data. Specifically, it is not clear if this stream includes water from the larger glacier to the east and east-northeast that feeds this glacier ( Fig. 10a) -in our delineation it does not. Furthermore, several glaciers to the north-northwest and detached from the glacier near the stream gauge appear to drain into a lake that then drains under the glacier and then to the stream gauge. This 310 latter issue is observable in any cloud-free satellite imagery and does not need the basin delineations provided here to highlight the complexities of this field site. Nonetheless, RCM discharge estimates are only slightly more than double the observations.
The Kiattuut Sermiat gauge location may have been selected in part due to its accessibility -it is walking distance from the Narsarsuaq airport. The data may also suit their intended purpose well and there are likely many results that can be derived independent of the area or location of the upstream source water. However, if the location or area of the upstream contributions 315 are important, then gauge location should balance ease of access and maintenance with the ease with which the data can be interpreted in the broader environment.

GEM observations near Nuuk
Four Greenland Ecosystem Monitoring Programme (GEM) stream gauges are located near Nuuk with similar basin properties.
All are small (7.56 to 37.52 km 2 ), and 10 % to 25 % ice in the basin mask, but two of the four (Kingigtorssuaq (Fig. 11) and 320 Oriartorfik ( Fig. 12)) contain small glaciers contributing to observed discharge but no RCM ice cells cover those glaciers, and the remaining two (Teqinngalip (Fig. 13) and Kobbefjord (Fig. 14)) have several small glaciers, but only one per basin has RCM ice coverage.
All four of these basins show some weak agreement. The maximum r 2 is 0.47 (Fig. 13c) and the minimum is 0.11 (Fig   11c), but we note that the worst agreement comes from a basin with no glaciers in the RCM domain, and that in all cases the
The Røde Elv basin is situated at the southern edge of Disko Island (Fig. 6). It has an area of 100 km 2 , of which 72 km 2 330 are land and 28 km 2 are ice (Fig 15a). The partial (last calendar year) discharge time series shows MAR, RACMO, and the observations all in approximately the same range but with high variability (Fig. 15b). Of the few samples here (n = 98), most are within the factor-of-five bands for MAR and a few more are outside the bands for RACMO (Fig. 15c). Mean discharge offset ranges from a ratio of 0.82 (RACMO low) to 1.85 (MAR low), with high discharge estimates slightly closer to observationsa 48 % and 77 % overestimate for MAR and RACMO respectively (Fig. 15d).

335
The Zackenberg basin in NE Greenland has an area of 487 km 2 , of which 378 km 2 (78 %) are land and 109 km 2 (22 %) are ice (Fig. 16a). The partial (last calendar year) discharge time series shows disagreement between MAR and RACMO that generally bound the observations (Fig. 16b). RACMO-derived discharge is consistently high for low discharge early in the year, but both discharge products fall mostly within the factor-of-five bands (Fig. 16c). For high discharge, mean modeled discharge is 9 % high (MAR) and 24 % low (RACMO), and has worst-case 5 to 95 % quantile range low by a factor of 0.29 (Fig. 16d).

340
The Qaanaaq basin in NW Greenland has an area of 13.2 km 2 , of which 2.2 km 2 (17 %) are land and 11 km 2 (83 %) are ice (Fig. 17a). The partial (last calendar year) discharge time series shows disagreement between MAR and RACMO that generally bound the observations (Fig 17b). Of the few samples (n = 82), MAR preferentially over-estimates and RACMO under-estimates discharge, but both generally within a factor of 5 (Fig 17c). The mean high discharge ratio is 1.26 (MAR) and 0.4 (RACMO) from Fig. 17d.

Uncertainty
The volume of data generated here is such that manually examining all of it or editing it to remove artifacts or improve the data would be time and cost prohibitive. A similar warning is provided with the ArcticDEM data used here. However, any ArcticDEM issues interior to a basin do not impact results here that are aggregated by basin and reported at the outlet.
ArcticDEM issues that cross basin boundaries should only impact a small part of the basin near the issue.

350
Uncertainty from RCM inputs and observations are considered external to this work, although they are still discussed (Sects.

and 4.3.4).
In this work, we introduce one new source of uncertainty -the routing model, which generates both temporal (runoff delay) and spatial (basin delineation) uncertainty.

Temporal uncertainty
The RCMs include a time lag between when water melts in the model and when it leaves a grid cell. RACMO retention occurs 355 only when there is firn cover (no retention when bare ice melts); MAR includes a time delay of up to 10 days that is primarily a function of surface slope (Zuo and Oerlemans, 1996;Yang et al., 2019). However, neither model includes a subglacial system.
Properly addressing time delays with runoff requires addressing storage and release of water across a variety of timescales in a variety of media: firn (e.g. Munneke et al. (2014); Vandecrux et al. (2019)), supraglacial streams and lakes (e.g. Zuo and terrestrial streams and lakes (e.g. van As et al. (2018)) and a variety of other physical processes that are not within the scope of surface mass balance (SMB) modeling. Runoff delay can be implemented outside the RCMs (e.g. Liston and Mernild (2012); Mernild et al. (2018)), but for this version of the product we assume that once an RCM classifies meltwater as "runoff", it is instantly transported to the outlet. Actual lags between melt and discharge range from hours to years (Colgan et al., 2011;van As et al., 2017;Rennermalm et al., 2013;Livingston et al., 2013).

365
Data released here includes no additional lag beyond the RCM lag, although a 7-day running mean (cf. van As et al. (2017)) is included in all of the results presented here except Fig. 6 which shows monthly summed data, and Fig. 4 which shows yearly summed data. When increasing the signal to noise by summing by year (Fig. 4 vs. Fig. 3), model results more closely match observations. scenario of k = 1.11 describing an over-pressured system. Because ρ w /ρ i = 1.09, setting k = 1.09 reduces Eq. 2 to h = z s , and is equivalent to an over-pressured system with surface routing of the water. In a limited examination comparing our results with k ∈ [0.8, 0.9, 1.0], we did not detect basins with large changes in upstream area. In addition, all time series graphics show the mean RCM discharge for k = 1.0, but the uncertainty among all three k values (not shown) is small enough that it is difficult to distinguish the three separate uncertainty bands -the difference between RCMs or between RCMs and observations is much 390 larger than uncertainty from the k parameter.

Basin uncertainty
The above issues are specific to ice basins. Land basin outlets do not change location, and the range of upstream runoff from different k simulations to a land outlet provides one metric of uncertainty introduced by the k parameter. This uncertainty among all three k values is small at ice margin outlets. It is even smaller at land outlets which act as spatial aggregators and increase the signal-to-noise ratio.

395
Below, we discuss the known uncertainties, ranging from least to most uncertain.
The basins presented here are static approximations based on 100 m DEM of a dynamic system. Land basin boundaries are likely to be more precise and accurate than ice basins, because the land surface is better resolved, has larger surface slopes, has negligible sub-surface flow, and is less dynamic than the ice surface. Even if basins and outlets seem visually correct from the 100 m product, the basin outline still has uncertainty on the order of hundreds of meters and will therefore include many minor 400 errors and non-physical properties, such as drainage basin boundaries bisecting lakes. However, all artefacts we did find are significantly smaller than the 1 km 2 grid of the RCM inputs. We do not show but note that when doing the same work with the 150 m BedMachine land surface DEM, some basins change their outlet locations significantly -draining on the opposite side of a spit or isthmus and into a different fjord than the streams do when observed in satellite imagery. We have not observed these errors in streams and basins derived from the 100 m ArcticDEM in a visual comparison with Google Earth, although they 405 may still exist.
Moving from land basins to subglacial ice basins, the uncertainty increases because subglacial routing is highly dynamic on timescales from minutes to seasons (e.g. Werder et al. (2013)). This dynamic system may introduce large spatial changes in outflow location (water or basin "piracy", Ahlstrøm et al. (2002) Near the ice margin quality issues exist. At the margin, many of the small basins (absorbed or not) may be incorrect because the bed uncertainty is larger relative to the ice thickness, and therefore uncertainty has a larger influence on routing. Minor 420 mask mis-alignments may cause hydraulic jumps (waterfalls) at the margin, or sinks that then need to be filled by the algorithm, and may overflow (i.e. the stream continues onward) somewhere at the sink edge different from the location of the real stream.
The solution for individual outlets is to visually examine modeled outlet location, nearby streams in satellite imagery, and the area of upstream catchments, as we did for the Leverett Glacier outlet (Sect 4.2.3). Alternatively, selecting several outlets in an area will likely include the nearby "correct" outlet. This can be automated and an effective method to aggregate all the 425 micro-ice basins that occur at the domain edge is to select the downstream land basin associated with one ice outlet, and then all upstream ice outlets for that land basin.

RCM uncertainty
In addition to the basin delineation issues discussed above, the runoff product from the RCMs also introduces uncertainty into the product generated here. The RCM input products do not provide formal time-or space-varying error estimates, but of course 430 do contain errors because they represent a simplified and discretised reality. RCM uncertainty is shown here with a value of ±15 %. The MAR uncertainty comes from an evaluation by the Greenland SMB Model Intercomparison Project (GrSMBMIP; Fettweis et al. (2020)) that examined the uncertainty of modelled SMB for 95 % of the 10767 in-situ measurements over the main ice sheet. 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 435 time scale. The RACMO uncertainty comes from an estimated average 5% runoff bias in RACMO2.3p2 compared to annual cumulative discharge from the Watson River (Noël et al., 2019). The bias increases to a maximum of 20 % for extreme runoff years (e.g. 2010 and 2012), so here we select 15 %, a value between the reported 5 % and the maximum 20 % that matches the MAR uncertainty. We display ±15 % uncertainty in the graphics here and suggest this is a minimum value for daily runoff data.

440
The 15 % RCM uncertainty is represented graphically in the time series plots when comparing to each of the observations.
It is not shown in the scatter plots because the log-log scaling and many points makes it difficult to display. In the time series plots, we show the mean value from the k = 1.0 scenario, and note that discharge from the other two k scenarios covered approximately the same range.

445
When comparing against observations, additional uncertainty is introduced because the stage-discharge relationship is neither completely precise or accurate. We use published observation uncertainty when it exists. Only two observational data sets come with uncertainty: Watson River and Qaanaaq. Similar to the RCM uncertainty, they are displayed in the time series but not in the scatter plots.

Mitigating uncertainties 450
Traditional uncertainty propagation is further complicated because it is not clear to what extent the three uncertainties (observational, RCM, and routing model) should be treated as independent from each other -all three uncertainties are likely to show some correlation with elevation, slope, air temperature, or other shared physical processes.
Many of the uncertainties discussed here can be mitigated by increasing the signal to noise ratio of the product. Because we provide a high spatial and temporal resolution product, this is equivalent to many signals, each of which has some uncertainty 455 (noise). Averaging results spatially or temporally, if possible for a downstream use of this product, will increase the signal to noise ratio and reduce uncertainty.
For example, because we provide basins for the entire ice sheet, total discharge is not subject to basin uncertainty. Any error in the delineation of one basin must necessarily be corrected by the inclusion (if underestimate) or exclusion (if overestimate) of a neighboring basin, although neighboring basins may introduce their own errors. Therefore, summing basins reduces the error introduced by basin outline uncertainty, and should be done if a downstream product does not need an estimate of discharge from a single outlet. This feature is built-in to coastal outlet discharge which is not as sensitive to our routing algorithm as ice margin outlet discharge because most coast outlets include a range of upstream ice margin outlets (e.g. Fig. 7 v. 9). Conversely, at the ice margin, outlet location and discharge volume is more uncertain than at the land coast. However, most runoff is generated near the ice margin and as runoff approaches the margin, there are fewer opportunities to change outlet location (Fig.   465 2).
Our coverage algorithm (Sect 3.3) only fills in glaciated regions that have at least some RCM coverage. When working with basins that have glaciated areas and no RCM coverage as in the case for all four of the GEM outlets near Nuuk, discharge could be approximated by estimating discharge from the nearest covered glaciated area with a similar climatic environment.
Temporally, errors introduced by this study's assumption of instantaneous discharge can be reduced by summing or averaging 470 discharge over larger time periods, or applying a lag function to the time series as done here and in van As et al. (2017).
Although a given volume of water may remain in storage long term, if one assumes that storage is roughly steady state, then long-term storage shown by, for example, dye trace studies, can be ignored -the volume with the dye may be stored, but a similar volume should be discharged in its place.

Quality control 475
The scale of the data are such that manual editing to remove artifacts is time and cost prohibitive. Here we provide one example of incorrect metadata. The elevation of each outlet is included as metadata by looking up the bed elevation in the BedMachine data set at the location of each outlet. Errors in BedMachine or in the outlet location (defined by the GIMP ocean mask) introduce errors in outlet elevation.
A large basin in NW Greenland has metadata outlet elevation > 0 (gray in Fig. 1) but appears to be marine terminating when 480 viewed in satellite imagery. Elsewhere the land vs. marine terminating color coding in Fig. 1 appears to be mostly correct, but this view only provides information about the sign of the elevation, not the magnitude (i.e. if the reported depth is correct). Ice outlets can occur above, at, or below 0 m. It is easier to validate the land terminating basins, which should in theory all have an outlet elevation of 0 m. That is not the case (Fig. 18). It is possible for land outlets to be correctly assigned an elevation > 0 m, if a land basin outlet occurs at a waterfall off a cliff (as might occur the edges of Petermann fjord) or due to DEM 485 discretisation of steep cells. However, most of the land outlets at elevations other than 0 are likely due to mask misalignment placing a section of coastline in a fjord (negative land elevation) or inland (positive land elevation). The bulk of land discharge (75 %) occurs within 0 ±10 m elevation, and 90 % within 0 ±30 m elevation (Fig. 18).

Other sources of freshwater
The liquid water discharge product provided here is only one source of freshwater that leaves the ice sheet and affects fjords 490 and coastal seas. The other primary freshwater source is iceberg calving and submarine melt at the ice/ocean boundary of marine terminating glaciers. A companion to the liquid water discharge product introduced here is provided by Mankoff et al. topography in BedMachine (i.e. correct basins modeled in this work) -here only Leverett Glacier (Sect. 4.2.3 & Fig. 9) 525 meets this criterion.
-For land basins, subglacial routing errors no longer exist, basins are well-defined, and errors are due to neglecting runoff delays or the RCM estimates of runoff.
-For ice basins, errors are dominated by basin uncertainty. Errors between similar-sized and neighboring basins are likely 530 to offset and may even cancel each other. Even so, a conservative treatment might consider errors between basins as random and reduce by the sum of the squares when summing discharge from multiple similar-sized and neighboring basins.

Product description
These data contain a static map of Greenland's hydrological outlets, basins, and streams and a times-series of discharge from 535 each outlet.
The output data are provided in the following formats: The stream data are provided as a GeoPackage standard GIS product and a metadata CSV that includes the stream type (start or intermediate segment), network, stream along-flow length, stream straight length, sinuosity, source elevation, outlet elevation, and a variety of stream indices such as the Strahler, Horton, Shreve, Hack, and other parameters (Jasiewicz and Metz, 2011). We note that the subglacial streams are unvalidated with respect to The data can be accessed with custom code from the raw data files. However, to support downstream users we provide a tool to access the outlets, basins, and discharge for any region of interest (ROI). The ROI can be a point, a list describing a polygon, or a file, with units in longitude, latitude (EPSG:4326) or meters (EPSG:3413). If the ROI includes any land basins, an option can be set to include all upstream ice basins and outlets, if they exist. The script can be called from the command line 555 (CLI) and returns CSV formatted tables, or within Python and returns standard Python data structures (from the GeoPandas or xarray package).
For example, to query for discharge at one point (50.5°W, 67.2°N), the following command is issued: python ./discharge.py --base ./freshwater --roi=-50.5,67.2 --discharge, where discharge.py is the provided script, ./freshwater is the folder containing the downloaded data, and --discharge 560 tells the program to return RCM discharge (as opposed to --outlets which would return basin and outlet information). The program documentation and usage examples are available at http://github.com/mankoff/freshwater (Mankoff, 2020b).
Because the --upstream option is not set, the --discharge option is set, and the point is over land, the results of this If the upstream option is set, two additional columns are added: One for each of the two RCM ice domains. A maximum of six columns may be returned: 2 RCM times (1 land + 1 ice + 1 upstream ice domain), because results are summed across all outlets within each domain when the script is called from the command line (summing is not done when the script is accessed from within Python).

570
If the --outlets option is set instead of the --discharge option, then results are a If the script is accessed from within Python, then the discharge option returns an xarray Dataset of discharge, without aggregating by outlet, and the outlets option returns a GeoPandas GeoDataFrame, and includes the geospatial location of all outlets and outline of all basins, and can be saved to GIS-standard file formats for further analysis.

Conclusions
We provide a 100 m spatial resolution data set of streams, outlets, and basins, and a 1 day temporal resolution data set of 580 discharge through those outlets for the entire ice-sheet area from 1958 through 2019. Access to this database is made simple for non-specialists with a Python script. Comparing the two RCM-derived discharge products to 10 gauged streams shows the uncertainty is approximately plus-or-minus a factor of five, or half an order-of-magnitude, or +500%/-80%, when comparing daily discharge for single or few basins.
Because of the high spatial (individual basins) and temporal (daily) resolution, larger uncertainty exists than when working 585 over larger areas or time-steps. These larger areas and times can be achieved through spatial/temporal aggregating or by implementing a lag function.
This liquid freshwater volumetric flow rate product is complemented by a solid ice discharge product .
Combined, these provide an estimate of the majority of freshwater (total solid ice and liquid) flow rates from the Greenland ice sheet into fjords and coastal seas, at high temporal resolution and process-level spatial resolution (i.e. glacier terminus for 590 solid ice discharge, stream for liquid discharge).
This estimate of freshwater volume flow rate into Greenland fjords aims to support further studies of the impact of freshwater on ocean physical, chemical, and biological properties; fjord nutrient, sediment, and ecosystems; and larger societal impacts of freshwater on the fjord and surrounding environments.
7 Code and data availability 595 The data from this work is available at doi:10.22008/promice/freshwater (Mankoff, 2020a).
The code and a website for post-publication updates is available at https://github.com/mankoff/freshwater (Mankoff, 2020b) where we document changes to this work and use the GitHub Issues feature to collect suggested improvements, document those improvements as they are implemented, document problems that made it through review, and mention related. This version of the document is generated with git commit version 831ce52 .  Table 1). Grey band shows 5 to 95 % prediction interval. Red band shows 5 to 95 % prediction interval when removing the GEM stations near Nuuk (Table 1) that have small glaciers not included in the RCMs (5341 days remain).    (Jordahl et al., 2020), numpy (Oliphant, 2006), x-array (Hoyer and Hamman, 2017), and Matplotlib (Hunter, 2007) packages. The entire work was performed in Emacs (Stallman, 1981) using Org Mode All code used in this work is available in the Supplemental Online Material.
Author contributions. KDM produced this work -wrote the code and the text. XF and BN supplied RCM inputs. APA, WIC, DVA, and RSF helped with discussions of methods, quality control, or writing. KK and SS supplied Qaanaaq data. KL provided GEM data.