{'fraction_non_alphanumeric': 0.03237241820733118, 'fraction_numerical': 0.026858029098249674, 'mean_word_length': 4.655297229569817, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Increasing surface runoff from Greenland’s firn areas

At high elevations of ice sheets, melting snow generally percolates and refreezes, so does not contribute to the shrinking of the ice sheet. Here, we systematically map the runoff area of the Greenland ice sheet, using surface rivers visible on satellite imagery. Between 1985 and 2020, the maximum runoff elevation rose by 58−329 metres, expanding the runoff area by 29% (-8%/+6%). Excess melt beyond the refreezing capacity of pores in snowfall has created near-impermeable ice slabs that sustain surface runoff even in cooler summers. We show that two surface mass balance models over-estimate the runoff area by 16−30%. Once restricted to our observed areas they indicate that 5−10% of recent runoff likely comes from the expanded runoff area. Runoff from higher elevations is sensitive to projected warming as further increases in the runoff limit will increase the runoff area disproportionately.

Two end-members of Greenland Ice Sheet firn behaviour under exposure to large amounts of meltwater have been identified. Firn aquifers, where meltwater percolates to several metres below the ice surface and is stored as a liquid or discharged 15,16 , are found in regions of high snowfall and are principally located along the steep south-east margins [17][18][19] (Fig. 1a). Towards the other end of the spectrum, in regions with less snowfall, percolating meltwater increasingly refreezes in firn pore space, initially forming ice lenses several centimetres thick during individual melt seasons 21,22 . When surface melting consistently exceeds the replenishment of pore space by accumulation, ice lenses can fuse together into low-permeability ice layers tens of centimetres thick, and subsequently homogeneous ice slabs several metres thick and kilometres wide, rendering the underlying pore space essentially unavailable for further refreezing 6 . Modelling indicates that while around 40 % of meltwater refroze in the years up to 2011, this buffering capacity might decrease rapidly through the 21 st century 23 . By reducing the buffering capacity of firn, ice slabs may raise the runoff limit to higher elevations.
Firn cores and ground penetrating radar measurements along the 'EGIG' line (70.5°N; Fig.  1) during the 2000s show several ice layers at depth 24,25 . Cores acquired in 2012 and 2013 at 1,840−2,360 m asl on the K-Transect slightly farther south ( Fig. 1) revealed that the buffering capacity of firn at up to ~1,900 m asl was near-exhausted due to the presence of metres-thick, contiguous near-surface ice slabs, rendering deep pore space difficult to access 14 . Furthermore, optical satellite imagery from summer 2012 show that surface rivers caused visible runoff from ~1,850 m asl in this region for the first time in the observational record. At an ice-sheet-wide scale, airborne radar observations made in 2010−2014 indicated that metres-thick ice slabs were widespread along the western and northern margins of the ice sheet 26 , but did not show whether runoff had occurred.
Here we quantify changes in the visible runoff limit (hereafter simply 'runoff limit') around the Greenland Ice Sheet by using surface hydrological features in remotely-sensed imagery as indicators of runoff (see Methods). Specifically, open river channel flow is an efficient means of evacuating water from the ice sheet surface, with stream velocities of 0.4 to 2.6 m s -1 in the ablation zone 27 . Networks of rivers develop over the ice sheet surface each summer, sourcing water from the inter-channel regions, and can flow for tens of kilometres before draining into terminal moulins 28 or directly off the ice. We processed 35 years of Landsat near-infrared 30 m imagery acquired each July and August from 1985 to 2020, totalling >25,000 scenes. We examined all regions where we would expect visible surface runoff to occur: the west margins north of 61.5 °N, the north margins, and the east margins north of 75°N. Conversely, in the high-snowfall south and south-east regions of the ice sheet, sub-surface firn aquifer systems dominate 18 , while the central-east margins are topographically complex and rise more steeply than the west (Fig. 1), inhibiting surface river development that would be visible in Landsat imagery.
We identified hydrological features using a sequence of high-pass Discrete Fourier Transform filtering, directional filtering, thresholding and morphological operations adapted from 29 (Extended Data Figs. 1 and 2). We then used the highest-elevation hydrological features to define the position of the runoff limit at Landsat's native 30 m resolution. We aggregated the retrievals to identify the annual runoff limit at ~1 km intervals. Exploiting the high spatial density of our observations, we express the uncertainty of each annual runoff limit elevation as the Median Absolute Deviation (MAD) of the observation cluster to which it belongs (Methods and Extended Data Fig. 3). We then mapped the corresponding annual runoff area ±1 MAD. To examine trends in runoff limit behaviour we calculated means over 100 km-wide zones and by geographic region. To examine net ice-sheet-wide changes we computed the mean runoff limits and areas across two multi-year reference periods: 1985−1992 and 2013−2020. During these periods we also manually mapped the runoff limits of two small regions in the NW and NE where rapid changes in slope direction prevented automatic identification (see Methods).

Rising visible runoff limits
Between 1985−1992 and 2013−2020, the runoff limit rose along the west margins of the ice sheet (on average by 242 m) and along most of the north (194 m) and north-east (59 m) margins (Fig. 1a). The area increased by between 21% (−1MAD) and 35% (+1MAD), with a non-centred median of 29% (~47,400 km 2 ) (Fig. 1b). Landsat acquisitions are intermittent, which may lead to the maximum runoff limit being missed in space and/or time. Nonetheless, the expected runoff limit elevation was observed multiple times by Landsat in most summers except for an ice-sheet-wide data gap in 1996 (see Methods and Extended Data Fig. 4). We further estimated each year's maximum likely runoff area using the mean seasonality of observed runoff limit behaviour (see Methods and Extended Data Fig. 5). The maximum likely runoff area falls within the ±1MAD bounds of the size of the increased runoff area (Fig. 1b). Overall, there were no net decreases in the runoff limit. From 1985 to 2020, the S, SSW and SW saw the largest sustained runoff limit rises (9.2−10.6 m a -1 ; Fig. 1d), while the rise in the NE was much smaller than other regions (2.4 m a -1 ).
During the late 1980s and 1990s, the runoff limit often varied substantially between years (Fig. 1c). This variability is especially clear farther north where more satellite observations are available (Extended Data Fig. 4). From the 2000s onward, the runoff limits rose and inter-annual variability (i.e. the standard deviation of the runoff limits) reduced on average by 20 m during 2000−2020 compared to 1985−1999, except in the S. Within each region, the runoff limit was generally consistent across all of the 100 km zones, indicating that the response of firn to meltwater production is relatively homogeneous at regional scales.

Sustained excess melting enables runoff
In light of the substantial rise in the runoff limit since 1985, we examined the melt conditions and their probable impacts on near-surface firn density. We considered the ability of pore space within annual accumulation to accommodate summer melt volumes. Excess melt (M e , Methods eq. 1) is the annual amount of meltwater which exceeds the storage capacity of the pore spaces in that year's accumulation by percolation and refreezing 7,26 . Conceptually, if M e is positive but firn pore spaces below that year's accumulation are accessible then meltwater will percolate into those pore spaces, densifying the underlying firn. Relatively impermeable ice layers will develop, inhibiting further percolation. Thus, sustained positive M e is indicative of firn with a densified stratigraphy conducive to surface runoff. We calculated M e over the the expanded runoff area (pink in Fig. 1a) using accumulation, melting, rainfall and air temperatures from the regional climate model RACMO 31 (see  Methods). M e increased on average by 174 mm w.e. a -1 between 1985−1992 and 2013−2020, equivalent to +34%, varying from +9% (S) to +102% (NE). M e shows no simple relationship with the runoff limit. We therefore searched for points in time when the relationship between Me and the runoff limit changed. We found significant break-years (p <0.01) in all regions ( Fig. 2a; Extended Data Fig. 6), which we used to split each region's time series into 'early' and 'late' periods (light and dark shading respectively in Fig. 2b-h).
Next we fitted least-squares regression models between M e and the runoff limit for each period, which are generally significant (p <0.05) although not in the N and NE (Fig. 2b-h).
After the break-year, the runoff limit depended less on M e (denoted by shallower regression gradients in most regions) and reached higher elevations even when M e was as low as typical 1980s/1990s amounts, especially in the SW, NW and NNW regions. In the SW, the 2003 break (Fig. 2a,d) corresponds with an abrupt 80% increase in proglacial runoff which has previously been attributed to a shift in atmospheric circulation 32 ; the high runoff limit probably amplified these atmospheric effects.
We assessed the importance of sustained positive M e in enabling higher runoff limits by examining 10-year averaged excess melt M e 10y 26 . In most regions, higher runoff limits only occurred once M e 10y had increased (shown by circle sizes in Fig. 2b-h), despite previous single years of high M e . This indicates that runoff from higher elevations is unlikely to occur solely in response to high M e . Conversely, once M e 10y increases then higher runoff limits can be sustained even in individual years of low M e . We interpret these findings as evidence that sustained high M e enables firn densification, causing relatively impermeable, metres-thick near-surface ice slabs to develop as mapped previously 26 . The expanded runoff area corresponds strongly with the locations of these ice slabs (Extended Data Fig. 7). Lower runoff limits observed during the 1980s and 1990s suggest that ice slabs developed as recently as the 2000s. Cores from the K-Transect show that once ice slabs have developed then they continue to thicken during meltwater infiltration and refreezing over several years of below or at median melt 33 , indicating that, once formed, they are relatively resilient features. Furthermore, these cores indicate that anomalous firn densification has already commenced above the present visible runoff limits 14 , priming the sub-surface for further increases in the runoff limits.
Our results suggest that northern regions, particularly the NE, are primed for large increases in the runoff limit. NE M e 10y has increased since the break-year of 2007, so far with relatively small impacts on the runoff limit, but airborne radar surveys of this area already show evidence of firn densification and formation of ice slabs extending beyond the expanded runoff area (Extended Data Fig. 7). Farther north, firn temperatures are colder (~-24 °C at Camp Century, 1,886 m asl (Fig. 1a) 34 compared with ~-10 °C at KAN_U on the K-Transect, 67°N, 1,840 m asl 14 ), causing water to refreeze at shallower depths 22 . This will enable ice slabs in the north to develop with less M e 10y than farther south. Our interpretation suggests that northern regions have less buffering capability against excess melting, a finding complemented by surface mass balance modelling which suggests that firn refreezing capacity in the NE has declined sharply in recent years 31 . Higher M e 10y in the north during the last decade suggests that these margins are committed to rising runoff limits in the coming years.

Models show large differences in runoff
We examined how two regional climate models which calculate the ice sheet's surface mass balance, MAR v3.11 35 and RACMO2.3p2 31 , perform in the expanded runoff area. Both models include multi-layer snow schemes which represent firn densification and distinguish between retention by refreezing versus runoff. These models have two particular limitations relevant to our comparison with the observed runoff limit: (i) they are not able to model lateral flow, and (ii) they do not distinguish between sub-surface runoff versus that which occurs over the surface. Our observations therefore cannot be compared exactly to modelled runoff area, but nevertheless enable a valuable assessment of model performance in these areas where in-situ observations are rare to non-existent 36 .
We estimated the modelled annual runoff area by applying a threshold of 10 mm w.e. a -1 (chosen following a sensitivity analysis; see Methods and Extended Data Fig. 8). Both models almost always yield larger runoff areas than observed: on average, MAR by 30% and RACMO by 16% ( Fig. 3a; region-specific comparisons in Extended Data Fig. 9). This means that the models produce runoff over the entire observed runoff area and beyond. Whilst MAR's runoff area trend relative to observations is essentially stationary, RACMO's is non-stationary, showing better agreement with observed runoff area from 2002 onwards.
The causes of this require further examination.
Next, we calculated first-order estimates of the modelled runoff volumes from above the 1985−1992 runoff limit. We calculated annual runoff (i) unbounded by our annual observations, and (ii) up to each annual runoff limit ±1MAD. When unbounded by our observations, cumulative runoff from our study regions since 1985 was 548 Gt according to RACMO and 941 Gt according to MAR. The unbounded volumes are especially large in major melt seasons such as 2012 and 2019; in these years, modelled runoff volumes from above the 1985-1992 runoff limit are often very high, peaking at as much as 20 Gt in the SW (Extended Data Fig. 9).
When we considered only the runoff area up to our observed limits ( Fig.3b-d) then cumulative runoff since 1985 fell to 232 Gt (RACMO) and 422 Gt (MAR), a factor-two reduction. Considering our observed area − 1MAD, cumulative runoff volumes are 190 and 320 Gt respectively, while at +1MAD they are 264 and 491 Gt, so uncertainty in modelled runoff volumes caused by the modelled runoff area is much greater than that originating from our observations. Nonetheless, the expanded runoff area has clearly been an emergent source of mass loss since the early 2000s (Fig. 3c,d).
Our comparison of the observed and modelled runoff areas shows that total modelled runoff volumes from the percolation zone are almost certainly over-estimated, particularly by MAR. It is likely that our observed runoff limits are indicative of the true runoff limit. First, field measurements of sub-surface matrix runoff in firn 37 and firn aquifers 15 suggest that sub-surface runoff is limited to within a few kilometres of the visible runoff limit (see Methods). Second, our observations agree well with field observations of meltwater retention by thermal tracking in several boreholes on the EGIG line 21 (Fig. 1). The boreholes showed that the transition from full meltwater retention to full runoff occurs over a zone. 21 did not quantify the change in retention through the zone, but they did conclude that during 2007 and 2008 the runoff zone spanned ~1,350 to ~1,500 m. Our observations show that during 2007 and 2008 the runoff limit in this region reached 1479−1490 m. They therefore agree well with the upper limit of the runoff transition zone estimated from thermal tracking, which provides further support that the actual runoff limit is located near to the visible runoff limit.
Overall, our observations reveal that Greenland's visible runoff limit rose by on average 194 m between 1985−1992 and 2013−2020. Rising runoff limits were preceded by sustained high excess melting, implicating recent ice slabbing in driving the runoff limit upward. Runoff now regularly occurs from higher elevations even in years of relatively little excess melting.
With our observations we showed that surface mass balance models overestimate highelevation runoff. Considering only runoff from the expanded visible runoff area, recent contributions likely total 12−15 Gt a -1 (RACMO) or 20-28 Gt a -1 (MAR). Ice-sheet-wide mass losses during 2012−2017 were 244 Gt a -12 , of which our estimates constitute 5 or 10%. This is similar to 11±4% of runoff originating from the expanded runoff area in the K-Transect region during 2012 according to a firn core analysis 14 .
The large uncertainties in modelled runoff from higher elevations show that modelling of melting and firn processes in these areas requires improvement ahead of projected warming. The rise in the runoff limit elevation since 1985 has transitioned an additional ~13% of our study area between 500 and 2200 m to runoff. Another rise of the same magnitude would see a further ~17% of this area transition to runoff, highlighting that the ice sheet's hypsometry 38 will cause future rises to yield disproportionately larger runoff volumes.

Extracting visible runoff limits from Landsat imagery
We use surface hydrological networks as indicators of the location of the visible runoff limit. In the percolation zone, rivers and saturated slush fields are sufficiently large to be detectable in Landsat  We acquired Landsat Level-1 imagery from the public Google Cloud Landsat dataset. We processed all images acquired over the ice sheet margin areas of interest during July and August each year with a Landsat Level-1 product 'Scene Cloud Cover' score of less than 50%, yielding >~25,000 scenes. East margins south of 75°N were excluded from analysis on the basis of (a) previously-mapped aquifers 18 , (b) steep surface slopes with crevassing preventing widespread surface drainage networks from developing, and (c) small melt areas mainly restricted to topographically complex valley glaciers. We gap-filled Landsat 7 OLI SLC-off data by interpolation using the 'fill no-data' functionality in the Python Rasterio library 40 .

Imagery pre-processing
To detect hydrological networks and particularly surface channels in NIR images, we characterise them by their Gaussian-like cross-sections and longitudinal continuity following approaches proposed by 29,41 . Here we outline our processing workflow developed within the Python 'eo-learn' library (Sentinel Hub, Sinergise) and which we apply to each Landsat image (numbered steps hereafter correspond to those in Extended Data Fig. 1), commencing with the Level-1 NIR imagery (1). (2) For TM imagery only we apply a 3×3 pixel median filter to reduce noise and speckling. (3) We apply spectral discrimination with a high-pass Discrete Fourier Transform (DFT) filter which ramps up between 1/10 and 13/1000 m -1 , delineating active rivers as dark, linear features with abrupt bright channel banks that manifest as high-frequency information. (4) As noise tends to remain following DFT filtering, we denoise with a non-local means filter with a 150×150 m search window and 60×60 m patch size.
We then proceed to hydrological network identification. (5) We use a 12-direction Gabor filter with a kernel width of 1.5 pixels to preferentially retain linear features. (6) We apply a path opening filter with length 1.2 km (40 pixels) to lengthen river channel continuity and further suppress noise. (7) Finally, we calculate a global threshold as T = x + 0.5σ (where x is the image mean and σ is the image standard deviation) to produce a binary hydrology map of each image.
We apply a cloud mask to the binary hydrology map by additively combining the Landsat cloud and cloud shadow masks from the Quality Assessment product with additional estimates of cloud cover based on Landsat shortwave-infrared reflectance (R SW ). We apply a threshold of R SW > 0.08 then dilate the result by 1.5 km. This approach provides additional redundancy in cloud masking without excessively reducing data coverage. Its application during the final stage of pre-processing prevents false-positive hydrological features from being generated along 'high frequency' mask edges.
To verify the suitability of 30 m NIR imagery for extracting high-elevation hydrology maps we ran our pipeline on both NIR and 15 m panchromatic bands for several Landsat 8 scenes. We found no major differences between the maps yielded by these two resolutions. Previously, 28 found that the drainage densities of river networks in coincident Landsat panchromatic and WorldView-2 0.5 m panchromatic imagery show good agreement, and indicate that Landsat panchromatic data discriminate rivers especially well at higher elevations. WorldView-1 0.5 m imagery acquired around KAN_U indicates that rivers can grow to at least 10 m wide less than 1 km from their channel heads 14 , so it is very unlikely that rivers can extend long distances while maintaining sizes below our detection threshold. Finally, our use of the same spatial resolution throughout the time series ensures that we always capture the same minimum size of hydrological feature. Thus, while our use of 30 m imagery may miss the smallest river channels, we are confident that the majority of the active channel network is mapped consistently.

Identification of visible runoff limits
Following image pre-processing, candidate hydrological features such as small isolated supraglacial lakes remain which are not representative of the runoff limit. We therefore dilate the binary hydrology mask by 200 m to consolidate feature continuity, then use image object segmentation to label each contiguous candidate (step 7 in Extended Data Fig. 1). We remove candidate features if they are (a) located at elevations greater than 2,100 m according to GIMPDEM 20 (the maximum elevation at which we manually identified connected hydrological features in recent satellite images, as opposed to isolated supraglacial lakes), (b) have an area less than 0.5 km 2 , (c) have a minimum NIR reflectance >0.2 (indicating that there are no water pixels present in the feature), (d) when skeletonized (i.e. reduced to their medial axis as a proxy for feature length) have a length less than ~1.5 km, or (e) are near-perfectly circular (which would be strongly indicative of a supraglacial lake). We also require features located at elevations greater than the 90th percentile of all features in the image to have an area of at least 1.5 km 2 . This yields a filtered set of features (corresponding to step 8 in Extended Data Fig. 1).
Next, the filtered image is rotated in the direction of the surface slope derived from GIMPDEM 20 , such that the elevation decreases from left to right. Then, in each row of the rotated image, the leading edge of the first candidate hydrological feature encountered is picked as the location of the runoff limit, yielding a set of up to ~ 8,000 (x, y, z) coordinates which describe the position and elevation of the runoff limit for the image under examination. By applying this process to >25,000 images we were able to make >33 M individual retrievals.

Annual visible runoff elevation
We used the multiple Landsat images acquired each melt season to derive robust annual estimates of the runoff limit. This occurred in three principal phases.
First, we stacked all the filtered candidate hydrological features identified across the whole time series to produce a reference map at 30 m resolution of hydrological feature incidence. This takes advantage of the propensity for hydrological features to re-occur in the same locations each melt season. Following manual visualisation, we converted the reference map to a binary representation of hydrological features using a minimum threshold of 65 occurences of a hydrological feature in a given pixel to be retained, which ensured that the highest contiguous hydrological features recorded in extreme melt years such as 2012 were included. Next, we filled holes in the map whose areas were less than 1 km 2 , then skeletonized all retained hydrological features to generate a proxy hydrological channel network. We generated a K-d nearest-neighbour lookup tree of coordinates corresponding to all the features in the proxy network, which we then used to filter out any individual runoff limit retrievals located >3 km away from the proxy hydrological network.
In the second phase, we divided the ice sheet margins into ~100 km-wide zones oriented approximately parallel to the surface slope. Each zone is 100 km wide only at 2,000 m asl, widening or narrowing towards lower elevations in order to conform to the ice sheet margin geometry. On an annual basis, in each zone we used K-Means clustering to group runoff limit retrievals together according to their elevation and date of measurement (Extended Data Fig. 3a). This enabled us to estimate measurement uncertainty as a function of cluster median absolute deviation (MAD). We initialised the K-Means algorithm with n clusters, where n denotes the number of discrete 7-day periods during which observations were made in July and August, e.g. if observations were made in the first week of July and the last week of August then two clusters were specified. All individual runoff limit retrievals were labelled with their cluster membership (Extended Data Fig. 3b).
In the third phase, we defined the runoff limit each year using our cluster-labelled retrievals. First, we split each zone into 1 km-wide slices (S) (Extended Data Fig. 3c). In each S, we examined all individual runoff limit retrievals made each year. We identified the cluster in the S with the highest median elevation (Extended Data Fig. 3d) then used the coordinates of this median elevation to define the annual visible runoff limit of the slice, Ύ S , with the uncertainty prescribed as the cluster's MAD.
The runoff limit in each 100 km zone was calculated as the median of all S, i.e. median (Ύ S1 … Ύ S100 ) and corresponds to the values in Fig. 1c. Regional annual runoff limits were calculated as the mean of all member zones and are shown in Fig. 1d. Observations along all margins during 1996 and in the SSW and S prior to 1991 were very sparse (Extended Data Fig. 4) so we excluded these observations from further analysis.
Where topography and slope direction change rapidly through space, our automatic workflow could not reliably find the runoff limit. The runoff limits of two small areas therefore had to be mapped manually: (1) between the NNW and N regions, and (2) between the southern limit of the NE region and 75°N (Fig. 1a). In these areas we mapped the runoff limit during each of our two reference periods using Landsat NIR 25 th percentile reflectance mosaics created from all July-August images collected during each reference period. We generated these mosaics using Google Earth Engine (Google, Inc.). Then, in each mosaic, we manually delineated the runoff limit by selecting the locations in each image which corresponded to the channel heads of all visible hydrological features.

Estimating the annual maximum runoff limit
Landsat imagery acquisitions are temporally sparse. The satellite orbit repeat time is several days and cloud cover often prevents observation of the ice sheet surface. Thus, these relatively sparse observations may not always capture the true (highest) annual runoff limit.
Conversely, 36 years of acquisitions can be combined to yield an average seasonal elevation change in the runoff limit through July and August. In each region we first normalised all Ύ S retrievals by subtracting the annual runoff limit, then we combined all 36 years of data into a seasonal series by calculating the mean daily relative runoff limit on each day of July and August. Next we fitted a 2 nd degree polynomial curve which describes the relative runoff limit according to the day of year (Extended Data Fig. 5, right panels). In the NW we fitted the polynomial curve only to observations up to day 230 to remove the impact of late-season outliers.
We used these region-specific polynomial curves to estimate the maximum likely annual runoff limit in each member S (see example in Extended Data Fig. 3e). We first set the zeroelevation-change point of the polynomial curve to correspond to the date of the maxmimum observed runoff limit, then we calculated the maximum likely runoff limit by adding the residual change-in-elevation to the observed runoff limit. Thus, in years when observations are made very close in time to the maximum annual runoff limit then the maximum likely runoff limit will nearly equal the maximum observed runoff limit. Conversely, if the highest observation is some time away from the seasonal runoff limit maximum then the maximum likely runoff limit can be substantially higher than that year's observation.
On an annual basis and by region (Extended Data Fig. 5), we see no systemmatic evidence of temporal trends in the difference between the observed runoff limit and the maximum likely runoff limit that might otherwise invalidate our trend analysis (Fig. 1d). Differences are lowest farther north. In the NW, surface slopes are generally steeper in the vicinity of the runoff limit compared to other regions, which manifests as greater differences between the observed and maximum likely runoff limits. The small S region displays the largest differences and are explicable by less dense observation coverage (Extended Data Fig. 4a).

Calculation of visible runoff area
We distinguish between calculating (i) the runoff area from observations of the runoff limit, which is straightforward, versus (ii) the areas which correspond to the runoff limit ±1MAD and the maximum likely runoff limit, which requires interim steps to convert the elevation range into (x, y) coordinates.
In the second case, we utilised 1-D centre-line elevation profiles of each S extracted from GIMPDEM. For each Ύ S we looked up the nearest (x, y) coordinate corresponding to the runoff limit ±1MAD or the maximum likely runoff limit. Finally, we converted these coordinates to runoff area polygons using the same spline approach outlined above.
We note that in moving from the runoff limit ±1 MAD to the runoff area ±1 MAD, we implicitly introduce the effect of ice sheet hypsometry, which means that unlike elevation the area MAD is not equal in both directions.

Observation availability by elevation
We checked for the possibility that the maximum observed runoff limit is an artefact of the maximum elevation at which clear-sky imagery were available. For each year, at each 30 m pixel on a common Polar Stereographic projection we counted the number of cloud-free observations made, then aggregated the counts to 1 km resolution using mean resampling. We divided each region into 100 m elevation bins and for each bin we calculated the mean number of observations made each July-August, to be compared against the identified runoff limit (Extended Data Fig. 4). We see no evidence that the observed runoff limit is systematically affected by sparse data availability at higher elevations.

Excess melting
We calculated the amount of annual excess melt M e in mm water equivalent (w.e.) following 7 modified by 26 to include rainwater: where L is the latent heat refreezing capacity of ice (J kg -1 ). Values for accumulation (C, kg m -2 ), melting (M, kg m -2 ), firn temperature, derived from the mean annual air temperature (Tf, °C) and rainwater (R, kg m -2 ) were taken from the RCMs. We prescribed the density of refrozen ice (ρ r ) as 873 kg m -3 . The density of fresh snow accumulation (ρ c ) was calculated according to a geographically-based parameterisation used in surface mass balance models 42 and ranged from 300 to 380 kg m -3 . We first calculated M e for each RCM cell. Then, in each region, we calculated mean annual M e in the area between the 1985−1992 and 2013−2020 runoff limits. We also applied a 10-year rolling mean to the M e time series in each region to yield M e 10y, such that a value for 2020 is the mean of values from 2011−2020 inclusive. Fig. 2 shows M e calculated from RACMO, in recognition of its closer correspondence with the observed runoff area compared to MAR. However, we note that the break-years and regression coefficients calculated using MAR M e are sufficiently similar to those from RACMO that our interpretation would not change if we used MAR-derived values instead.

Break-years in visible runoff limit behaviour
We undertook a break-point analysis on each runoff limit time series in each region. First, we split each time series (Fig. 1d) into two for all possible break-years in the range (1995, 2010). We checked each split time series for normality with a two-tailed D'Agostino-Pearson test (Extended Data Fig. 6a,b), under the null hypothesis (H 0 ) that each time series was normal. We accepted H 0 when p > 0.05. Next, we used two-tailed t-tests to evalulate the significance of each possible break-year (Extended Data Fig. 6c), under H 0 that there is no difference between the two populations divided by the chosen break-year. We rejected H 0 when p < 0.05. Each region therefore had several possible significant break-years. We chose the year with the highest significant T-statistic (Fig. 2a).
To quantify the different behaviours before versus after the break-year we under-took ordinary least-squares regression on each part of the time series (Fig. 2b-h).

Runoff volumes
We estimated annual runoff volumes from each region using modelled runoff output by the regional climate models (RCMs) (i) MAR v3.11 forced with ERA-5 43 , run at 10 km resolution 35 , and (ii) RACMO2.3p2 forced with ERA-Interim 44 (1985−1989) and ERA-5 (1990ERA-5 ( −2020, run at 5.5 km resolution then statistically downscaled to 1 km resolution 31 . Recall that our observations of the runoff area quantify the visible extent of runoff, whereas RCMs (1) do not model lateral flow between adjacent grid cells and (2) do not explicitly distinguish between the portion of runoff visible on the surface versus that which initially runs off via the sub-surface before entering the visible surface drainage network. Whilst meltwater flow through snow has been relatively well studied (e.g. 45,46 ), flow through firn has received limited attention. Furthermore, the boundary conditions which govern the lateral transition from sub-surface (matrix) to visible (open channel) flow remain unknown. A direct comparison between observed and modelled runoff area is therefore not a realistic goal. Nevertheless, while the RCM runoff schemes are relatively simple, it is important that they approximately capture the expansion of the runoff area captured by our visible runoff limits, which in the absence of observations of runoff volume from these areas provide an important means of model validation. Otherwise, future scenarios of Greenland mass balance would be subject to further uncertainties.
We applied a series of thresholds to the RCM runoff fields, ranging from 1 to 20 mm w.e. a -1 (Extended Data Fig. 8). This threshold range is designed to check the sensitivity of total runoff area computed from RCM runoff to relatively small but non-zero quantities of runoff. It is not necessarily indicative of RCM visible runoff area, which as highlighted above would require explicit consideration of (i) upstream contributing area to define the flux of runoff through any one model cell and (ii) whether that runoff flux becomes visible at the surface, for instance as a function of the sub-surface water table height.
RACMO is more sensitive than MAR to the selected threshold, with ~90−95 % of the 1 mm w.e. a -1 area retained at 20 mm w.e. a -1 compared to ~98 % for MAR (Extended Data Fig. 8). This is partly a function of RACMO's 1 km spatial resolution, which permits a more spatially gradated reduction in runoff quantity than is possible at the 10 km resolution of MAR. Nonetheless, the runoff area simulated by MAR for the same threshold is substantially larger than RACMO's, exceeding the uncertainty in the area calculation associated with the choice of threshold. We adopted 10 mm w.e. a -1 as the threshold underlying the data shown in Fig. 3a as a compromise between screening out areas with a very low runoff contribution (at 10 mm w.e. a -1 , a 10 km 2 grid cell yields only ~0.001 Gt a -1 ) versus the conceptual basis that all runoff modelled by RCMs is assumed to exit the ice sheet.
We estimated the contribution of the expanded runoff area by integrating RCM runoff occurring above the 1985−1992 runoff limit bounded by each year's runoff limit ±1MAD (Fig. 3b,c). For comparison purposes, we also computed the maximum runoff contribution from above the 1985-1992, ignoring the observed runoff limit (Extended Data Fig. 9). MAR has a relatively coarse spatial resolution compared to the variability of the runoff limit, so for the MAR cells intersecting the runoff limits we computed the runoff volume contribution from the pixel's fractional area that was located above the 1985−1992 limit or below the specific year's limit. We applied a latitudinal scaling correction for area 47 prior to integration.
Given that our observations are of the visible runoff limit, we note that runoff volumes estimated from the runoff area are likely be conservative because they do not account for runoff from higher elevations. The runoff volumes estimated from between the 1985−1992 limit and the annual limits +1MAD could be more realistic estimates of total runoff because they implicitly capture sub-surface runoff. However, runoff which commences in the sub-surface above the visible runoff limit is unlikely to form a major contribution. For example, 37  circle shows the estimated maximum runoff limit in 2019 (see e,f). d Density histogram of all runoff limit retrievals made in SW-1 during 2019 (colours correspond to the clusters in panels a and c, and saturation shows observation density). The runoff limit picked for S in panel c is shown by the black circle. e The polynomial adjusted to intercept 2019's runoff limit (black circle and in panel c) estimates that the maximum likely runoff limit was slightly higher (red circle and in panel c).  for the part of the time series after the specified break date. c Two-tailed t-test of whether each runoff limit time series can be separated into early versus late periods by the specified break-year. p < 0.01 indicates that the break-year is significant. Each chosen break-year has the largest T-statistic at p < 0.01, and the early and late periods can be considered normal (as panels a and b show that H 0 could not be rejected).
Extended Data Figure 7. Correspondence of visible runoff limits with ice slab locations. 1985-1992 (blue) and 2013-2020 (pink) runoff areas overlaid with ice slab locations (black outlines) mapped by 26 using airborne Accumulation Radar. Grey lines show the flight lines of the Operation IceBridge flights on which the Accumulation Radar was flown.
Extended Data Figure 8. Estimates of runoff area from regional climate models.
The sensitivity of modelled runoff area depending on the chosen threshold value above which annual runoff contributes to runoff area, over the range 1 to 20 mm w.e. a -1 . Results in the main text use 10 mm w.e. a -1 . a, c Estimates for MAR. b, d Estimates for RACMO. a, b Runoff area at each threshold. c, d For the threshold values 5, 10 and 20 mm w.e. a -1 , the percentage difference to the area calculated at the 1 mm w.e. a -1 threshold. . c Elevation of the runoff limit in each 100 km zone, grouped by region. Hatching indicates less than 25% spatial coverage. d Regional runoff limits ±1MAD (solid lines and shading) and corresponding linear regression weighted by MAD (dashed lines and text labels).    a Comparison between the observed and modelled runoff areas. Modelled runoff area includes all pixels in the model domain which contribute >10 mm w.e. a -1 . Values >100% indicate that the modelled runoff area is larger than observed. b Annual runoff volume (gigatonnes, Gt) from above the 1985−1992 runoff limit up to each year's observed runoff limit. c, d Cumulative annual runoff volume since 1985 (Gt) from above the 1985-1992 runoff limit up to each year's observed runoff limit, by region, c by MAR and d by RACMO.
Tedstone and Machguth Page 28 Nat Clim Chang. Author manuscript; available in PMC 2022 December 16.