1. Introduction
Surface meltwater runoff from the Greenland ice sheet has increased since the early 1990s (van den Broeke and others, Reference van den Broeke2016). This is due to intensified melting (Fettweis and others, Reference Fettweis, Tedesco, van den Broeke and Ettema2011; Hall and others, Reference Hall, Comiso, DiGirolamo, Shuman, Box and Koenig2013; Trusel and others, Reference Trusel2018) but also because the area of the ice sheet drained by surface rivers increased by 29% between 1985 and 2020, with concomitant rises in the maximum visible runoff limit, that is the elevation up to which visible surface drainage networks develop (Tedstone and Machguth, Reference Tedstone and Machguth2022). The newly expanded runoff areas account for 5–10% of recent mass losses from the Greenland ice sheet and correspond strongly with where metres-thick low-permeability ice slabs are found in the firn (Tedstone and Machguth, Reference Tedstone and Machguth2022). MacFerrin and others Reference MacFerrin(2019) define an ice slab as an ice horizon which is thicker than 1 m, laterally widespread and contiguous (unlike lenses or layers). Slabs form when large amounts of meltwater percolate and refreeze in the firn, consistently exceeding the replenishment of pore space by accumulation over several years (de la Peña and others, Reference de la Peña2015; Machguth and others, Reference Machguth2016; MacFerrin and others, Reference MacFerrin2019). Slabs inhibit vertical percolation and refreezing, forcing surface meltwater to instead flow laterally over the slab and hastening the transition of percolation facies from a melt-retention regime to a melt-runoff regime (Machguth and others, Reference Machguth2016; Mikkelsen and others, Reference Mikkelsen2016). Supraglacial streams and rivers form on top of these ice slabs and transport meltwater downstream (Tedstone and Machguth, Reference Tedstone and Machguth2022), sometimes flowing tens of kilometres before they finally drain either into the ice or over the margin (Holmes, Reference Holmes1955; Poinar and others, Reference Poinar, Joughin, Das, Behn, Lenaerts and Broeke2015; Yang and Smith, Reference Yang and Smith2016).
Ice slabs are relatively widespread features on the Greenland ice sheet. Airborne radar campaigns by NASA’s Operation IceBridge Accumulation Radar (hereafter OIB AR) between 2010 and 2018 observed ice slabs covering at least 60 400–73 500 
$\mathrm{km^2}$ (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). In comparison, satellite-based passive L-band radar data suggested that they covered 76 000 
$\mathrm{km^2}$ between 2015 and 2019 (Miller and others, Reference Miller, Culberg, Long, Shuman, Schroeder and Brodzik2022), while retrievals made using Sentinel-1’s C-band radar increased their ice-sheet-wide coverage further to as much as 148 800 
$\mathrm{km^2}$ in 2016–17 (Culberg and others, Reference Culberg, Michaelides and Miller2024). Finally, modelling suggests that slabs could be even more widespread, covering 230 000 
$\mathrm{km^2}$ (Brils and others, Reference Brils2024).
According to airborne radar data, ice slabs had already started to develop in the early 2000s (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023; Brils and others, Reference Brils2024). Between 2010 and 2018, OIB AR flights along repeat transects showed ice slabs both expanding to higher elevations and thickening, especially during the extreme melt summers of 2010 and 2012 (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). At the uppermost limits of slabs, firn replenishment has been observed since 2012, burying nascent slabs and likely making them unavailable to support meltwater runoff (Culberg and others, Reference Culberg, Schroeder and Chu2021; Rennermalm and others, Reference Rennermalm2021; Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023; Rutishauser and others, Reference Rutishauser2024).Meanwhile, slab thickening at lower elevations did not occur uniformly, instead producing slabs of heterogeneous thickness over scales of hundreds of metres as ice accreted both on top of the slabs and below them (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023).
It has been proposed that under strong melt conditions the ice slabs may limit the maximum elevation of the runoff limit (Mikkelsen and others, Reference Mikkelsen2016; Machguth and others, Reference Machguth, Tedstone and Mattea2022; Tedstone and Machguth, Reference Tedstone and Machguth2022), yet the minimum slab thickness which can sustain runoff remains unknown. Furthermore, current state-of-the-art regional climate models (RCMs) used to calculate the surface mass balance of the ice sheet overestimate the runoff area (Tedstone and Machguth, Reference Tedstone and Machguth2022). Reducing this uncertainty would help to better constrain the ice sheet’s surface mass balance. Here, we first investigate the minimum average ice slab thickness which enables visible surface runoff, by comparing OIB AR ice thickness data to visible runoff limit retrievals from Landsat during the extraordinary melt season of 2012 (Tedesco and others, Reference Tedesco2011; Reference Tedesco2013; van Angelen and others, Reference van Angelen, van den Broeke, Wouters and Lenaerts2014). We use 2012 because the melt was so strong that most if not all of the snow and firn on top of the ice slabs was removed; on the K-Transect in south-west Greenland the ice slab was exposed to the surface (Machguth and others, Reference Machguth2016). Second, we analyse measurements of ice slabs thickness in 2017–18 alongside visible runoff limits during the large melt season of 2019 (Tedesco and Fettweis, Reference Tedesco and Fettweis2020) to consider heterogeneity in slab thickness. We examine the role that the surface hydrological network plays in producing this heterogeneity by transporting and coalescing water into localised areas which presumably subsequently refreeze.
2. Data and methods
2.1. Visible runoff limits
We used annual visible runoff limits extracted from Landsat imagery with a posting of 1 km parallel to the ice sheet margin (Tedstone and Machguth, Reference Tedstone and Machguth2022). In brief, these visible runoff limits were derived from the channel head locations of the highest-elevation streams and rivers which developed on the ice sheet surface each year. They are available ice-sheet-wide, except (i) areas underlain by firn aquifers and (ii) areas with particularly complex topography, namely the eastern margin south of 
$75^{\circ}\,\mathrm{N}$ and the western margin south of 
$\sim$$63^{\circ}\,\mathrm{N}$ and in between 75.9 and 
$76.7^{\circ}\,\mathrm{N}$ (Fig. 1c). For more details, we refer the reader to Tedstone and Machguth Reference Tedstone and Machguth2022.
Figure 1. Definition of zones around the visible runoff limit. (a) Normalized Difference Water Index (NDWI) map during 2012 (background), used to manually filter the maximum visible runoff limit retrievals (black and orange circles). The runoff limit zones are shown by shaded areas. Accumulation radar flight-lines in 2011 and 2012 are displayed as thin black lines, overlaid with thick coloured lines where they intersect runoff limit zones. (b) Ice slab thickness distribution in the different runoff limit zones. (c) Greenland ice sheet map showing the location of panel a (bold black box) and drainage basins from Rignot and Mouginot Reference Rignot and Mouginot(2012) (black lines).
Here we use the visible runoff limits from the 2012 and 2019 melt seasons, which were the first and second highest runoff years on record respectively (Nghiem and others, Reference Nghiem2012; Tedesco and others, Reference Tedesco2013; Tedesco and Fettweis, Reference Tedesco and Fettweis2020). We manually filtered these data further to remove any retrievals which we qualitatively evaluated to be located too low compared to the surrounding retrievals, as follows. At higher elevations, we compared the retrievals to annual mosaics of the Normalized Difference Water Index (NDWI) derived from the 
$\mathrm{10\mathrm{th}}$ percentile of Landsat’s red and blue bands we computed in the Google Earth Engine environment (Gorelick and others, Reference Gorelick, Hancher, Dixon, Ilyushchenko, Thau and Moore2017) (Fig. 1a). We removed any remaining points that were either not associated with visible surface water (possibly relating to cloud artefacts) or likely not connected to the downstream visible surface hydrological network (e.g. isolated supraglacial lakes).
2.2. Definition of zones around the visible runoff limit
While the visible runoff limit identified from satellite imagery constitutes a boundary line (Fig. 1a), the actual runoff limit lies some distance upstream as sufficient meltwater needs to coalesce before eventually appearing at the surface (Clerx and others, Reference Clerx2022). It is, therefore, instructive to consider a zone several kilometres wide, reflecting the transition from 100% retention to appreciable (i.e. visible) runoff (Holmes, Reference Holmes1955; Greuell and Knap, Reference Greuell and Knap2000; Humphrey and others, Reference Humphrey, Harper and Pfeffer2012; Clerx and others, Reference Clerx2022). To capture the ice slab thickness that is characteristic above and below the 2012 and 2019 visible runoff limits, we defined four zones around each year’s visible runoff limit (Fig. 1a). We hypothesise that within these several kilometres-wide zones, a change in ice slab thickness occurs which governs the transition from a meltwater retention regime to a meltwater runoff regime. The zone ‘at’ the visible runoff limit extends 0.5 km either side of the visible runoff limit. Above the visible runoff limit, we considered two zones. First, summertime observations at the runoff limit of SW Greenland indicate that water can flow laterally atop an ice horizon in a snow/slush matrix for ∼4 km before becoming visible at the surface (Holmes, Reference Holmes1955; Clerx and others, Reference Clerx2022). Hence, we defined an ‘in-between’ zone as extending between 0.5 and 4.5 km above the visible runoff limit. Second, we defined an ‘upstream’ zone as extending from 4.5 to 9.5 km above the visible runoff limit. Finally, the zone ‘downstream’ extends between 0.5 and 5.5 km downslope of the visible runoff limit.
2.3. Ice slabs extent and thickness
We use the ‘maximum likely ice slabs’ dataset derived from OIB AR between 2010 and 2018, which was first presented in Jullien and others Reference Jullien, Tedstone, Machguth, Karlsson and Helm(2023). The AR operates at 550–900 MHz, has an average horizontal sampling resolution of 16 m and a vertical resolution of ∼0.65 m in snow and firn (Rodriguez-Morales and others, Reference Rodriguez-Morales2014). The original ice slabs dataset by Jullien and others Reference Jullien, Tedstone, Machguth, Karlsson and Helm(2023) consists of ice slab thickness retrievals from airborne campaigns between 2010 and 2018, in which an ice slab is at least 1 m thick and no thicker than 16 m. They derived the ‘maximum likely’ dataset by using the maximum radar signal strength threshold that was indicative of slabs according to comparison with ground-penetrating radar data acquired coincident with an OIB AR flight.
To compare the ice slab extents with 2012’s visible runoff limit, we used the existing 2010–12 ice slabs extent in Jullien and others Reference Jullien, Tedstone, Machguth, Karlsson and Helm(2023). To analyse the ice slab thickness in the ‘downstream’, ‘at’, ‘in-between’ and ‘upstream’ zones, we used the ice slab thickness retrievals from Jullien and others Reference Jullien, Tedstone, Machguth, Karlsson and Helm(2023), retaining both (i) ice thicker than 16 m (maximum thickness of 20 m) and (ii) thinner than 1 m, and (iii) we assigned 0 m ice thickness to those parts of the flight-lines where no ice was identified. An example of the ice slabs thickness extraction in the different zones is shown in Figure 1b. We used OIB AR flight-lines acquired in 2011–12 and 2017–18 in the vicinity of the 2012 and 2019 visible runoff limits, respectively. To assess the variance in ice slab thickness, we calculated the ratio between the median absolute deviation (MAD) and the median slab thickness in each zone. Conceptually, as the value of this ratio increases, the slab thickness becomes more variable.
2.4. Constraining the minimum ice thickness which supports visible runoff
To constrain the minimum ice thickness required to support visible runoff, the upslope expansion of the visible runoff limit must be limited by sub-surface properties rather than by melt production (Machguth and others, Reference Machguth2016; Reference Machguth, Tedstone and Mattea2022). This could have been the case in both 2012 and 2019 as summer melting was exceptional in these years (Tedesco and others, Reference Tedesco2013; Tedesco and Fettweis, Reference Tedesco and Fettweis2020). However, (i) the visible runoff limit in 2019 was located slightly lower than in 2012 (Tedstone and Machguth, Reference Tedstone and Machguth2022), (ii) accretion of ice at the top of ice slab in the intervening years meant that the slab beneath 2019’s visible runoff limit was thicker than in 2012 (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023) and (iii) partial firn replenishment on top of ice slabs took place following summer 2012 (Rutishauser and others, Reference Rutishauser2024), burying some of the highest slabs at depth (Rennermalm and others, Reference Rennermalm2021). In contrast, in summer 2012, the progression of the visible runoff limit was very likely stopped by the structure of the underlying firn rather than by insufficient melt or abundant pore space, at least on the K-Transect (Machguth and others, Reference Machguth, Tedstone and Mattea2022).
We further consider the suitability of the 2012 season for this analysis by examining the annual melt-over-accumulation (MoA) ratio. Theoretical considerations of meltwater refreezing in firn show that when the MoA ratio is larger than ∼0.7 then runoff will occur (Pfeffer and others, Reference Pfeffer, Meier and Illangasekare1991; Brils and others, Reference Brils2024). MacFerrin and others Reference MacFerrin(2019) found that ice slabs have formed where this threshold was exceeded for a decade or more, although in situ observations on the K-Transect in south-west Greenland showed that meltwater runoff can commence at a lower MoA ratio of 0.6 (Braithwaite and others, Reference Braithwaite, Laternser and Pfeffer1994). Here, we examine the MoA ratio during 2001–11 and 2012, computed using melt and accumulation outputs from the regional climate model MAR v3.14 (Fettweis and others, Reference Fettweis2020).
2.5. Idealised ice slab thickness
The thickness of ice slabs can vary by several metres over horizontal distances of only a few hundred metres (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). To investigate the mechanisms that may explain the variability in ice slab thickness, we compared observed ice slab thickness with an idealised ice slab transect whose thickness is assumed to be solely determined by the local vertical percolation and refreezing of meltwater (i.e. a function of local melt intensity only).
On our idealised slab transect, the ice thickness decreases linearly between the lowest and highest elevations, modulated by a noise term. It is based on the characteristics of a radargram acquired during spring 2018 in SW Greenland (transect D in Jullien and others Reference Jullien, Tedstone, Machguth, Karlsson and Helm(2023); also shown here in Fig. 5f). The ice thickness along the transect was computed as follows. First, we applied a rolling average over a 3 km window to remove noise from the radargram. We tested windows of 1–4 km wide, at intervals of 1 km, choosing 3 km as this was found optimal to smoothen the noise in the radargram yet preserve the larger-scale patterns of ice slab thickness variability. Second, we set the ice thickness minima and maxima of the idealised transect to the same values as the thinnest and thickest average ice thickness of the observed ice slab (i.e. Fig. 5f) after smoothing. Next, we applied random noise to mimic the noise in the observed radargram. The range in the noise that we apply is the maximum of the absolute first discrete difference of ice thickness determined from a 5 km-long section of the radargram which exhibits uniform ice thickness. Finally, we assess the spatial variability in ice thickness by calculating the local coefficient of variation (standard deviation divided by mean at each sampling point) after applying the 3 km moving average, which is the parametric assessment of the MAD/median ratio described previously.
2.6. Frequency of surface hydrology development
To assess the relationship between the thickness of ice slabs and the surface hydrological network, we use a Landsat-derived Greenland ice sheet-wide 1985–2020 hydrological reference map by Tedstone and Machguth Reference Tedstone and Machguth2022. Briefly, it consists of the count of visible surface hydrology features which occurred between 1985 and 2020, at 30 m resolution. The reference map was produced by applying a series of filtering operations to each constituent individual satellite scene to prioritise the detection of linear channel-like features extending multiple kilometres over the ice sheet surface. Landsat 7 Enhanced Thematic Mapper (ETM+) Scan Line Corrector (SLC)-off data were gap-filled by interpolation. Several cloud-free acquisitions were generally made at any given location each summer, especially from 1999 onwards following the launch of Landsat-7 (Extended Data fig. 4 in Tedstone and Machguth Reference Tedstone and Machguth2022). The filtered satellite scenes were then stacked together to form the reference map. We divided the map by the ice sheet’s drainage basins (Rignot and Mouginot, Reference Rignot and Mouginot2012) and then normalised the map in each basin following a minimum–maximum normalisation approach to scale the frequency of surface hydrology between 0 and 1. The larger the frequency of surface hydrology, the more frequently a hydrological feature developed at a given location between 1985 and 2020.
3. Results
3.1. Ice slab thickness around the visible runoff limit
In 2012, the visible runoff limits were strongly co-located with the upper limits of ice slabs (Fig. 2), with a median separation of 0.6 km (Supplementary Table 1). In a few locations, the 2012 visible runoff limit extended slightly above the mapped ice slab area, but never beyond 63 m in elevation or 11 km in distance (maximum in CW Greenland, Supplementary Table 1). These discrepancies may be related to insufficient flight-line coverage leading to an inaccurate or incomplete mapping of ice slabs. MoA > 0.7 extended well beyond the visible runoff limit (Fig. 2), providing further evidence that 2012’s maximum visible runoff limit was not limited by meltwater availability; it is also clear that the majority of slab areas experienced modelled MoA > 0.7 throughout the decade preceding 2012. Thus, there are multiple lines of evidence which support the finding that 2012’s maximum visible runoff limits were controlled by the presence of ice slabs.
Figure 2. Maximum visible runoff limits in 2012 were controlled by ice slab extent. (a–d) The number of years between 2001 and 2011 when the MoA ratio computed from MAR model outputs exceeded 0.7 (background), the MoA = 0.7 contour in 2012 (green line), the 2010–12 ice slabs extent (white) and the 2012 visible runoff limit (light blue). (e) Elevation differences between the 2010–12 ice slab maximum elevation and the 2012 runoff limit in each region. A negative difference indicates that the ice slab was located at a lower elevation compared to the runoff limit. (f) Context map showing the location of panels a–d, and the MoA = 0.7 contour in 2012 (green lines).
We next analyse ice slab thickness which supported surface meltwater in the different zones around the visible runoff limit (Fig. 3). In 2012, the ice slab thickness in each basin was significantly larger ‘downstream’ of the visible runoff limit compared to ‘upstream’ (Welsch’s t test, p < 0.01). Overall, the ‘downstream’ ice slabs had a median thickness of 4.7 m. The distribution of the ice slabs thickness sampling shows that 75% of OIB retrievals in this zone were thicker than 2.6 m (Fig. 3). ‘At’ the visible runoff limit, the ice slabs remained thick (median 3.5 m). However, the ice thickness was more variable ‘at’ the visible runoff limit (MAD/median = 0.6) compared to ‘downstream’ (MAD/median = 0.5) (Supplementary Table 2). In the ‘in-between’ zone, ice slabs were already present as highlighted by a median ice thickness of 2.1 m, but with substantially more variability (MAD/median = 1.0). Finally, ‘upstream’ of the visible runoff limit was mostly free from ice slabs (median ice thickness of 0 m, Supplementary Table 2). The variability in slab thickness ‘upstream’ of the visible runoff limit cannot be assessed as the median and MAD were both 0 m apart from in the NW (Supplementary Table 2).
Figure 3. Ice slab thickness measured by OIB AR in 2011–12 in each zone around the 2012 visible runoff limit. Letter-values plot follows Hofmann and others Reference Hofmann, Wickham and Kafadar(2017): the tallest box contains the median (black vertical glyph) and spans the 25th–75th percentiles. The next box on both sides cover together with the first box 75% of the distribution, the following ones 87.5% and so on.
Given that meltwater is present on the ice sheet surface up to the visible runoff limits employed in this study, we can first define the conservative minimum thickness associated with visible surface hydrology over at least several kilometres as the median ice thickness ‘at’ the visible runoff limit, as on average 3.5 m. We note that the thinnest slab retrievals ‘at’ the visible runoff limit, i.e. the 
$\mathrm{0\mathrm{th}}$ percentile, are not necessarily indicative of the minimum ice thickness that enables visible surface runoff, primarily because the algorithm to extract slabs from radar data has an overall performance of 88%; in practise, this means that it mis-identifies some slabs as porous firn (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). Secondarily, despite the 1 km posting of the visible runoff limit, it is possible that some ice thickness measurements were acquired over hydrologic catchment divides between these points or where there were gaps in the retrieved visible runoff limit.
The differences in slab thickness between the ‘upstream’ and the ‘downstream’ zones give further insight into the minimum slab thickness which can support visible runoff, on the basis that the shift from probable complete retention to the onset of visible runoff occurs within these bounds. In all regions except the NW, the 
$\mathrm{25{th}}$ percentile of ice slab thicknesses in the ‘downstream’ zone is thicker than or similar to the 
$\mathrm{75{th}}$ percentile in the ‘upstream’ zone (Fig. 3). Furthermore, the median ice thickness ‘at’ the visible runoff limit is thicker than or similar to the 
$\mathrm{25{th}}$ percentile in the ‘downstream’ zone in all regions. We, therefore, define the smallest minimum thickness associated with visible surface hydrology over at least several kilometres as 2.8 m, i.e. corresponding to the 
$\mathrm{25{th}}$ percentile of the ‘downstream’ zone averaged over all regions.
‘Upstream’ from the visible runoff limit, 2–5% (0% in the NE, 42% in the NW) of the ice slabs sampled by OIB AR are thicker than the conservative minimum thickness, and 3–10% (0% in the NE, 49% in the NW) are thicker than the smallest minimum thickness. The ice slabs thickness distribution in the ‘in-between’ zone is sometimes closer to the distribution of those in the ‘upstream’ zone (in the SW, NO and NE) and other times closer to the distribution of those in the ‘at’ zone (CW and NW).
3.2. Expected ice slab thickness variability as a function of local melt
By considering an idealised ice slab transect, we can examine if an ice slab’s thickness is determined only by local melt (whose quantity we take to decrease with increasing elevation). Our idealised slab has maximum and minimum ice thicknesses of 16.5 and 4.6 m, with variability in thickness given by the introduction of random noise ranging between +3.0 and −3.0 m (Fig. 4a; Methods). The idealised transect has a median local standard deviation of 2.0 m, which is 0.1 m larger than the observed transect from which the idealised slab transect was derived. This indicates that the noise introduced in the idealised transect represents an upper bound for ice slab thickness heterogeneity.
Figure 4. Idealised ice slab transect. (a) Vertical ice slab extent along the transect relative to the surface. (b) Corresponding local ice thickness and local coefficient of variation.
On the idealised transect, the local standard deviation stays constant while the coefficient of variation gradually increases towards the highest elevations because the average ice thickness decreases (Fig. 4b). We compared our idealised ice slab transect with the 2018 radargram from which the idealised transect is derived (shown in Fig. 5f). Namely, we compared the coefficient of variation of the observed ice slab (black line in Fig. 5a) with the coefficient of variation that is expected from the idealised transect: for a given observed ice slab thickness, we extracted the corresponding coefficient of variation from the idealised transect (orange line in Fig. 5a). Within the section of the radargram where we extracted the noise signal for the idealised transect (between 6 and 11 km), the coefficient of variation extracted from the idealised transect is similar to the observations. In contrast, there are kilometres-long segments where the observed ice slab’s coefficient of variation is larger than the idealised value, indicating that the local melt intensity cannot be the only driver of ice thickness variability (Fig. 5a: 5–7 km, 14–21 km, 25–26.5 km).
Figure 5. Ice slab thickness, surface hydrology, topography and strain rates along a transect in SW Greenland. (a) Ice slab thickness (blue), actual coefficient of variation (black) and coefficient of variation derived from the idealised transect (orange) along the 2018 transect (shown in panel f). (b) Frequency of surface hydrology during 1985–2020 (background), and maximum visible runoff limits corresponding to panels c–f. (c–f) 0–20 m depth radargrams showing ice slabs in the firn from 2003 to 2018 (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). Runoff limit locations are displayed at the top of the radargrams. (c) The ice slab in the 2003 radargram is located in between sharp bright transitions where the signal strength is locally smoother. The 2005 runoff limit was the highest before 2010. (d–f) The ice slab in the radargrams from 2010 onwards is identifiable in darker grey, while porous firn appears in lighter grey. (g) Surface elevation (background) and 5 m contours (black lines). (h) Principal strain rates from Poinar and Andrews Reference Poinar and Andrews(2021). (g–h) The 2018 ice slab thickness along the OIB AR transect is displayed in shades of grey. Vertical dashed lines delineate sectors i–iv.
3.3. Impact of surface hydrology on slab thickness
In light of our proof that variability in ice slab thickness cannot be ascribed only to local melt, we investigated the degree to which surface hydrology can explain the variability that we observe. Figure 5a and b show that locations with a larger ice slab thickness coefficient of variation than expected often experience visible surface runoff (except in sector ii). Elsewhere, the ice slab thickness is equally or more homogeneous than could be expected from the idealised transect.
Radargrams acquired from 2003 to 2018 along this transect enable us to relate changes in the firn structure through time with visible surface hydrology (Fig. 5b–f). Note that the radargram in 2003 appears differently to subsequent years because of different radar setup (see Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). First, the visible runoff limits in 2005, 2010, 2012 and 2019 were always located a few kilometres downstream of the upper limits of visually contiguous ice slabs (Fig. 5c–f). Second, some of the largest increases in slab thickness coincide with areas where visible surface hydrology frequently developed. At 5 and 13 km along the transect, rivers formed frequently on top of thick ice slabs which were present as early as 2003. In sector iii and the higher elevations of sector iv, the ice thickened more uniformly through time, associated with slush fields of larger spatial extent than the rivers (Fig. 5c–f). In sector iv and occasionally upstream, while there was likely no ice slab in 2003, by 2018 there was a widespread ice slab, with locally thicker areas corresponding with the locations of ephemeral slush fields and rivers (Fig. 5d,f).
Away from areas of visible surface hydrology, the ice slab generally thickened by relatively homogeneous ‘top down accretion’ (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023) (except in sector ii, Fig. 5b,d–f). However, in sector ii, while the firn was mostly ice-free in 2010, by 2018 a thick ice slab had developed despite an absence of visible surface hydrology (Fig. 5b). ArcticDEM v3 (Porter and others, Reference Porter2018) shows that this sector is located on a slope downstream of a local plateau (Fig. 5g). The runoff limit from 2010 onwards was often located upstream of this sector, indicating that large amounts of meltwater were available at this elevation (Fig. 5d–f). Radargrams (Fig. 5d–f) show deep vertical structures which we interpret as the result of accretion of ice due to meltwater refreezing at depth, a phenomenon previously documented by Culberg and others Reference Culberg, Chu and Schroeder(2022). Positive winter-time principal strain rates (Poinar and Andrews, Reference Poinar and Andrews2021) indicative of extensional ice flow dominate the sector (Fig. 5h). This suggests that crevasses are likely opening roughly perpendicular to the west-to-east ice flow (Doyle and others, Reference Doyle2014). Thus, meltwater generated upstream of sector ii likely flows over existing sub-surface ice layers until it reaches crevasses which permit percolation deep into the firn column.
To examine the applicability of these processes across the rest of the ice sheet, we calculated the ice slab thickness coefficient of variation for all OIB AR measurements acquired during 2017–18 in the 5 km-wide zone ‘downstream’ of the 2019 visible runoff limit, using the same 3 km moving window described previously (Supplementary Fig. 1 and Supplementary Table 3). We use the 2019 visible runoff limit because the ice slabs thickness dataset from 2017–18 contains 39% more measurements in the ‘downstream’ zone than the dataset acquired during 2011–12 (Supplementary Tables 2 and 3). This most recent dataset also better corresponds to the time period over which the frequency of surface hydrology map was computed (until 2020) and captures several summers with strong melt (particularly 2012 and 2019, Nghiem and others, Reference Nghiem2012; Tedesco and Fettweis, Reference Tedesco and Fettweis2020). We discarded ice slab thicknesses lower than 1 m (corresponding to 4.3% of the dataset) which would otherwise artificially increase the coefficient of variation.
The inter-quartile range of the thickness observations shows more variability than the idealised coefficients of variation would indicate (0.17–0.40 versus 0.16–0.28). We investigate this variability by examining the local ice slab thickness anomaly, calculated by subtracting thickness at each OIB AR sampling point from the median thickness along its parent transect within the ‘downstream’ zone. Figure 6 shows that this thickness variability is driven at least partly by surface hydrological features. In all regions except the NE, the ice thickness anomaly as expressed by the inter-quartile range is more positive at locations with the higher frequencies of surface hydrology. In contrast, the NE has a reduced inter-quartile range at higher frequencies but is nonetheless skewed to more positive thickness anomalies as shown by the tails. We discuss the possible reasons for this contrast in Section 4.
Figure 6. Boxenplots of ice slab thickness anomaly separated by frequency classes of surface hydrology presence, in the ‘downstream’ zone. The bounds used for scaling the frequency of surface hydrology following min–max normalisation in each region are quantile 0.01 and quantile 0.99. N refers to the count in each class. Outlier points have been omitted for clarity.
4. Discussion
4.1. Constraining the minimum ice thickness which supports runoff
Our thickness threshold is extracted from ice slab thickness retrievals identified in radargrams acquired by OIB’s AR. We discuss here potential sources of uncertainties stemming from the use of these data. First, we note that the AR has a vertical resolution of ∼50–60 cm in firn and ice (Rodriguez-Morales and others, Reference Rodriguez-Morales2014; MacGregor and others, Reference MacGregor2021), which makes the ice slabs thickness retrievals accurate to at best half a metre. Jullien and others Reference Jullien, Tedstone, Machguth, Karlsson and Helm(2023) identified ice slabs in OIB AR data using radar strength thresholds derived from comparison between coincident OIB AR and ground-based radar transects. While the overall accuracy of the identification was high (88%) on the reference OIB AR radargram, this approach yields somewhat higher uncertainty in the identification of ice slabs in other OIB AR radargrams. Jullien and others Reference Jullien, Tedstone, Machguth, Karlsson and Helm(2023), therefore, performed a visual expert comparison between individual radargrams and their corresponding ice slab products, to exclude erroneous slab thickness retrievals. Here we use their maximum likely ice slab thickness retrieval (see Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023), which maximises the retrieval of ice slabs. While we cannot exclude that some ice slab thicknesses might be overestimated, we expect it to be of second order importance.
For the extreme melt season of 2012, we showed the ice slabs ‘downstream’ from the runoff limit are significantly thicker than they are ‘upstream’. A widespread MoA ratio larger than 0.7 extending far beyond the visible runoff limit (Fig. 2) suggests that ice in the firn column was insufficient to support surface meltwater runoff from higher elevations despite abundant melt supply. Our results, therefore, suggest that at first order the presence of ice slabs governs whether surface meltwater runoff can occur once MoA exceeds ∼0.7. The thickness of the ice slabs controls the maximum possible elevation of the visible runoff limit; we find that the minimum ice thickness that enables visible surface meltwater lies between 2.8 m (smallest minimum thickness) and 3.5 m (conservative minimum thickness). In the different regions ‘downstream’ from the runoff limit, between 60 and 73% of the ice slabs are thicker than the minimum conservative thickness (20% in the NE), and 69–82% are thicker than the smallest minimum thickness (34% in the NE). As the distribution of ice slabs thickness indicates that slabs are effectively contiguous ‘downstream’ from the runoff limit, we, therefore, argue that the thresholds we derived correspond to the ice slab thickness range that allows them to support widespread visible surface meltwater over scales of at least several kilometres.
‘Upstream’ from the runoff limit, less than 10% of the ice slabs sampled by OIB AR are thicker than both the smallest and conservative minimum thicknesses (all regions except the NW). This suggests that in general, very little of the upstream area was underlain by slabs thick enough to support visible meltwater runoff in summer 2012. Although the slab thickness is heterogeneous in the ‘in-between’ zone, on an ice-slabs-wide basis more than 25% of the observations in this zone are already thicker than the smallest minimum thickness (Fig. 3). We, therefore, suggest that this ‘in-between’ zone is already primed to support visible runoff under projected warming.
We can contextualise our minimum thicknesses using in situ measurements at KAN_U in south-west Greenland. Here, in spring 2012, firn core measurements showed the presence of three ice layers of at least 1–1.5 m thick each within a 6 m section (Machguth and others, Reference Machguth2016) while contemporaneous OIB measurements indicated a 2–3 m thick ice slab (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). The uppermost extent of the ice slab was roughly located at KAN_ U during this time (Machguth and others, Reference Machguth2016; MacFerrin and others, Reference MacFerrin2019). There is clear evidence that in summer 2012 KAN_ U was located at the visible runoff limit (Machguth and others, Reference Machguth2016; Mikkelsen and others, Reference Mikkelsen2016). The maximum elevation of visible meltwater plateaued at KAN_ U during summer 2012, despite substantial melting that could have enabled runoff from ∼300 m higher elevations if a near-surface aquitard had been present to support it (Machguth and others, Reference Machguth, Tedstone and Mattea2022). Thus, the ice slab thickness and surface hydrology at KAN_ U during summer 2012 support our smallest and conservative minimum estimates of 2.8–3.5 m thick as a reasonable range for the minimum contiguous ice thickness which can support widespread runoff.
We found significant differences in the ice slab thickness between the ‘upstream’ and ‘downstream’ zones (Fig. 3), which suggests that the transition from complete retention to visible runoff (presumed to be indicative of at least partial runoff) occurs over distances of less than 15 km. Firn cores acquired in spring 2013 (Machguth and others, Reference Machguth2016) and spring 2021 (Clerx and others, Reference Clerx2022) along the K-transect show that the structure of the shallow firn differs widely between the runoff limit and several kilometres further inland. Both these studies showed that areas which had recently experienced visible runoff were underlain by ice slabs >5 m thick (e.g. Core_1_2013 in Machguth and others Reference Machguth(2016); Core FS2 in Clerx and others Reference Clerx(2022)). Meanwhile, porous firn was present ∼17 km and ∼18 km further inland (110–130 m higher in elevation), which was subjected to the vertical percolation and refreezing of meltwater (Core_3_2013 in Machguth and others Reference Machguth(2016); Core FS4 in Clerx and others Reference Clerx(2022)). Finally, based on field observations at the runoff limit on the western ice sheet margin at around 
$70^{\circ}\,\mathrm{N}$, Humphrey and others Reference Humphrey, Harper and Pfeffer(2012) suggested that the transition from complete runoff to complete retention occurs over a ∼20 km wide zone. Our conclusions are in broad agreement with their findings.
4.2. Ice slabs, surface hydrology and topography
The transect analysis in Figure 5 highlights that local melt intensity alone cannot explain ice thickness spatial variability. Ice slabs with locally higher thickness are associated with more frequent occurrence of visible surface hydrology (Fig. 6; see also Supplementary Fig. 2). This finding supports the refreezing of meltwater away from where it was produced, thus indicating that lateral transport does not necessarily lead to runoff. Refreezing downstream of meltwater production locations can occur by ‘top down accretion’ (Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023) or at depth (Culberg and others, Reference Culberg, Chu and Schroeder2022; Jullien and others, Reference Jullien, Tedstone, Machguth, Karlsson and Helm2023). Both mechanisms locally increase the slab thickness (Fig. 5 and Supplementary Fig. 2). However, as observed in Figure 5, ‘at-depth’ accretion can occur at locations distant to visible surface hydrology with extensive strain rates. This may explain why the differences in thickness anomalies between less versus more frequent visible hydrology are relatively small: meltwater may transit via the sub-surface out of sight before refreezing (Clerx and others, Reference Clerx2022), generating local thickness anomalies in locations without visible surface hydrology.
A further mechanism which may be responsible for the wide spread in thickness anomaly in areas with infrequent visible hydrology is as follows. While the surface hydrological network often re-establishes in the same locations each year (e.g. Holmes, Reference Holmes1955), in contrast ice slabs are multi-annual features which get advected downslope through time (Jullien and others Reference Jullien, Tedstone, Machguth, Karlsson and Helm(2023)). In some cases, we can, therefore, find thicker slabs located down-glacier of the most frequent surface hydrological features, especially in locations where ice flow is faster.
Topography is also likely to play a role in slab thickness. At the runoff limit, ice slabs are generally thinner in the CW compared to the other regions. The steeper surface slopes in the CW might route meltwater more efficiently downstream, such as suggested by Gascon and others Reference Gascon, Sharp, Burgess, Bezeau and Bush(2013) on the Devon Ice Cap. In the NO, ice slabs are thicker in the four different zones around the runoff limit compared to other regions, while in the NW only slabs in the ‘upstream’ zone are thicker than other regions (Fig. 3). This might be explained by the fact that the 2012 visible runoff limit in most of the NW and the NO was only the second or fifth highest visible runoff limit respectively since 1985, unlike in the CW and the SW where it was the highest (Tedstone and Machguth, Reference Tedstone and Machguth2022). Finally, the distribution of ice slabs thickness in the CW and the NW deviates slightly compared to the other regions. While the visible runoff limit position in the SW and NO is relatively steady (Fig. 2) as a function of shallow surface gradients and simple topography, the visible runoff limit in the CW and NW exhibits a somewhat ‘saw-tooth’-like pattern which likely originates in steeper gradients and more complex topography. Consequently, the zones around the visible runoff limit capture more variability in surface topography, which could introduce uncertainty in the sampled ice slabs thickness. Finally, comparing lakes mapped by Dunmire and others Reference Dunmire, Banwell, Wever, Lenaerts and Datta(2021) with ice slab extents shows that numerous surface and subsurface lakes are observed in the slabbed areas of the NW and CW, while few are present in the SW and the NO. Culberg and others Reference Culberg, Chu and Schroeder(2022) showed that supraglacial streams and lakes in the NW can drain through crevasses into relic firn beneath ice slabs. Therefore, further work is required to understand relationships between ice slab thickness, meltwater routing and surface topography at process scales.
4.3. Thin ice slabs in the north-east
Around the visible runoff limit, ice slabs in the NE are substantially thinner than in other regions (Fig. 3, Supplementary Fig. 1). Our results suggest that in the NE, 1 m thick ice slabs can already support surface meltwater runoff, which contradicts our findings from elsewhere on the ice sheet. The statistics in the NE are likely hampered by uncertainties in the location of 2012’s visible runoff limit. There is a strong gradient in slab thickness around 2012’s visible runoff limit, with thicknesses exceeding three metres 5 km downstream, and thicknesses which are close to the ice slab minimum thickness (1 m) within several kilometres upstream (Supplementary Fig. 2).
The uncertainties in identifying 2012’s visible runoff limit in the region NE appear to be related to the north-east Greenland ice stream. The multi-annual surface hydrological network has ‘smeared’ structures in areas of high ice velocities (Supplementary Fig. 2). Slightly north of the ice stream, where ice velocities are more similar to values typical for Greenland’s ice slab areas, then an arborescent stream network similar to other regions develops. Given the diffuse pattern of surface hydrology in the NE which resulted in a noisy automatic retrieval of the visible runoff limit, we had to manually remove several runoff limit outliers in this study by reference to satellite imagery (Methods), yet its exact position remains open to interpretation and is likely the main reason for the identification of the very thin ice slabs which appear to support runoff.
Regardless of uncertainties in the 2012 runoff limit, Supplementary Figure 2 clearly shows that the surface hydrology and thickness of ice slabs are related. With the exception of one OIB profile, we found a good spatial agreement between the upper boundary of ice slabs and 2012’s visible runoff limit, and also that surface streams as well as lakes are in most cases underlain by thicker slabs than in their direct vicinity.
4.4. Implications for modelling meltwater runoff
The surface mass balance of the Greenland ice sheet is commonly estimated using RCMs which incorporate a coupled 1-D (vertical) surface energy balance model to determine exchanges between the atmosphere, ice sheet surface and sub-surface, and are run at spatial resolutions of several kilometres (e.g. Fettweis and others, Reference Fettweis2020). Tedstone and Machguth Reference Tedstone and Machguth2022 showed that two models, MAR and RACMO, overestimate the runoff area, which suggests that parameterisations of physical processes in the vicinity of the runoff limit require improvement.
We propose that an improved RCM parameterisation of the transition from vertical percolation to lateral runoff could be centred at first order around the vertical firn column becoming impermeable to further percolation as the modelled slab thickens from 2.8 to 3.5 m. Beyond 3.5 m thick then meltwater in excess of the storage capacity of the overlying pore space can be considered to run off laterally instead. We consider that such a parameterisation would be valid for horizontal model resolutions of several kilometres, whereby the minimum slab thickness required to support lateral runoff implicitly captures the necessary slab contiguity across the grid cell. We stress that this does not mean that slabs only become impermeable beyond 2.8 m thick. Their 1-D impermeability threshold is likely much smaller (e.g. Clerx and others, Reference Clerx2022) but, according to our findings, is unlikely to be contiguous enough to support widespread visible runoff (in other words, the water will flow laterally over scales of only several metres until it encounters accessible pore space).
We caution, however, that a runoff criterion based on slab thickness alone is likely too simplistic. Our findings show that ice slabs are thicker at hydrologically-active locations. This suggests that although meltwater enters a supraglacial hydrological network, at least some of this water may refreeze atop or beneath the ice slabs in these locations rather than run off. Focusing on refreezing at the top, recent quasi 2-D modelling of lateral meltwater flow in sluch fields suggests that it may flow at least several kilometres before it gets refrozen (Clerx and others, Reference Clerx, Machguth, Tedstone and van As2024), raising the question of whether such a process can be captured within existing 1-D surface schemes. Refreezing on to ice slabs but also beneath them at places, therefore, needs further attention if we are to have confidence in modelled estimates of runoff from higher elevations. Finally, we emphasise the importance of modelling high-elevation surface melt accurately. This remains at outstanding challenge as RCMs generally use the skin layer formulation in their surface energy balance scheme, which is liable to bias surface temperatures and thereby over-estimate melting Covi and others Reference Covi, Hock and Reijmer(2022).
5. Conclusions
We have shown that the presence of ice slabs controlled the maximum possible location of the visible runoff limit in 2012, enabling us to use 2012 to determine the minimum slab thickness which supports visible surface runoff according to ice thickness retrievals by airborne radar. The shift from a retention-dominated regime to a regime where meltwater can potentially run off over the ice sheet surface is associated with slabs 2.8–3.5 m thick (smallest and conservative minimum thicknesses, respectively). We infer that this range of minimum thicknesses is indicative of ice slabs that are spatially contiguous enough to support lateral meltwater runoff. Nonetheless, it is important to recognise that in most of the recent summers the maximum position of the visible runoff limit has been controlled by melt magnitude, not slab extent. Ice slabs only support runoff when melting in a given melt season is sufficient to exceed the refreezing capacity of snow and firn pore space located above those slabs, otherwise any meltwater will simply be refrozen in situ.
Once slabs have formed, subsequent local increases in ice slab thickness are largely associated with areas that undergo visible surface runoff, indicating that some meltwater refreezes downstream of its production area. This highlights that although ice slabs support runoff, they can also mediate total mass losses by supporting refreezing in overlying (and sometimes underlying) pore space. Their ability to support refreezing in the future will depend on whether annual melting routinely exceeds annual accumulation.
Over scales of several kilometres, we consider ice slabs whose thickness is larger than our contiguity threshold to be effectively impermeable to meltwater percolation, thereby supporting surface meltwater runoff. Implementing such a threshold in the firn evolution schemes used in RCMs can provide a first step to improve the parameterisation of the transition from firn retention-dominated regimes to firn regimes where some meltwater runs off. We would expect that such an approach will improve the partitioning between meltwater refreezing and runoff, a key aspect when estimating future mass loss rates from the Greenland ice sheet under projected warming (Intergovernmental Panel on Climate Change (IPCC), 2023).
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.25.
Data availability statement
Cleaned maximum visible runoff limit retrievals in 2012 and 2019, the modified 2010–18 ice slabs thickness dataset and the connectivity map from Tedstone and Machguth Reference Tedstone and Machguth2022 can be downloaded at: https://doi.org/10.5281/zenodo.14237814. 2002–03 radargrams, ice slabs products and ice slabs extent from Jullien and others Reference Jullien, Tedstone, Machguth, Karlsson and Helm(2023) are accessible at https://zenodo.org/records/7505426. Runoff limit retrievals and excess melt calculations from Tedstone and Machguth Reference Tedstone and Machguth2022 are accessible at https://zenodo.org/records/6472348. 10 m resolution mosaics from the ArcticDEMv3 (Porter and others, Reference Porter2018) are accessible at: https://data.pgc.umn.edu/elev/dem/setsm/ArcticDEM/mosaic/v3.0/. Winter time principal strain rates map from Poinar and Andrews Reference Poinar and Andrews(2021) are accessible at: https://ubir.buffalo.edu/xmlui/handle/10477/82127. Ice sheet drainage basins from Rignot and Mouginot Reference Rignot and Mouginot(2012) can be downloaded at: http://imbie.org/imbie-3/drainage-basins/. The scripts used to perform the analysis for this study can be found at https://github.com/jullienn/IceSlabs_SurfaceRunoff.
Acknowledgements
This work was funded under the European Research Council award 818994—CASSANDRA. AT was also funded by Swiss National Science Foundation grant TMSGI2_218095—FlowState. We acknowledge the Polar Geospatial Center from the University of Minnesota’s College of Science and Engineering for the use of ArticDEM. We acknowledge the use of data and data products from CReSIS generated with support from the University of Kansas, NASA Operation IceBridge Grant NNX16AH54G, NASA Grant NNX10AT68G, NSF Grants ACI-1443054, OPP-1739003, and IIS-1838230, ANT-0424589, Lilly Endowment Incorporated and the Indiana METACyt Initiative. We thank NB Karlsson for helpful discussions in the initial phase of this work.
Author contributions
Conceptualisation, methodology, software, resources: NJ, AT.
Formal analysis, visualisation: NJ, AT, HM.
Writing—original draft: NJ.
Writing—review and editing: AT, HM.
Supervision: AT, HM.
Funding acquisition: HM, AT.
Competing interests
The authors declare that they have no conflicts of interest.
Open access
