Concurrent superimposed ice formation and meltwater runoff on Greenland’s ice slabs

Domain: earth_science
Published: 2025-05-14
Editorial outcome: major_revision

Abstract
Rivers and slush fields on the Greenland Ice Sheet increasingly develop in locations where the accumulation zone hosts near-impermeable ice slabs. However, the division between runoff versus retention in these areas remains unmeasured. We present field measurements of superimposed ice formation onto slabs around the visible runoff limit. The quantity of superimposed ice varies by proximity to visible surface water and the surface slope, highlighting that meltwater can flow laterally before refreezing. We use heat conduction modelling and radar observations of autumn wetness to show that in our field area in 2022, 65% of superimposed ice formed during summer and the rest during autumn in the relict supraglacial hydrological network. Overall, 84% of melt around the visible runoff limit refroze. Ice-sheet-wide we estimate that slabs refroze 56 gigatonnes of melt (26-69 gigatonnes according to slab extent) between 2017 and 2022. Slabs are thus both hotspots of refreezing and emerging zones of runoff.

Introduction
In the accumulation area of the Greenland Ice Sheet where the surface is underlain by firn, meltwater generally percolates into the firn pore space where it refreezes in-situ, and therefore does not runoff. Surface mass balance (SMB) modelling indicates that in the past several decades 40-50% of ice-sheet-wide melt has refrozen in the firn 1–3. This makes firn an important buffer against future sea level rise 4,5. However, in the last three decades surface melting has increased more quickly than accumulation 6.

In the lower accumulation area the firn pore space has decreased anomalously as a result of intense meltwater percolation and refreezing 4,7, leading to metres-thick near-impermeable ice ‘slabs’ forming just below the ice sheet surface 8,9. Meanwhile, the visible runoff area—the area drained by surface river networks—has expanded by 29% in the past four decades, with concomitant increases in the elevation of the visible runoff limit 10. This expansion has mainly occurred along the ice sheet’s western and northern flanks in areas now underlain by ice slabs 11. According to synthetic aperture radar observations, ice slabs are estimated to underlie ~150,000 km2 of the ice sheet, though they may occupy as little as ~60,000 km2 or as much as ~220,000 km2 (ref. 12).

In-situ hydrological observations around the runoff limit on top of an ice slab show that surface meltwater percolates through the overlying snowpack and firn then flows laterally through a slush matrix on top of the ice slab before it reaches a surface stream or river. This matrix is up to several decimetres thick and supports flow speeds ranging from 1.3 to 15.1 m h-1 (ref. 13)—roughly three orders of magnitude lower than subsequent open river flow of 1–3 m s-1 (ref. 14). A slush matrix on top of a cold ice slab is therefore a potential source for superimposed ice formation onto the slab surface as its wet pore space gets refrozen from the bottom up 15.

On the Greenland Ice Sheet, studies of refreezing have generally focused on vertical infiltration and refreezing in porous firn above the runoff limit 7,16,17, sometimes considering decimetre-scale horizontal flow fingers 18. Field measurements north of Jakobshavn Isbrae during the 2000s were suggestive of slow water migration at depth beneath multiple ice layers but in the absence of ice slabs, indicative of a gradient between complete refreezing and complete runoff 5,19. Measurements of superimposed ice formation specifically have predominantly been made elsewhere such as Svalbard 20,21 or the Canadian Arctic e.g. 22,23. Thickening of Greenland’s ice slabs during years in which they supported surface river development can be observed in repeat airborne radar measurements 11, yet the contribution of refreezing in ice slab areas to SMB has not been quantified.

Here we present measurements of superimposed ice formation (hereafter SIF) on top of ice slabs on the K-Transect (67° N), south-west Greenland, acquired during two years in which the visible runoff limit intersected the area. We estimate the impact of SIF on net water discharge from the hydrological catchment and on SMB. Then, in conjunction with vertical heat flow modelling and synthetic aperture radar observations of sub-surface wetness, we make an assessment of the importance of SIF to SMB across all regions of the ice sheet underlain by ice slabs.

Results
We measured SIF in the 2021 and 2022 melt seasons through repeat in-situ observations at stakes located 1760–1825 m above sea level (a.s.l.) (Fig. 1a-c; Methods; Suppl. Fig. 1), more than 200 m higher than the long-term equilibrium line altitude of ~1550 m a.s.l 24. The average surface slope is ~0.3° 13. Across the elevation range spanned by our study area, melt derived from a surface energy balance model forced by in-situ observations 25 was relatively uniform, with 0.33–0.46 m w.e. in 2021 and 0.26–0.33 m w.e. in 2022 (Methods). In 2017 the area was underlain by an ice slab with an average thickness and standard deviation of 8.7 ± 3.9 m 11. The area developed surface slush fields and rivers or linear drainage pathways suggestive of runoff every summer between 2011 and 2023 except 2017 and 2018 10.

Between April/May 2021 and 2022 (hereafter ‘2021’) we measured SIF at three sites (FS1–3), and between April/May 2022 and 2023 (hereafter ‘2022’) we made measurements at eleven sites along two transects. One transect was perpendicular to the main visible drainage pathway through the area (and co-located with the line of summer hydrological measurements by ref. 13) while the second transect was normal to the first and ran west-east, parallel to the general surface slope.

SIF was spatially heterogeneous and corresponded with the locations of surface water presence near to the peak of the melt season and the surface slope, as follows (Fig. 1d). (i) Considering the largest measurements of SIF (BH0, BH1, FS2 in 2021), all were located within tens of metres of a major river or linear drainage pathway and on very shallow slopes. (ii) The other sites with relatively high SIF (FS1, BH5, 3u3) were further away from linear drainage pathways and did not necessarily have surface water nearby. At FS1 there was 0.21 and 0.21 m w.e. SIF in 2021 and 2022 respectively, even though a slush field was only visible nearby in 2021. Similarly, site 3u3 also recorded substantial SIF in 2022 despite being 4 km away from visible surface water, but Sentinel-2 imagery from recent high-melt years (e.g. 2019, 2023) shows that 3u3 lies within a few hundred metres of the axis of the same linear drainage pathway which flows past FS2. We therefore interpret the larger SIF at these sites as indicative of a higher local water table as the near-zero slope is conducive to ponding of water supplied by sub-surface flow inside the snowpack from the surrounding area 13. (iii) Considering the remaining measurements with less SIF (C1/C2, 122, FS3, 231, FS2 in 2022), all were located at least hundreds of metres away from visible surface water, i.e. more distal to the main drainage axis of the catchment, and all (apart from FS2) had steeper slopes compared to the other measurement sites. The role of a steeper slope in reducing SIF is particularly evident at site 231.

Considering only measurements from 2022 and excluding those sites proximal to a linear drainage pathway (BH0 and BH1), average SIF was 0.16 ± 0.05 m w.e. (± 1 standard deviation; n = 9). These sites had an average climatological surface mass balance (cSMB) calculated from repeat firn cores (Methods) of 0.53 ± 0.06 m w.e., of which SIF constituted 30 ± 9%. In comparison, BH0 and BH1 had a larger cSMB of 0.77 ± 0.09 m w.e. yr-1 (n = 2), of which SIF constituted 72 ± 3%. As melt is near-uniform across the area of our observations, these measurements reinforce that meltwater for SIF is supplied not only by local melt production but also by lateral water flow.

We now focus on the mechanisms of SIF onto the ice slab’s uppermost surface. We installed a 12 m thermistor string in a borehole at FS2 in May 2022, suspended from a mast frozen into the sub-surface (Methods). It recorded the warming of the overlying snowpack during the early melt season, becoming isothermal on 6 July which enabled meltwater delivery to the ice slab surface (Fig. 2a, b). We then observed subsequent ice slab warming to at least 6 m depth, indicating refreezing of the surface water onto the top of the slab surface which warmed the slab below as latent heat is released. We term this mode of SIF ‘bottom-up’ as the cold content comes from the underlying slab (Fig. 3, Summer).

To calculate how much SIF can occur during summer, we modelled SIF potential as a 0 °C slush matrix refreezing onto a cold impermeable ice slab due to heat conduction, a purely thermal approach that assumes infinite and instantaneous melt availability (Methods). We initialised our model with the vertical temperature profile measured on 6 July, the first day that the ice slab surface reached 0 °C. We modelled SIF and the slab temperature profile for the duration of the melt season, assuming a constant 0 °C at the slab surface.

Our calculations yield 0.15 m w.e. summer bottom-up SIF by the end of the 2022 melt season on 18 Sep (Fig. 2c), agreeing well with our measured total 0.16 ± 0.05 m w.e. yr-1 SIF away from surface hydrological features. Calculated temperature profiles show overall good agreement with measurements (Fig. 2d), with a slight warm bias: at the end of the melt season on 18 September, the modelled temperature profile was on average 0.4 °C warmer than observed. The observed 18 September profile was reproduced most closely by the model on 3 September, which suggests that the model overestimated SIF by 0.02 m w.e., perhaps due to continuing to simulate SIF during the melt hiatus in late August 2022. Nonetheless, overall our model reproduces the observed seasonal slab warming driven by latent heat released during SIF. These results are evidence that away from surface hydrological features, 94% of SIF (at 0.16 m w.e. yr-1; 68-100% considering the standard deviation of our observations of ± 0.05 m w.e. yr-1) can be explained by direct refreezing of local surface melt during summer once it has percolated through the overlying snow/slush matrix.

While summertime SIF was generally sufficient to accommodate local meltwater generation and percolation in 2022, our in-situ measurements also indicate that much larger quantities of SI formed in specific areas. We argue that these larger amounts of SI can be explained by lateral melt transport through the slush matrix (as observed by ref. 13) and localised ponding which provides residual meltwater supply for SIF to continue in specific locations in autumn, beyond the end of the melt season.

To map the duration of residual wetness in the uppermost decimetres to metres of the firn column after the melt season, we used satellite C-band Synthetic Aperture Radar (SAR) measurements of Horizontal-Vertical polarisation backscatter intensity, de-biased relative to winter time and then thresholded at -2.17 dB (Methods). This approach is similar to previous studies which have used Sentinel-1 to detect firn aquifers in Greenland 26 and snow melt in Antarctica 27. Five days after the end of melting, 60% of the study area remained wet; 30% was still wet after 18 days (Fig. 4a, b). By 66 days after melting, 99% of the area was refrozen.

In autumn freezing air temperatures predominate (Fig. 4c). These cool the snowpack from above, so ‘top-down’ SIF presumably occurs in parallel with bottom-up SIF as slab cold content continues to be consumed (Fig. 3). We modelled these two modes as follows: we calculated bottom-up SIF by continuing our summer model run into autumn, until residual wetness disappeared. Separately, we estimated top-down SIF by heat conduction away from the slush through the snowpack to the cold atmosphere (Methods). Based on our thermistor string measurements (Fig. 2), the snowpack thickness above the slush was approximately 0.6 m at the start of autumn. As this is the most sensitive parameter in our model, we assessed top-down SIF in our study area occurring beneath snowpacks of 0.6 m, 0.3 m (0.5×) and 0.9 m (1.5×) thickness. We then used our retrievals of residual wetness duration to identify the SIF contribution from each of the two autumn modes.

Autumn bottom-up SIF proceeds at a slower rate than in summer, reflecting the diminished ice slab cold content and with little sensitivity to the ± 1 °C uncertainty in the initial temperature at 10 m depth 28 (Fig. 4c). Meanwhile, the start of top-down SIF is delayed by the requirement to first refreeze irreducible water content in the overlying snowpack. It then commences at a similar rate to bottom-up SIF, subsequently diverging as a function of snowpack thickness, especially in the case of a thin snowpack which cools very quickly in comparison to a thicker pack. Whether a slush matrix that is thick enough to produce ~0.4 m w.e. top-down SIF over 60 d can co-exist with a sufficiently thin snowpack is unclear as we lack constraints from the field, but observations made in the area during summer 2020 13 suggest that such conditions are unlikely.

We summed our pixel-by-pixel estimates of summer and autumn SIF introduced above to derive total SIF over our field study area (Fig. 5a, b). Considering melt of 0.26–0.33 m between 1700 and 1800 m a.s.l., areas with surface hydrological features (e.g. around BH0) refroze substantially more meltwater than they produced. Our model-derived values of total SIF compare reasonably well with our in-situ measurements (Fig. 5c, R2 0.66, p < 0.01, n = 10). We note that the lack of modelled variability around 0.2 m w.e. is likely to be associated with the 12-day resolution of the Sentinel-1 observations. Considering the outlier of site 231, calculated SIF was roughly 0.3 m w.e. too high due to prolonged autumn wetness. This site was located in the margins of the wet subsurface area (Fig. 5b) but was not visibly wet and had a locally steeper slope (Fig. 1c, d); the wetness conditions observed by SAR over tens of metres therefore likely obscured the most local patterns of wetness.

We examined summed total SIF in the first 10 km of FS2’s upstream hydrological catchment, roughly corresponding to the plausible limit of sub-surface hydrological routing based on SAR residual wetness (Fig. 4a). In addition to our best-estimate scenario (Fig. 5a), we examined model sensitivity using 1.0 °C warmer (colder) initial 10 m depth temperatures based on ref. 28 firn temperature uncertainties for our study area, and 1.5× (0.5×) snowpacks, respectively (Fig. 5d). Across these scenarios, summer SIF is always the largest constituent of total SIF, while autumn bottom-up SIF is the next-largest component. Although the range in autumn top-down SIF associated with snowpack thickness is very large, it is a relatively small component of total SIF apart from in the thinnest snowpack scenario. The uncertainty in autumn SIF introduced by the radiometric sensitivity of the Sentinel-1 SAR sensor is low compared to our other sources of uncertainty (Fig. 5d thin vertical bars).

Considering the likely partitioning between retention and runoff at the visible runoff limit, our best-estimate scenario yields 84% of total melt refrozen as SI, of which 65% occurred during summer. The remaining 16% of melt likely ran off via the surface river network passing by FS2 (Fig. 1c), though we cannot rule out subsequent refreezing at lower elevations. Considering the stronger SIF scenario (Fig. 5d, right), 105% of total melt forms SIF, making this scenario unrealistic for our field area (SIF in excess of melt; precludes runoff, converse to indications from optical satellite imagery), while in the weaker SIF scenario (Fig. 5d, left) then only 71% of melt is refrozen.

We have established that refreezing onto ice slabs consists of at least two mechanisms: bottom-up SIF driven by slab cold content, and top-down SIF during autumn. A third mechanism is implicated by localised ice slab fractures that can enable melt to refreeze on the underside of the slabs 29 which we do not consider further here. Now we constrain at first order the importance of SIF onto the surface of ice slabs.

First we consider ice-slabs-wide autumn residual wetness. We used passive microwave radiometry 30 to identify the last day of the melt season (Methods). We used our SAR approach introduced above to identify the last wet day (Methods), yielding the annual duration of residual wetness (Fig. 6). Like in our field study area, residual wetness duration is effectively insensitive to the radiometric sensitivity of the Sentinel-1 SAR sensor, apart from minor sensitivity in 2022 (Suppl. Fig. 2). Residual wetness primarily varies according to the surface drainage network structure. On the western flanks remnants of slush fields and surface rivers are visible (Fig. 6b, c). Undrained lakes are particularly evident towards the lower slab extents. In the north and north-east, residual wetness is less pervasive and mainly consists of isolated wet patches (Fig. 6d, e), while in the south-east some residual wetness is detected around the upper margins of firn aquifers (Fig. 6a; ref. 12). The duration and extent of residual wetness around the annual visible runoff limit varied substantially by year (Suppl. Fig. 3), with more extensive and longer-duration wetness generally associated with a higher visible runoff limit. The western margins exhibited the most extensive wetness, extending over up to 80% of the area around the visible runoff limit in the several days following the end of the melt season, but as little as ~10% in years when the visible runoff limit was lower than average. Extensive wetness can persist for up to roughly 30 days; beyond this, wetness tends to persist only in more isolated patches.

Next, we estimated summer and autumn SIF across the ‘most likely’, ‘minimum’ and ‘maximum’ ice slab extents 12. For summer we used our SIF model in its bottom-up configuration, initialised with spatially-variable 10 m depth July firn temperatures (T 10m) 28 (Methods). To assess the uncertainty associated with the initial temperature profiles, we also ran the model for the uncertainty bounds of T 10m. We determined the start of summer SIF by using the regional climate model MAR 3 to identify when surface melt first became available for SIF (Methods). In autumn, we used our wetness duration retrievals (excluding those areas which remained wet at the end of the calendar year, which are generally undrained supraglacial lakes) in conjunction with our model to estimate bottom-up and top-down SIF (Methods). We also considered the possibility of SIF due to top-down atmospheric cooling during colder summers with prolonged melt hiatuses. We found that at first order, this process can be neglected compared with the other modes of SIF production, so we do not consider it further (Supplementary Discussion).

Between 2017 and 2022 over the ‘most likely’ ice slab extents, summer top-down SIF occurred at an average rate of 7.2 (6.4–8.2) Gt yr-1 (Fig. 7b; bracketed values calculated from T 10m uncertainty bounds). On average autumn bottom-up SIF occurred at 1.2 (1.0–1.4) Gt yr-1 (Fig. 7c). Autumn top-down SIF is harder to assess because it is highly sensitive to the thickness of the snowpack overlying the slush matrix, yet this thickness cannot be adequately retrieved from MAR (Methods). In lieu we computed autumn top-down SIF under three fixed snow thicknesses: 0.5, 1.0 and 1.5 m, yielding on average 2.8, 0.9 and 0.4 Gt yr-1 respectively (Fig. 7d), so its maximum contribution is much less than half the average SIF which occurs in summer. Considering all SIF modes, total SIF occurred at 4.3, 9.3 or 11.6 Gt yr-1 across the minimum, most likely and maximum slab extents respectively when using predicted T 10m and a 1.0 m autumn snowpack.

Total SIF varied considerably between years (Fig. 7e). Over the ‘most likely’ slab extents, considering the high-melt year of 2019 31 SIF was divided between 5.3 Gt during the melt season and 2.3 Gt during autumn (assuming a 1.0 m snowpack), totalling 7.6 Gt. In contrast, 2022 was a more moderate melt year but with a large melt pulse in early September; 10.0 Gt SIF formed during the melt season and 2.3 Gt in autumn, totalling 12.3 Gt, 57% more than in 2019. This difference is caused mainly by summer SIF, which in 2019 was reduced because of the complete melting of seasonal snowpack up to elevations approaching the visible runoff limit. This caused some of the SIF which had formed earlier in the season to be re-melted, and highlights that SIF on ice slabs can only occur under moderate melt conditions. If melt removes all accumulation, exposing the slab itself to direct melting, then retention will be depleted, such as in the extraordinary melt year of 2012 32. Such a reduction in slab retention during strong-melt summers will be accompanied by ice layer formation at elevations above existing slabs, laying the foundations for future slab expansion 11,33.

We quantified SIF by ice sheet drainage basin (Fig. 7f). Slabs along the western flank refroze the most meltwater, which is unsurprising given their large spatial extent. Despite sometimes extensive residual wetness, autumn refreezing accounts for at most one third of SIF in most basins, except the south-east, where autumn SIF apparently dominates; here, ref. 12 found possible ice slabs along the upper margins of firn aquifers, which may be proto-aquifers with a similar tendency to refreeze residual water slowly over the autumn.

Methods
Surface hydrological catchments
We used SAGA GIS 41 to delineate the hydrological catchments of our field sites FS1, FS2, and FS3 from Arctic DEM 42 to the uppermost limits of the underlying ice slab 11. The surface catchments of FS1 and FS3 are very small (1.1 and 0.1 km2 respectively), so we considered the catchment area upstream of FS2 only, which has an area up to the ice slab limits of 175 km2.

Measuring in-situ superimposed ice formation
In April 2021 (before the melt season) we dug pits and buried metal meshes to set marker horizons at three locations (Fig. 1a). In April 2022 we returned to measure ice formation on top of them by digging down to the new ice slab surface and then using an auger to drill into the ice until metal chips appeared at the top of the hole. This method is subject to uncertainties which are likely several centimetres of underestimation, in particular at FS3 where we had difficulties determining the April 2021 slab surface without destroying it due to the presence of mixed firn and ice lenses on top of the slab.

Between April 2022 and April 2023, we instead employed a repeat shallow core method at eleven locations. Using a Kovacs Mark II coring system, we drilled a 9 cm diameter core from the surface down through the snow and any firn (if present), until approximately 0.5 m into the slab (Suppl. Fig. 1a). We measured the depth of the borehole independently of the retrieved core length to check for core collapsing. One third of the cores had short collapsed sections, on average 6 cm long (max. 10 cm) and always situated in brittle depth hoar. We noted the core stratigraphy and measured density at 10 cm resolution. For the collapsed sections, we used the average measured density of depth hoar. We then fixed a bamboo stake in direct contact with the ice at the base of the borehole and measured the length of the stake which remained exposed (Suppl. Fig. 1b). One year later we visited each site, that is the location of the stake which had moved due to ice flow (~60 m yr−1 in the study area). We measured the distance from the top of the stake to the snow surface. Then we repeated the drilling in close vicinity of the stake (<1 m; Suppl. Fig. 1c) and we again noted stratigraphy, measured densities at 10 cm resolution and measured the depth of the hole independently of the core. Key parameters in the density measurements such as core diameter were checked throughout logging.

We used the top of each stake to align the two density profiles (Suppl. Fig. 1d) to (i) quantify SIF in metres as well as metres water equivalent, and (ii) quantify local climatic surface mass balance (cSMB). The latter was calculated by directly subtracting the water equivalent measured in spring one and spring two. For the purpose of the subtraction, water equivalent was summed up over the depth range from the snow surface in year two down to the ice slab surface in year one.

We drilled repeat cores, which indicated that our measurements were repeatable and that metres scale variability in SIF was negligible (Supplementary Discussion). We also considered the impact of springtime densification on SIF as a potential source of systematic errors (Supplementary Discussion).

Measuring slab temperatures
At FS2, in May 2022 we installed a 12 m digital temperature chain with 34 sensors, with spacing increasing from 0.15 m at the top of the chain to 1 m at the bottom (beaded stream, inc). The chain was installed in a borehole of ~11 cm diameter drilled by a Kovacs Mark II coring system and suspended from a boom attached to a mast frozen into the sub-surface. This enabled us to capture not only slab temperatures but also seasonal snowpack temperatures with respect to a fixed datum. Measurements were logged on a D605 logger (beaded stream, inc) every 15 minutes until recording stopped unexpectedly on 18 September at 04:00 UTC.

Melt available for runoff or SIF in field area
We calculated meltwater production in our study area using a surface energy balance/surface mass balance (SEB/SMB) model forced by local observations from three weather stations, including KAN_U, supplemented with MODIS albedo 25. The SEB/SMB model calculates surface energy and mass balance components in 100 m elevation bins.

The model’s sub-surface scheme (0.25 m resolution) consists of a seasonal snowpack overlying firn and ice layers generated since the model initialisation in spring 2009. The model calculates surface melt refreezing in underlying layers with pore space and cold content. At the start of the melt season, meltwater refreezes locally, warming and densifying the snowpack, and no runoff occurs. Once the overlying snowpack is isothermal then water percolates down to the underlying modelled ice layers, which are defined as impermeable when >1 m thick. They are present throughout our study’s elevation range, causing instantaneous runoff. We therefore estimate the meltwater available for SIF or runoff as that produced between the start of modelled runoff and the end of the melt season, less the refreezing at 0.25 m depth which captures diurnal refreezing as the near-surface cools at night.

We assessed model performance by comparison with observations of surface lowering at FS1, FS2 and FS3 which were each equipped with a CR1000X datalogger (Campbell Scientific; CS) and a CS SR50A sonic ranger mounted on a horizontal boom roughly 1–2 m above the snow to measure surface height changes. In 2022 at 1800 m a.s.l., the model produced 0.25 m w.e. melt available for runoff or SIF. This is in excellent agreement with our observations: using a fixed summer-time density of 504 kg m-3 (ref. 13) to convert surface lowering to water equivalent, 0.23 to 0.27 m w.e. melt occurred in the same time period.

To interpolate melt available for runoff or SIF in 2022 over our hydrological catchment we used linear regression between elevation and melt modelled in the 100 m elevation bins between 1600 and 1,900 m (R2 = 0.99, p < 0.01, n = 4) to map melt production as a function of elevation.

Modelling superimposed ice formation
There have been limited previous efforts to simulate SIF by slush refreezing on top of an ice slab, principally addressing only bottom-up refreezing and neglecting any impact of lateral flow of meltwater 15,43,44. The model by ref. 45 addresses the effect of lateral meltwater flow in SIF but lacks top-down refreezing and uses an experimental flow-line approach. Most regional climate models (RCMs) applied to the Greenland Ice Sheet do not have an ability to simulate slush and they are all restricted to 1D simulations along the depth axis. Lateral flow of meltwater and ponding, important in autumn SIF, therefore cannot be simulated. Furthermore, RCM simulated firn structure and firn temperatures in the lower accumulation area, essential to reliably simulate SIF, are of limited accuracy 28,46. In summary there are currently no models suitable to simulate SIF on an ice-sheet-wide scale.

We aimed to make a first-order estimate of Greenland-wide SIF onto ice slabs and to assess the relative importance of bottom-up and top-down SIF. We do not claim to capture the full complexity of the processes involved. Here we describe the model formulation and its application to our field study area. Its ice-slabs-wide application is described separately in ‘Ice-slabs-wide wetness and superimposed ice formation’, below.

Our 1-dimensional heat flow model computes SIF based on the amount of latent heat evacuated from slush. The heat flux follows the temperature gradient between the slush (0 °C) and the cold slab in the case of bottom-up refreezing, or between the slush and cold snow during top-down refreezing (Fig. 3). We assume that heat conduction is the only mechanism of heat flux. Heat conduction is calculated numerically in time steps of 200 seconds and along the depth axis in 10 cm layers. We use a density of 920 kg m-3 for ice and a density of 400 kg m-3 for snow with a porosity of 0.4 and 5% initial irreducible water content. This value of irreducible water content is located between the classical value of 7% indicated by ref. 47 and values used in regional climate models, such as MAR where it varies with depth from 7% at the surface to 2% at 1 m below the surface. Thermal conductivity was set to 2.25 W m-1 K-1 for ice and 0.5 W m-1 K-1 for snow at 400 kg m-3 (refs. 48,49). We used a specific heat capacity of ice of 2090 J kg-1 K-1 and latent heat of 334,000 J kg-1.

To model bottom-up SIF we assume that there is always slush at 0 °C on top of the impermeable ice slab. Slush is continuously refrozen as its latent heat is released and conducted to the ice slab. Temperatures at the bottom of the model domain are allowed to fluctuate without prescribing a fixed bottom heat flux; at every time step, the temperature change of the bottom-most grid cell equals the temperature change in the cell above. We specified our bottom-up model domain and initialisation according to thermistor string measurements acquired at FS2, where the slab is >12 m thick 13. The model domain was 12 m deep, with the interface between slush and ice slab located at 0 m depth. The model was initialised with a measured firn temperature profile from 6 July 2022 14:15 UTC, the moment where the thermistor closest to the ice slab surface reached ~0 °C. After initialisation, the model ran independently of measured firn temperatures. We then partitioned bottom-up SIF into its summer and autumn components according to the last day of the summer melt season and the last day of residual wetness (according to SAR acquisitions, below).

We model top-down SIF separately to bottom-up SIF, for two major reasons. First, the thickness of the slush layer is unknown, especially during autumn when both processes of SIF are active at the same time. By modelling the two processes separately, we avoid the need to prescribe the slush layer thickness. Second, the snowpack overlying the layer of slush is also poorly constrained (Fig. 3). In our field study area we evaluated the sensitivity of the model results to snowpack thickness overlying slush by running our calculations for the observed snowpack thickness of 0.6 m as well as snowpacks 0.5× and 1.5× thinner and thicker respectively. Initial irreducible water content was taken to be the same as for bottom-up SIF. The top model boundary was taken to be the snow surface where the temperature therefore varied according to the surface temperatures calculated by the SEB/SMB model. The bottom boundary was fixed at 0 °C. Model runs started following the last day of surface melting and ended on the last day of residual wetness.

Remote sensing of autumn residual wetness
Theoretical radar backscatter response to sub-surface meltwater. Horizontally-transmitted, vertically-received (HV) backscatter intensity observations (hereafter σHV0) made by 5.4 GHz radar are highly sensitive to the dielectric contrast between water and snow/ice and penetrate up to a few metres of ice under dry conditions 50. σHV0 is therefore indicative of water content in the uppermost metres of the near-surface snow and firn.

To further verify the suitability of our approach to detect sub-surface wetness we carried out a suite of radiative transfer (RT) modelling experiments in the Snow Microwave Radiative Transfer (SMRT) model 51 using the sticky hard sphere (SHS) model (Supplementary Methods). Our results indicate that variability in the liquid water content of the sub-surface layer dominates the reduction in the modelled σHV0 (Supplementary Discussion).

Synthetic Aperture Radar data. We used Synthetic Aperture Radar (SAR) backscatter intensity data acquired by the European Space Agency Sentinel-1 platform. Data availability prior to 2017 is insufficient for this study. Since 2017/2018 the acquisition strategy favours more frequent data acquisition. We used Interferometric Wide mode σHV0 data, delivered as a ground-range-detected (GRD) product at 25 m resolution. Acquisitions were more frequent from 2017 until the end of 2021 when Sentinel-1B failed. All data were pre-processed following standard procedures (Supplementary Methods).

Detecting residual water presence. We detected residual water in the upper several metres of the snow/firn/ice column following the end of the melt season until the end of the calendar year by examining successive autumn SAR acquisitions co-normalised by subtracting winter reference scenes of the corresponding satellite orbit to yield σbHV0 (ref. 27). This approach reduces variations in σHV0 introduced by observation geometry, terrain fluctuations and melt conditions between images acquired at different times. Moreover, it allows to account for spatial variations in surface and subsurface heterogeneities (e.g. surface roughness, sastrugi, ice lenses, ice pipes, percolation pipes) that affect the spatial variations in σHV0. Therefore, we assume that these heterogeneities are temporally non-varying, such that changes in the cross-polarization backscatter can be attributed solely to changes in dielectric permittivity rather than surface roughness and volume scattering 52.

We generated a winter reference scene for each orbit by calculating the median of all data acquired between 1 March and 16 April each year between 2017 and 2022. Next, we co-normalised scenes acquired between July and December each year using the corresponding reference scene. When the snow/firn column is completely frozen the σbHV0 distribution is zero-centred with little variance (Suppl. Fig. 4, green curve). Conversely, when water is present in the column it can be detected as an anomaly relative to background σbHV0 (Suppl. Fig. 4, purple curve.).

To identify the threshold at which σbHV0 pixels can be considered wet we used the 1st percentile of σbHV0 in May 2021 (i.e. prior to surface melting) across all ice slabs as mapped by 12 (Suppl. Fig. 4, green curve). We added the Sentinel-1’s radiometric sensitivity of -1 dB, yielding a wetness detection threshold T wet < -2.17 dB (Suppl. Fig. 4, blue vertical line). To evaluate uncertainty in our wetness durations we varied our threshold by ± 0.3 dB, which is the stated 1σ precision of the Sentinel-1 sensor 53.

Ice-slabs-wide wetness and superimposed ice formation
Ice slab extent. We used observationally-derived ice slab extents at 500 m resolution 12. This data product distinguishes between the ‘most likely’, ‘minimum’ and ‘maximum’ slab extents. We estimated SIF over each of these extents during the period 2017 to 2022.

Identifying the end of the melt season. Passive microwave radiometry (PMR) is widely adopted to detect surface melting over the Greenland Ice Sheet. We used Ka band (37 GHz) horizontally polarized brightness temperatures measured on a daily basis by the SSMI/S-F17 sensor and binned at a spatial resolution of 3.125 km 54. At this frequency PMR interacts with the uppermost few centimetres of the ice sheet surface and is thus indicative purely of surficial water status. Surface melting is detected for a pixel when the measured brightness temperature is higher than a threshold value. The threshold is dynamically computed for each pixel and for every year according to the simulations of a microwave emission model 30. We considered all data up to 30 September each year. In each pixel we identified the last day of surface melting each year.

Autumn residual wetness. We defined the residual wetness duration as between the end of surface melting (PMR) and the last date of residual wetness (SAR).

We analysed ascending and descending orbits of Sentinel-1 separately. We preferentially used data from ascending orbits to mitigate against conflicting geo-location errors present in the Google Earth Engine version of the Sentinel-1 data (Supplementary Methods). This approach limits the temporal resolution to the relative orbit of the chosen orbit direction: 5 days when Sentinel-1B was active, otherwise 11 days. We multi-looked to 500 m resolution, matching the ice slab extents dataset.

We excluded two short periods of SAR data as follows: 26 Nov 2021 and 28–31 Oct 2022. During both these periods temperatures rose briefly to ~0 °C, causing identification of surface melting which was not indicative of multi-day sub-surface wetness. The snowpack had already undergone several weeks of atmospheric cooling, making it highly likely that any surface melt refroze in the upper snowpack before it could reach the slab surface.

Bottom-up SIF. We calculated the slab-wide amount of bottom-up SIF in summer and autumn by running our 1-D heat flow model over the 10 km grid of the regional climate model MAR and then resolving the end of residual wetness using SAR measurements at 500 m resolution.

To simplify computation over all ice slabs and because the slabs extent dataset by ref. 12 does not also provide slab thickness, we assumed that slabs are 12 m thick like at FS2. This value is roughly in the middle of the observed ice-sheet-wide slabs thickness range: well-developed slabs have a median thickness of 15 m, while slabs in development or initiation have a median thickness of around 5–11 m 11. We checked the sensitivity of our model to slab thickness by repeating the model run at FS2 (Fig. 2) with slab thicknesses of 5 and 15 m. Over this thickness range, SIF at the end of the melt season on 18 September varied by < 2 mm w.e. and so we conclude that prescribing slab thickness of 12 m has a negligible impact on our results.

We specified the initial slab temperature profile in each grid cell based on observation-derived monthly 10 m firn temperatures (T 10m) interpolated across the ice sheet at 0.1° resolution 28. We used this dataset to scale our initialisation profile measured at FS2 so that T 10m agreed with the closest value T 10m local in July each year. Temperatures at all depths were then multiplied by T 10m local / T 10m FS2 which preserves the condition of 0 °C at the top of the ice slab. To assess the sensitivity of bottom-up SIF to the initial temperature profile, we ran the model using T 10m ± their associated per-pixel uncertainty estimate.

The simulation period was individual to each grid cell. The simulation commenced when meltwater was inferred to percolate to the top of the slab for the first time. This date is logically always some time after surface melt begins, so we cannot use PMR for this purpose. Instead we used MAR v3.14.0 37 forced by ERA-5 55 to identify the date on which modelled runoff from beneath the seasonal snowpack commenced. MAR simulates solid ice layers below seasonal snowpack in all ice slab areas and modelled runoff therefore begins once meltwater has warmed the snowpack and percolated down to the ice layer.

In principle the bottom-up SIF attributable to summer ceases at the end of the melt season determined by PMR. However, we also needed to account for the possibility that the snowpack melts away before the end of the melt season, exposing the slab at the surface, which can be relevant towards the lower elevations of slabs. If a cell became snow-free according to MAR then we halted summer-time bottom-up SIF accordingly and subtracted any subsequent melting from the already-formed SI.

We calculated autumn bottom-up SIF by allowing the bottom-up model to continue into autumn beyond the end of surface melting for as long as residual wetness determined by SAR (at 500 m resolution) was present.

Top-down SIF. Like bottom-up SIF we ran the top-down SIF model over the MAR grid. We initialised the model with an isothermal 0 °C snowpack with 5% irreducible water content then forced it with MAR v3.14.0 surface temperatures at the topmost snow layer in each numerical time step. As soon as subzero temperatures reached the snow-to-slab interface then top-down SIF could begin.

In each cell the simulation started the day after the end of the melt season according to PMR. We ran the simulation until the end of the year. Then, we used the residual wetness duration measured by SAR to determine the end of top-down SIF at 500 m resolution. The thickness of the overlying snowpack was held constant, which means that snowfall events were not included. We considered using MAR end-of-season snowpack thicknesses but found them to be spatially incoherent, so instead we investigated SIF under 0.5, 1.0 and 1.5 m of overlying snowpack to provide bounds on model sensitivity.

Top-down autumn SIF is likely to be more complicated than our first-order simplified approach. For instance, although rarely observed, slush fields can suddenly appear at the surface during cold periods, which may have similarities to aufeis formation in Arctic terrain 56–58. Further consideration is beyond the scope of the present study.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript presents estimates of superimposed ice formation on top of ice slabs along the K-Transect in SW Greenland. In this study, the authors present field measurements that suggest local ice formation on ice slabs can retain ¾ of local annual surface melt. They estimate that ice formation along ice slabs facilitated the retention of ~32 gigatonnes of melt between 2017 and 2022, roughly ~4% of the Greenland SMB.\nOverall, the manuscript is well-written, thorough, and presents important observational evidence for surface melt refreeze onto ice slabs, which plays a role on the overall SMB of the GrIS. I have several comments, most of which are relatively minor; I believe that their consideration will overall improve the manuscript.\nI have organized the comments below, grouped by section within the manuscript.\nResults:\nPg 3 line 77: Because local melt is near uniform across the study area (presumably as a result of their modeling results, but this is not made explicit), they hypothesis that measured spatial variation in SIF is indicative of inputs not only from local melt production, but also lateral water flow. This seems entirely reasonable given the observed increases in SIF for sites proximal to major surface rivers. However, the inclusion of external references to reinforce this claim would be beneficial.\nPage 4 line 79: “Finally, beyond 2–4 km away from visible surface runoff, we found that SIF began to increase again, as illustrated by 0.25 m w.e. at site 3u3 (Suppl. Fig. 1).” – With what degree of confidence can this claim be made, given that it is reliant upon a single point wise measurement of SIF, and the measured value is within 2 standard deviations of the reported average SIF from 2022?\nPage 7 line 139: “In contrast, larger amounts of SIF can only be explained by lateral transport of meltwater which provides additional water for SIF to continue in specific locations beyond the end of the melt season.” – Similar to first comment; Isn’t it possible that this elevated SIF can also be explained by spatial heterogeneity in melt that is not captured in the SMB model? I tend to agree with the authors that the more plausible explanation is lateral transport of meltwater, but to what degree can spatial heterogeneity of melt be ruled out?\nAn important point to mention in the refreezing during Autumn section is the role that volume scattering due to surface and subsurface heterogeneities (e.g. surface roughness, sastrugi, ice lenses, ice pipes, percolation pipes) plays. Implicitly, the bias-removal approach that is employed here *should* account for this process, provided the geometry of these features is assumed temporally non-varying, such that changes in the cross-polarization backscatter can be attributed solely to changes in dielectric permittivity rather than surface roughness and volume scattering. I would suggest that this point be elaborated upon in the methods section, perhaps with reference to some relevant literature (e.g., Rignot 1995, Backscatter Model for the Unusual Radar Properties of the Greenland Ice Sheet).\nMethods:\nPage 16, line 396: Regarding the model producing 0.25 m w.e. melt; is this merely a consequence of the subsurface model discretization being 0.25 m? If not, what is the uncertainty of this model-derived estimate of w.e. melt?\nPage 18 line 469: Section on the linear regression between measured SIF and mean crosspol backscatter: given that an uncertainty in the in situ SIF values are provided, it would be worthwhile to propagate this uncertainty through the linear regression for a nominal uncertainty estimate of the linearly-regressed area-wide SIF estimates, particularly since the R2 in the linear regression isn’t particularly strong.", "reviewer_2": "Reviewer #1 (Remarks to the Author):\nThis manuscript presents estimates of superimposed ice formation on top of ice slabs along the K-Transect in\nSW Greenland. In this study, the authors present field measurements that suggest local ice formation on ice\nslabs can retain ¾ of local annual surface melt. They estimate that ice formation along ice slabs facilitated the\nretention of ~32 gigatonnes of melt between 2017 and 2022, roughly ~4% of the Greenland SMB.\nOverall, the manuscript is well-written, thorough, and presents important observational evidence for surface\nmelt refreeze onto ice slabs, which plays a role on the overall SMB of the GrIS. I have several comments,\nmost of which are relatively minor; I believe that their consideration will overall improve the manuscript.\nI have organized the comments below, grouped by section within the manuscript.\nResults:\nPg 3 line 77: Because local melt is near uniform across the study area (presumably as a result of their\nmodeling results, but this is not made explicit), they hypothesis that measured spatial variation in SIF is\nindicative of inputs not only from local melt production, but also lateral water flow. This seems entirely\nreasonable given the observed increases in SIF for sites proximal to major surface rivers. However, the\ninclusion of external references to reinforce this claim would be beneficial.\nAs the reviewer infers, we take melt as a smoothly varying field of decreasing melt with increasing elevation,\nas calculated by a surface energy balance model forced by in-situ observations. Fig 5d now illustrates how\nmelt decreases with increasing elevation. We also now refer to findings from Noel et al. 2016 in the first\nparagraph of our Results. They showed that melt rates show very little spatial variability in this region.\nIt is possible that, once visible slush fields have developed, then local melt rates might increase further as a\nfunction of their reduced albedo. However, these slush fields can only initially develop by the transport and\nponding of melt that is essentially uniform in space considering the smooth surface topography over 100s\nmetres scales and homogeneous snow surface cover (e.g. Holmes 1955; Noel et al., 2016; Clerx et al., 2024/in\npress). In other words, their origin does not lie in albedo. Nonetheless, we now note this possibility of\nsubsequent local melt enhancement in the first paragraph of the Discussion.\nPage 4 line 79: “Finally, beyond 2–4 km away from visible surface runoff, we found that SIF began to\nincrease again, as illustrated by 0.25 m w.e. at site 3u3 (Suppl. Fig. 1).” – With what degree of confidence can\nthis claim be made, given that it is reliant upon a single point wise measurement of SIF, and the measured\nvalue is within 2 standard deviations of the reported average SIF from 2022?\nOn reflection, we agree, so we have removed this statement from the main text.\nPage 7 line 139: “In contrast, larger amounts of SIF can only be explained by lateral transport of meltwater\nwhich provides additional water for SIF to continue in specific locations beyond the end of the melt season.” –\nSimilar to first comment; Isn’t it possible that this elevated SIF can also be explained by spatial heterogeneity\nin melt that is not captured in the SMB model? I tend to agree with the authors that the more plausible\nexplanation is lateral transport of meltwater, but to what degree can spatial heterogeneity of melt be ruled out?\nPlease see our response to the previous comment about spatial variability in melt.\nAn important point to mention in the refreezing during Autumn section is the role that volume scattering due\nto surface and subsurface heterogeneities (e.g. surface roughness, sastrugi, ice lenses, ice pipes, percolation\npipes) plays. Implicitly, the bias-removal approach that is employed here *should* account for this process,\nprovided the geometry of these features is assumed temporally non-varying, such that changes in the cross-\npolarization backscatter can be attributed solely to changes in dielectric permittivity rather than surface\nroughness and volume scattering. I would suggest that this point be elaborated upon in the methods section,\nperhaps with reference to some relevant literature (e.g., Rignot 1995, Backscatter Model for the Unusual\nRadar Properties of the Greenland Ice Sheet).\nThanks for this suggestion. We have elaborated in the Methods (‘Detecting residual water presence’)\naccordingly.\nMethods:\nPage 16, line 396: Regarding the model producing 0.25 m w.e. melt; is this merely a consequence of the\nsubsurface model discretization being 0.25 m? If not, what is the uncertainty of this model-derived estimate of\nw.e. melt?\nNo, this is not a consequence of the sub-surface model discretisation. In each time step, the model’s grid\nmoves up/down after the occurrence of accumulation/ablation; van As et al. (2017) contains full details.\nConcerning uncertainty: the model does not come with uncertainty bounds because of the complex and non-\nlinear nature of the modelling. This makes validation a crucial aspect of SEB modelling, as is done in\nnumerous previous studies. Concerning this model specifically, see e.g. van As et al. (2017) cited in\nmanuscript. In this manuscript we therefore evaluate how sensible its melt estimates are by reference to in-situ\nsurface lowering measurements, which show excellent agreement (see Methods).\nPage 18 line 469: Section on the linear regression between measured SIF and mean crosspol backscatter:\ngiven that an uncertainty in the in situ SIF values are provided, it would be worthwhile to propagate this\nuncertainty through the linear regression for a nominal uncertainty estimate of the linearly-regressed area-\nwide SIF estimates, particularly since the R2 in the linear regression isn’t particularly strong.\nThe linear regression approach has been removed in this re-submission.\n\nReviewer #2\nThis manuscript describes field measurements of superimposed ice formation on top of ice slabs on\nGreenland’s K-transect in the summers of 2021 and 2022. Using repeated firn cores, time domain\nreflectometry, and thermal modeling, the authors show that a significant portion of annual melt refreezes on\ntop of the ice slabs due to ice slab cold content and cold surface temperatures during autumn. In particular, the\nauthors show that refreezing occurs even where water flows laterally over the surface, potentially for kms,\nsuggesting that surface mass loss may not scale directly with visible surface runoff during the melt season.\nThe authors then use Sentinel-1 SAR measurements to establish relationships between their field\nmeasurements of SIF and radar-inferred surface wetness measurements. They use these relationships to scale\ntheir SIF measurements first to the area surrounding their field site and then to the entire ice sheet. The\nauthors conclude that SIF accounts for about 4% of the contemporary GrIS SMB, making it a small but non-\ntrivial contributor to the overall ice sheet mass balance.\nOverall, this paper addresses an important caveat to the popular notion that ice slabs lead to complete surface\nrunoff. In particular, the field measurements are novel and convincing and suggest that superimposed ice\nformation (SIF) needs to be accounted for in order to understand how quickly surface melting becomes mass\nloss as melting advances inland on the ice sheet. However, the remote sensing analyses used to scale these\nfield measurements to the entire ice sheet are subject to large uncertainties, many of which are not well\nquantified in the current manuscript. I am particularly concerned that the approach to estimating pan-\nGreenland SIF relies on a highly uncertain relationship between SIF and nCDD in a single ~30 km x 30 km\nregion which cannot be validated anywhere else on the ice sheet, as well a background SIF rate that is\nassumed to not vary at all across the ice sheet. This opens the potential for large, uncharacterizable regional\nbiases to affect the pan-Greenland SIF estimates.\n‘[the authors use a] background SIF rate that is assumed to not vary at all across the ice sheet’. In fact, the\nbackground SIF rate used for the ice-slabs-wide analysis varies in two ways:\n(1) As in our first submission, our treatment of the conditions which permit SIF across the ice sheet is\nexplicitly spatial. At a given MAR model grid node, SIF can only occur during the melt season and\nwhen the slab is covered in (porous) snow which provides the basis for a slush matrix overlying the\nslab.\n(2) Previously, we used three different fixed initial temperature scenarios to model bottom-up summer\nSIF. In this re-submission, we instead use the state-of-the-art ice-sheet-wide 10 m depth\ntemperatures data product published by Vandecrux et al. (2024, TC) to initialise our bottom-up\nSIF modelling. We also run our model at +/- the stated uncertainty of the Vandecrux product.\nWe consider this to be a major improvement compared to the previous version of our manuscript.\nConcerning the remote sensing analyses: we concur with the reviewer’s concerns that these were subject to\nlarge uncertainties and had the potential for large, uncharacterizable regional biases to affect the pan-GrIS SIF\nestimates. In conjunction with similar concerns raised by reviewer 3, we have therefore made major\nrevisions which replace the previous statistical approach completely with physical heat flow modelling\nto make our upscaling calculations. In summary:\n(1) In areas which remain wet after the end of the melt season according to SAR, we now allow our\nsummertime “bottom-up” calculations to continue into autumn until SAR indicates that the area is\ndry, yielding the additional bottom-up SIF which occurs in autumn. We partition the contributions of\nsummer vs autumn bottom-up SIF using the last day of melting from passive microwave radiometry.\n(2) In these same areas, we now also explicitly address “top-down” refreezing as a result of negative\nsurface temperatures, beginning once the melt season has finished. These calculations are very\nsensitive to the thickness of the snowpack overlying the slush matrix; we discuss this further in the\nmanuscript.\n(3) Any pixels suspected of being lakes (i.e. still ‘residually wet’ at the end of the year) are now\ncompletely excluded from our calculations.\n(4) We have recalculated the SAR sensitivity threshold using data acquired in all areas labelled as\n‘minimum likely’ slabs by Culberg et al. (2024, TC), with little resulting change (see below).\n(5) Over our field area, we now show the uncertainty in our SIF calculations associated with the stated\nsensitivity of the Sentinel-1 SAR sensor used to determine wetness duration. We show that this source\nof uncertainty is second order, i.e. wetness duration does not change appreciably within this\nsensitivity range. We also checked this on an ice-slabs-wide basis where again it is a second-order\nsource of uncertainty (see Supplementary Figure 4).\nMajor Comments:\n[1] I am a bit confused about the primary purpose of modelling of the rate of superimposed ice formation. The\nauthors state that the reason they conduct the modelling is to determine the degree to which ice slab warming\nis due to latent heating. If this is indeed the main goal, then the argument would be better supported by\ncomparing the model run presented in the paper to a model run forced by the measured surface temperature\nand which only considers heat conduction (no latent heating). The difference between those two scenarios\nwould provide a stronger quantification of the role of latent heating in warming the ice slabs.\nHowever, the primary conclusion that seems to come out of this section is that refreezing due\nto ice slab cold content can only account for measured SIF rates in regions without additional\nadvection of meltwater. If the goal is instead to accurately model local SIF, then it seems like the authors\nmight be better served by using their fully coupled SEB-firn model instead of this simplified model that\nimposes a fixed meltwater supply that may or may not be consistent with the actual melt supply and timing\nduring the season when measurements were made.\nWe agree that we should have been clearer here. Our main goal was neither to super-accurately model local\nSIF nor to determine the degree to which ice slab warming is due to latent heating. Rather, it was to provide a\nconstraint on how much SIF can occur during the melt season under the plausible assumption of a continuous\nmeltwater supply at the slab surface. Then, SIF in excess of this logically occurs during autumn, i.e. after\nactive surface melt has finished. In other words, the primary purpose of this modelling approach is to partition\nthe SIF that we measured into the melt and post-melt seasons, thereby setting the scene for a slabs-wide first-\norder estimate of SIF later in the study. With our revised modelling approach we believe that this is\nmotivation is much clearer.\nAn ideal solution would indeed be to use a fully coupled SEB-firn model. However, no such model that is\nsuitable for our application exists. We are dealing with the problem of ponded meltwater in a snow matrix\natop an impermeable and cold barrier, which is an extreme end-member of what can be considered ‘firn’. We\nunderstand from firn modelling colleagues that even when ponding subroutines are available then they often\nremain disabled for want of validation. We now describe these issues in the Methods section ‘Modelling\nSIF’. Thus, it is more informative to simplify the problem as far as possible, i.e. by considering a relatively\nthin snowpack with a continuously wet slush matrix on top of an ice slab (conditions which were met for\nmuch of the melt season according to our in-situ measurements and further bolstered by observations in Clerx\net al., 2022) and then using classical heat flux theory to determine the rate of refreezing at first order.\nNew in this re-submission is an alternative configuration of the model to estimate SIF caused by\natmospheric cooling (‘top-down’). As we describe in the paper, an issue here is that we lack good bounds on\nthe thickness of the snowpack overlying the slush matrix. We wanted to use end-of-season snowpack\nthickness simulated by MAR but found that the outputs were highly spatially incoherent, which is a further\ncaution against 'blind’ use of RCMs assuming that they always deliver the best results. The reviewer may be\ninterested to see a paper with several of the same co-authors as this submission currently in Discussion, which\ndescribes some of the difficulties with RCM behaviour at the runoff limits (https://doi.org/10.5194/egusphere-\n2024-2750).\nOverall, we contend that our chosen approach is much simpler to interpret than the outputs of a fully coupled\nSEB-firn model. In turn, we hope that our findings can be used to advance ponding schemes in firn models.\nWe have made narrative changes to the relevant parts of the manuscript accordingly.\n[2] Uncertainty in the regression between measured SIF and the multi-temporal mean 𝜎̂ 𝐻𝑉 can be better\nquantified and mapped into uncertainty in the study site SIF map (Fig. 4b). One option would be to use a\nBayesian regression that can quantify the range of linear models that are consistent with the data. This would\nmake it possible to produce an uncertainty map showing the range of the plausible SIF for each pixel in Fig.\n4b given the range of plausible models. But the key point is that there should be corresponding uncertainty\nmap for Fig 4b that account for uncertainty in the SIF-𝜎̂ 𝐻𝑉 relationship and the uncertainty in the background\nSIF rate.\nSimilarly, better uncertainty quantification is needed for the SIF-nCDD regression. Particularly given the\nrelatively poor R2 value for this regression, I strongly feel that it is important to quantify the uncertainty in the\nmodel fit and clearly demonstrate how that uncertainty maps onto uncertainty in the final estimates of regional\nautumn SIF that are presented in Figure 7.\nThis analysis has been removed from our re-submission.\n[3] For setting the SAR wetness backscatter threshold, I would like to see some sensitivity tests from different\nregions of the ice sheet showing how the 1st percentile of 𝜎̂ 𝐻𝑉 prior to the melt season varies regionally. For\na reliable pan-Greenland mapping, it is very important to establish whether conditions in SW Greenland are\nconsistent across the ice sheet. Other studies that have used Sentinel-1 data to quantify either melting or the\nextent of wet snow facies over Greenland have found that performance with a single threshold can vary\nsignificantly across the ice sheet (Hu et al., 2022; Li et al., 2023). Regional bias due to a fixed wetness\nthreshold should be quantified and incorporated into assessments of overall uncertainty.\nWe thank the reviewer for raising these concerns. First, we comment on the literature citations:\n(1) The findings of Hu et al. 2022 are not as relevant to our study as they first appear. Hu et al. did not\nnormalise summertime retrievals with respect to winter, unlike in our current approach, and so\ntherefore a considerable source of variability remains.\n(2) Li et al. 2024’s methodology is more similar to ours. They emphasise that SAR-based mapping of\nwetness at a set threshold is only reliable above the snowline, but that there is no appreciable regional\nvariability above the snowline. As we consider only snow-covered surfaces above the ELA, the\ndifficulties that Li et al. experienced with selecting a threshold applicable to melt detection over bare\nice is not relevant here.\nWe recap that our threshold of -2.1 dB was quite similar to Liang et al. 2021, who rather than the 1st percentile\ntake the less conservative 5th percentile of dry snow as their threshold: -2.66 dB. Both Liang et al., 2021 and\nLi et al. 2024 employ a more conservative value of -3 dB, +/- 0.3 dB, but they had the express aim of\ncapturing surface melting conditions which they validated by reference to in-situ air temperature\nmeasurements. This is not a relevant validation metric in our case as we are seeking to identify residual\nwetness after the end of active melt, which also indicates that a 3 dB threshold would be too conservative for\nour purposes here.\nNotwithstanding this comparison, we have undertaken another examination of our threshold choice. Now, we\ncompute histograms by sampling all areas within the ‘minimum likely ice slab’ classification of Culberg\net al. (2024), rather than just the western flank of the ice sheet as previously. The 1st percentile value\nhas decreased slightly to -1.17 dB (i.e. 0.07 dB more conservative than in our previous submission).\nIn the revised manuscript we adopt the new threshold of -1.17 dB, +- 0.3 dB indicative of the\nradiometric uncertainty of Sentinel-1’s sensor as per Geudtner and Torres (2012). We use the +/- to\npropagate an uncertainty in wetness duration through to our final calculations in our field area. On an\nice-slabs-wide basis, we provide a new Supplementary Figure which shows that the residual wetness\nduration/area is essentially insensitive to the radiometric uncertainty of the sensor.\n[4] I am concerned about how well the nCDD model can be generalized to map autumn SIF rates across the\nentire ice sheet. This is conceptually quite similar to the positive degree day model used for melting, which\nhas a spatially and temporally variable rate factor that typically needs to be tuned for individual sites. I have\nsimilar concerns for the nCDD model. Essentially, the authors’ scaling assumes that air temperature is the\nonly driver of refreezing and neglects all the other components of the subsurface energy balance, to include\nfirn and snow conditions, which might vary significantly across the ice sheet. It also assumes that the amount\nof melt that is refrozen scales directly with the length of time required to fully refreeze an area, which may be\ngenerally reasonable, but would absolutely be influenced by the transient accumulation rate, subsurface\ntemperatures, etc on the local scale. Without secondary validation sites, it seems basically impossible to know\nif the regression model coefficients will apply anywhere outside of southwest Greenland.\nThis analysis has been removed from our re-submission.\n[5] How are buried lakes accounted for when calculating autumn SIF? If I have understood the methods\ncorrectly, it seems like there has been no attempt to separate that mode of water\nimpoundment from superimposed ice formation. It seems problematic to apply a relationship\nbetween nCDD and refreezing in porous snow or a slush matrix to the refreezing of large open water bodies.\nSince buried lakes will have the highest and most persistent wetness values (see the authors’ own Figure 5,\nwhich is dominated by lakes), I am concerned that including lakes may significantly bias the results.\nThanks for raising this concern. Now, we screen out likely lakes by removing pixels which are still wet in\nthe last SAR acquisition each calendar year.\nMinor Comments:\nOverall, the section titles and introductory sentences could do a better job of explaining to the reader the end\ngoal of the analyses being described in the section. For example, the section “Modelling the rate of\nsuperimposed ice formation” would be easier to follow if it started with a sentence explaining why the authors\nwant to model this quantity. The section title could also better reflect the findings instead of the methods,\nfor\nexample, something like “Ice slab cold content a key factor in summer SIF”.\nThanks for this suggestion. We have made several changes accordingly.\nLine 34 – I am not sure I understand what the second half of this sentence is trying to express.\nWe have rephrased as we agree that this was very unclear previously.\nLine 92 – Why was there no spike in 𝜖 when this water reached and fully flooded the probe before refreezing?\nWe cannot know for sure, but as the section ‘Timing of SIF’ has been removed in response to reviewer 3\nconcerns, we have not investigated further.\nLine 185/296 – there are now estimates of ice slab extent from both modeling (Brils et al., 2024) and Sentinel-\n1 observations (Culberg et al., 2024) that provide a more up-to-date and comprehensive outline of ice slab\nextent on the GrIS.\nIndeed - though we find it a little concerning that literature published in 2023 can already be considered less\n‘up-to-date’! Nonetheless, given that we have made major changes to our upscaling approach, we have taken\nthe opportunity to replace the Jullien et al. (2023) outlines with those based on Sentinel-1 observations\nby Culberg et al. (2024). We now explore the impact that the uncertainty in ice slab extent makes on our final\nestimates of SIF.\nWe have moved the statements of ice slab extent previously found at this location into the Introduction.\nLine 223 – modeled temperature profile, for example, from IMAU FDM with ice slabs (Brils et al., 2024)\nwould at least provide a constraint on the shape and range of temperature profiles in ice slab regions.\nThanks for this suggestion. Rather than using modelled temperature profiles for initialisation with all the\nuncertainties in RCM sub-surface schemes, we take advantage of another study that has been published this\nyear: Vandecrux et al. (2024, TC), which mapped observed 10 m depth temperatures ice-sheet-wide using\nmachine learning techniques, presented along with an associated uncertainty product. We use the Vandecrux\ndata to set the bottom temperature of our initial temperature profiles and then scale our observed\nprofile from FS2 accordingly. We also run our calculations +/- the stated 10 m temperature uncertainty.\nPlease see our revised Methods section for more information.\nHow is the long-term equilibrium line in Figure 5 defined?\nThis was based on RACMO data 1981-2000. But in this re-submission we have modified this figure to only\nshow wetness in the Culberg et al. (2024) maximum ice slabs extent area.\nLine 297 – I found the context for the discussion future warming and peak refreezing a bit confusing. Is the\npoint simply that we should not project current SIF rates as a function of SMB into the future due to the\nuncertainty in climate forcing and hypsometry?\nWe have made efforts to improve the narrative of this paragraph – please see full text.\nLine 476 – where did n=4 come from when only 3 PROMICE weather stations were noted in the study area?\nWhich stations were used and how are they distributed in elevation across the study area? (e.g. is it fully\nrepresentative of the elevation distribution, or does this scaling do significant extrapolation outside the data\nbounds?)\nThe middle and upper weather stations largely encompass our study area, with the upper station being situated\nat 1840 m ASL. Above this elevation SEB model input variables are extrapolated, not interpolated, and\ntherefore calculations do not exceed 2000 m ASL so that extrapolation does not become unrealistic. The\nmodel uses these observations to compute in discrete 100 m elevation bins. We computed our regression\nacross four elevation bins (1,600, 1700, 1800 and 1,900 m asl). This was not clear previously, so thanks for\ndrawing this to our attention and we have updated the text accordingly.\nLine 506 – What determines if a cell becomes snow free? Is this from MAR, an albedo remote sensing\nproduct, or something else?\nAccording to modelled MAR snowpack. This is now noted in the text.\nLine 526 – wouldn’t a brief rise above 0C potentially lead to melting that could supply more water for SIF?\nWe have added a note that, given several weeks of prior atmospheric cooling of the snowpack, it is probable\nthat this melt would simply have been refrozen in the uppermost centimetres below the surface.\n\nReviewer #3 (Remarks to the Author):\nSummary\nThis paper uses a combination of field observations, satellite observations, and simple/conceptual modelling\nto address the important topic of meltwater retention in the percolation zone of the Greenland ice sheet –\nabout which there is some uncertainty. The paper builds on the finding that meltwater retention has led to the\nformation of ice slabs at depth the firn pack, and that the presence of these slabs disrupts the normal process of\nmeltwater percolation and redistribution.\nThe paper is quite challenging to comprehend as it contains several elements of work that are not easily linked\nand are each quite obscure. Despite several attempts to read and review, I could not form a clear narrative\nfrom the study and, more importantly, I was not persuaded by the arguments presented. Either the narrative\nneeds to be strengthened, or the arguments need to be made more compelling with further work.\nMajor\nI have identified five major concerns that I hope will act as a guide for the authors to reformulate their paper.\n1. Loose language. There is lots of loose language and (for me) too much hyperbole in this paper. The term\n“runoff” in particular is central to this work but is not well or consistently defined; it is interchangeably used\nto describe meltwater, mobile water, and water that has left the ice sheet. In my book it is the latter only.\nElsewhere the term “rivers” is used but this is informal and not helpful; it would be better to make an\nassessment of the intensity of meltwater flow in supraglacial channels and contrast the magnitude.\nAnd it is stated in the abstract and repeated elsewhere that the amount of runoff from the percolation zone is\nnot known when this is simply not true – the authors used for example regional climate models that make an\nestimate of this and there have been many studies on meltwater retention including the identification of ice\nslabs on which this work is predicated. I would recommend a thorough parse of loose and hyperbolic\nlanguage.\nWe have endeavoured to tighten up our use of language throughout the manuscript, particularly around\nthe term ‘runoff’ which we agree was not used as precisely as it should have been. We hope that these changes\nmeet with approval.\nThe reviewer said that they find ‘too much hyperbole’. The only example they provided is: “it is stated in the\nabstract and repeated elsewhere that the amount of runoff from the percolation zone is not known when this is\nsimply not true – the authors used for example regional climate models that make an estimate of this and there\nhave been many studies on meltwater retention including the identification of ice slabs on which this work is\npredicated.”. Concerning this specific example, we disagree: regional climate models make estimates which\nso far remain unvalidated for the percolation zone, especially where ice slabs are present. They are potentially\nsubject to major issues associated with invalid parameterisations. We currently have a manuscript co-written\nwith the MAR and RACMO teams under discussion which covers many of these themes in detail:\nhttps://doi.org/10.5194/egusphere-2024-2750. Nonetheless, to tighten the language we are happy to change\nfrom ‘unknown’ to ‘unmeasured’, ‘unvalidated’. We also note that we have added further justification\nfor our heat flux modelling approach (over RCMs) in the Methods section ‘Modelling SIF’.\nThe reviewer states that there have been ‘many studies on meltwater retention including the identification of\nice slabs:\n• The reviewer does not provide examples of these ‘many studies’.\n• Taking this comment literally, we think the reviewer is suggesting that there are already several\nstudies of meltwater retention on top of ice slabs. We are not aware that such studies exist, with the\nexceptions of already-cited studies of SIF atop glacier ice in the Canadian Arctic and Svalbard. In\nGreenland, we are only aware of slab mapping/identification exercises, e.g. MacFerrin et al. (2019),\nJullien et al. (2023), Culberg et al. (2024), which do not simultaneously assess meltwater retention\nprocesses.\n• There is arguably perhaps one study, Machguth et al., (2016), which focuses on meltwater retention\non a slab. In this study, the term ‘superimposed ice formation’ appears only once, in the 3rd-last\nparagraph of the Methods, only in reference to core alignment. This study undertook ‘barcode\nalignment’ on a couple of firn cores at KAN_U, which as this reviewer points out subsequently can be\nsubject to substantial uncertainties. But having re-read this study in detail for this response we have\nconcluded that the evidence of SIF they present is so limited that we now do not cite it in the\nIntroduction as an observation of Greenland SIF.\n2. Superimposed ice formation. The initial part of this work is based on in situ measurements that purport to\ndocument superimposed ice formation. I am not persuaded that the authors have demonstrated this effectively.\nThe main approach is to measure the depth to the slab surface two years apart at a variety of locations. This is\nnot an actual measurement of superimposed ice formation, it’s a crude approximation. The slabs could be\nmobile, because the site is not at the equilibrium line. More importantly, the authors are not able to resurvey\nthe same site and so their measurement is in fact a combination of spatial and temporal variation in slab\nsurface depth. Without knowing the former, it is impossible to deduce the latter. This calculation\nunfortunately underpins the entire study, and so it is quite fundamental to convince the reader that the\nmeasurements are robust. Having drilled many firn cores in this region I know that measuring the depth to a\ncommon horizon is a very challenging thing to do with precision, and that there are typically large natural\nvariations in the depth of common horizons in any case. I wonder if repeat GPR would have been a better way\nto do this entire experiment. Or drilling through the entire slab. Or analysing the core chemistry to isolate\nannual layers. Certainly the spatial variation in slab surface depth needs to be reliable quantified somehow\nbefore the temporal variation can be deduced.\nWe strongly disagree with the thrust of this comment that our in-situ observations are invalid. We highlight\nthat neither of the other two reviewers raised concerns of this nature. From our perspective there seems to\nhave been a misunderstanding of our field measurement protocol. At root, this seems to relate to the locations\nat which we make our measurements. We do resurvey the same site between years: these locations are exactly\nthe same each year with respect to the underlying stratigraphy, as we returned to a location marked with a\nbamboo stake. Each location had of course moved by several 10s of metres due to ice flow, but this process\ndoes not affect the near surface stratigraphy. To clarify this further, we have improved our explanation,\nincluding adding a supplementary figure which graphically illustrates the measurement process, and\nsupplementary discussion to discuss local variability in SIF more extensively.\nIn our judgement the suggestion of using GPR is a non-starter. If we were to make a GPR survey, we would\nmeasure the thickness of the snow layer, albeit with larger uncertainty than by drilling. We ask ourselves: how\ndoes GPR measure the thickness of the slab and in particular its annual growth? We could potentially detect\nthe bottom depth of the slab in subsequent years, but would be hampered by all the GPR-related uncertainties\nthat would reduce the precision of the outcome in such a study. And crucially, even if a change in thickness\ncould be measured beyond uncertainty, we still would not know whether the slab grew from accretion below\nor above. In contrast, this is precisely what our approach of returning to a marked surface location in year 2\nenables us to do.\nThe reviewer states “Having drilled so many firn cores in this region I know that measuring the depth to a\ncommon horizon is very challenging…”. We have no way of assessing what expertise this reviewer has, nor\nthe scope of the region they are considering. The reviewer may have drilled firn cores at high elevations, but\nperhaps not into ice slabs. Our cores are from the slabs, so are basically solid ice. In these cores, to say where\nsnow ends and firn starts is not always clear, but the horizon where the slab begins is always abundantly clear.\nFor more details, the core logs can be examined in our data repository.\n3. Spatial variation in SIF. A main conclusion of this work is that the degree of superimposed ice formation is\nhigher closer to surface meltwater channels than it is distal from them. This is I think the most original aspect\nof the study. However, the conclusion itself depends on how reliable the estimates of SIF (see point 2),\nbut\nalso on how strong the spatial correlation is. Looking at the data in Figure 1 I am not convinced that the\nspatial variation in SIF is significantly correlated with the distribution of meltwater channels; there seem to be\njust as many sites that buck the trend as there are that fit the hypothesis. The sites sampled are also quite\nrandomly distributed, which doesn’t help. Wouldn’t it make sense to survey a profile tangential to a channel,\nfor example? I think more work and quite possibly more data are needed to convince the reader that this\nconclusion is robust.\nAs the reviewer notes, “the conclusion itself depends on how reliable the estimates of SIF (see point 2)” are.\nAs we respond above, we are confident in our SIF estimates and believe that there has been a\nmisunderstanding here.\nConcerning spatial variation, Fig. 1 provides a high-level overview of the SIF magnitudes. Further evidence is\npresented in Supplementary Figure 1, which explicitly compares the quantity of SIF with the Euclidean\ndistance to surface drainage visible on optical satellite imagery. Our measurements were concentrated on two\ntransects, with a spacing that reflects the proximity to surface drainage. Sites FS1, 122, camp, 231, 3u3\nconstitute one transect. The sites shown in Fig. 1d constitute another transect, which is quasi-tangential to the\nmain surface channel. We now note these transect configurations in the main text. Concerning a profile\ntruly tangential to a channel: the location of the channel varies somewhat between years, so this is less feasible\nto deploy than it might first appear, and also what a transect equidistant of a stream is one year, would not be\nthe next. Concerning the measurement stations, these were installed to acquire data in three distinctively\ndifferent surface slope regimes and with contrasting surface hydrology in any given melt summer.\nPlease see also our response to reviewer 1 concerning possible spatial variability in melt.\n4. Timing of SIF. There is a section of work on the timing of superimposed ice formation deduced from\nmeasurements of the firn dielectric properties. I found this section very difficult to navigate and not\npersuasive; the only thing that seems to be consistent about the probe measurements is how inconsistent they\nare in time and space. It seems that very few measurements were actually made and yet apparently firm\nconclusions are drawn from them. This is very similar to point 3; it suggests to me that the field experiment\nwas again not up to the task and in hindsight other measurements would have been made. It’s also not clear to\nme how important it is to know the timing of superimposed ice formation – is that really material to this study\n? My feeling is this whole is a distraction and could easily be dropped.\nOn reflection, we agree that knowing the sub-seasonal timing of SIF is not important to our key messages and\nis challenging to interpret from our data. We have therefore removed this whole section as suggested.\n5. Modelling. The study next leaps into a series of “model” experiments designed to help scale up the field\nobservations. However, these experiments are tangentially linked to the field measurements at best. The first\n“model” is a physical 1D heat transfer model, the second is a linear regression between satellite radar\nbackscatter and measured superimposed ice formation, and the third is another linear regression between the\nmodelled superimposed ice formation and cumulative negative degree days derived (I think) from the 1D heat\ntransfer model. The heat flow model is very loosely related to the field observations as it is contrasted to\ntemperature profiles that are not directly related to the SIF. Moreover, we can’t really judge how well the\nmodel and observations agree from the data as presented.\nWe are grateful for these comments, which together with similar comments by reviewer 2 caused us to\nfundamentally rethink how we conduct our upscaling. We believe that we have now successfully\naddressed this concern and that our response to reviewer 2 answers this major comment.\nThe satellite radar backscatter is very loosely related to in situ superimposed ice formation; many other factors\naffect backscatter and in any case the correlation deduced is largely owing to the clustering of the\nmeasurements into two distinct populations, without which there would be no correlation.\nWe recognise that in our previous submission there was the potential to conflate radar backscatter with in-situ\nSIF, primarily because of the statistical way in which we attempted to conduct our upscaling; however, this\nwas not our intention. To be clear, we consider backscatter variations in autumn to tell us only about residual\nwetness in the near-surface. It is the disappearance of this backscatter that we attribute to SIF, which we now\nmodel physically instead of our previous problematic statistical approach. In the revised manuscript this is\nnow clearer because we talk only of residual wetness [duration] and not also in terms of sigmaHV0\nbackscatter strength like previously.\nThe degree of correlation between satellite derived SIF and negative degree days is relatively weak in any\ncase, half of the variation is not explained. I find this approach to be quite chaotic, highly empirical, and not\npersuasive. I would recommend the authors rethink their approach here entirely and find a better way to add\ncontext and value to their sparse field measurements. The authors do note that regional climate models with\ncoupled firn process models are another way to study meltwater retention and I agree this would be a better\nway to approach this problem.\nAs explained in our response summary, we have taken heed of these concerns about the upscaling\napproach and, we believe, now have a much better methodology with which to deliver our upscaled\ncalculations.\nConcerning the suggestion here that RCMs with coupled firn process models would be a ‘better way’ to\napproach the problem: we disagree, and also note that we do not state that they constitute a ‘better way’.\nAs we have tried to clarify in our re-submission, there are substantial uncertainties/problems remaining with how\nRCMs and even observation-based SEB models represent the sub-surface, which makes it desirable to instead\nadopt a simplified approach to modelling SIF. Please see in particular the Methods section ‘Modelling\nsuperimposed ice formation’.\nMinor comments\nLine 1: I’m not sure the word rivers is helpful here as they are generally associated with land and rain fed. can\nyou describe in glaciological terminology instead?\nThe term ‘rivers’ is used elsewhere in the literature. Non exhaustive examples include line 2 of Yang and\nSmith (2016, JGR); Smith et al. (2015, PNAS). We have decided to stick with this intuitive terminology rather\nthan ‘supraglacial streams and rivers’ at every mention.\nLine 4: I don't think its unknown; there are estimates. need to be a bit clearer here what you mean\nWe have changed to ‘unmeasured’.\nLine5: should probably be clear about what you mean by runoff; is it meltwater until it is mobile, or until it is\nknown that it will leave the ice sheet locus, or until it has left the ice sheet locus?\nHere, in-line with the definition of runoff, water that is assumed to leave the ice sheet locus – updated to\nrefer to ‘visible runoff’ which is then linked to previous studies in the Introduction.\nLine 6: I think you mean that up to 3/4 of the annual melt can be retained by SIF (i.e. the converse of what is\nstated), or that the process of SIF can lead to meltwater being retained\nWe are a confused by this comment – the text said that ‘SIF….can retain up to roughly three quarters of\nlocally-produced annual melt’, which doesn’t seem the be converse. Nonetheless, we have tried to clarify in\nour revised Abstract.\nLine 7: Is the heat flow modelling informed somehow by the field observations? or how else are they linked?\nThis comment refers to the abstract. We have added a note that the modelling is initialised by in-situ field\nobservations. The full information is available in the main text as before.\nLine 7: How are the SAR data used? would be helpful to state\nWe have added a note that they are used to derive residual wetness.\nLine 8: Quantify \"majority\"\nWe prefer not to do this in the abstract, as while the numbers are always a majority (i.e. greater than half),\nthe\nprecise value depends on parameter choices. The reader should read the main text for these details.\nLine 8: Quantify \"remainder\"\nSee previous comment.\nLine 9: This probably answers my earlier point about the definition of runoff\nYes, we believe so; see also modified abstract.\nLine 10: Using your model?\nWe are not sure what this comment refers to, sorry.\nLine 11: Does this mean total or local SMB? not sure SMB is a good metric since it is the sum of\naccumulation and ablation\nThis is a good point – we concur that SMB was a poor choice of metric. We now instead express this in\nterms of runoff estimates made using Cryosat-2 altimetry measurements. We estimate that recent SIF\nhas decreased recent runoff by 2-3 %.\nLine 12: Is there an uncertainty on 32 Gt?\nThis value has changed with our new upscaling approach; we now provide a ‘most-likely’ value (60 Gt)\ntogether with a range (27-75 Gt). The range is primarily a function of which ice slab extent we examine, but\nalso varies depending on 10 m firn temperatures and autumn snowpack thickness.\nLine 12: To call them hotspots this is relative to other areas, but it doesn’t seem that you have assessed those?\nConsidering refreezing in absolute terms (metres water equivalent), they are indeed hotspots compared to the\nrest of the ice sheet. Nowhere else is this combination of (a) meltwater supply and (b) a snow/firn matrix in\nwhich to refreeze onto the underlying ice met.\nLine 19: Decreased not densified\nChanged\nLine 25: Is this really true? i think there are estimates in the models you are using- what I guess you are saying\nis that they don't account for potential retention associated with the ice slabs. Thats very different to unknown\nWe have replaced with ‘unmeasured’.\nLine 37: Running water not runoff\nChanged\nLine 40: Radar returns is not a jargon. Measurements?\nWe have changed to “measurements” as suggested, but note that this term has been used in the literature\npreviously (see for example MacFerrin et al. 2019 Supplementary).\nLine 42: Technically the slabs themselves are superimposed ice and so yours are not the first observations\nWe are observing superimposed ice formation, i.e. the process of forming superimposed ice, as we state in the\ntext. Slab mapping has indeed already mapped superimposed ice.\nLine 43: What is the runoff limit?\nIn lines 23-24 we now define the visible runoff limit alongside the visible runoff area.\nLine 44: What is local runoff?\nChanged to ‘net water discharge from the hydrological catchment’.\nFig 1: Why is fs1 high on both years?\nThat there was no slush field at the surface doesn’t mean that there was no slush beneath the surface. Visible\nhydrological features do not tell the whole story, so we also need to look at the sub-surface. C-band SAR\nclearly shows that this area remained wet for a long time after the end of the melt season (Fig. 5a, see wetter\nfeature between FS2 and downstream later, between 1730 and 1760 m asl).\nLine 69: You are showing that the SIF varies by 10s of cm w e over short distances, so perhaps the smb does\ntoo, let alone surface features? In which case you experiment is not well constrained\nPlease see our responses to reviewer 1 concerning surface melt rates.\nLine 88: Why use e>3 when water is e=80?\nThe TDR measurements have been removed from this re-submission.\nLine 94: TDR probe\nThe TDR measurements have been removed from this re-submission.\nFig.3 : It’s very difficult in this figure to contrast the observed and modelled temperature profiles, which seem\nto be the only validation of the model. I suggest replotting them as e.g. a direct scatter plot with a colour\ncoding to illustrate depth so that we can see how well they compare (or otherwise)\nWe have redrawn this figure panel to make the comparisons between observed vs modelled temperature\nprofiles clearer. We have also added some text to the Results which addresses the offset between\nobserved and modelled profiles at the end of the melt season.\nLine 139: Isn't this obviously the case? What else would be causing the SIF. I don't think this is a grand\nconclusion\nWe are not exactly sure what `this’ refers to, as one sentence finishes and another starts at this line number.\nAssuming that ‘this’ refers to the sentence in L137-139, we do not make out that this is a grand conclusion.\nRather, this statement supported the sentence which followed, i.e. ‘In contrast, larger amount of SIF can only\nbe explained by lateral transport of meltwater….’. We have made various minor changes to better explain\nour model findings.\nLine 141: Just because your model agrees with SIF in some places doesn't mean you can conclude that\nsomething else must be happening elsewhere – that’s not a hypothesis you have tested by experiment. Is your\nmodel able to simulate lateral movement of meltwater? What other processes might lead to local variations in\nmelting that are not simulated by your model? More work required here\nThe majority of cores show rather similar SIF. We clearly show that simulated background refreezing agrees\nto SIF in several locations. It clearly does not agree in a few locations. We therefore looked for the reasons for\nthe disagreement. It is obvious that refreezing from only the melting directly above a given location cannot\nvary so enormously in space given how homogeneous the snow-covered surface is and the rather small\nelevation envelope that we worked within (see also our response to reviewer 1). We also observe the largest\nrefreezing rates where there is most slush. To us, it is logical that these locations should therefore have less\ncold content in the slab. Finally, it is unclear to us what the benefit of including lateral movement in our 1-D\nmodel would be.\nLine 144: It does not indicate this; it just indicates you have a match within 1sigma in some locations\nThis line says ‘Since our model indicates that the localised larger quantities of superimposed ice cannot be\naccommodated by slab cold content during the melt season alone…’.\nOur model does indicate this as the modelled summertime SIF is valid for most of our measurement sites, not\njust the one where we compare. Of course, we do not have more firn temperature data with which to undertake\na more thorough comparison. But in the context of how much larger some of SI measurements are (by\ndecimetres), then we consider a match within 1 sigma (0.06 m w.e.) sufficient for our assertions.\nLine 156: Correlation, not relationship\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nFig. 4: These data are clustered in 2 populations which is leading to most of the correlation. I don't think they\nare distributed well enough to develop an empirical parameterisation. why are they clustered in this way? I\nsuggest this is explored in more detail before jumping to conclusions\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nLine 174: This is all a bit contrived; your result is entirely dependent upon the relationship that you force to\nexist.\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nFig. 6: This isn't \"wetness\", its sigma nought. I think you are jumping to too many conclusions here\nAs we describe in the Methods, we derive a binary wetness indication by applying a threshold to detrended\nsigma0HV. Thus, Figure 6 did present wetness (duration), not sigma0. This approach is based on\npreviously published work, most notably Liang et al. (2021). Note that Fig. 6 is now Suppl. Fig. 5.\nLine 211: So only half of the variance is explained\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 219: I'm lost here; if you have a model of the firn pack why are you counting cNDD to model SIF?\nwhy\nnot just compute it from your actual model?\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 334: Such as?\nWe now summarise the specific issues (difficulty locating the slab surface at FS3; uncertainty related to\ndrilling into SIF until metal chips appeared).\nLine 336: Why are they minimum estimates?\nPlease see response to previous comment.\nLine 343: This is a poorly constrained experiment as you don't know the spatial variation in firn stratigraphy.\nhow can you be sure you are measuring superimposed ice?\nWe believe that this concern is addressed by our more detailed explanation of our field measurement\nprotocol and Supplementary Discussion; specifically, during our repeat measurements, we found very little\nspatial variability at metres scales.\nLine 358: So this is an underestimtae of the SIF?\nThis paragraph is about the cSMB, not about SIF. By definition (Cogley et al. 2011), SIF refers only to\nrefreezing on top of impermeable horizons, in this case the ice slab.\nLine 459: Why choose the 1st percentile? is this based on theory?\nPlease see our response to Reviewer 2.\nLine 461: Why likely?\nThanks for pointing this out – on reflection, we do not think we can quantify this statement that the threshold\nis ‘likely conservative’, so have removed it in this re-submission.\nLine 466: Extrapolate or scale up, not distribute\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 470: This isn't determining the existence of a relationship, it is assuming the existence of a relationship\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 486: Detect not retrieve\nChanged.", "reviewer_3": "Tedstone et al, 2024\nSummary\nThis paper uses a combination of field observations, satellite observations, and simple/conceptual modelling to address the important topic of meltwater retention in the percolation zone of the Greenland ice sheet – about which there is some uncertainty. The paper builds on the finding that meltwater retention has led to the formation of ice slabs at depth the firn pack, and that the presence of these slabs disrupts the normal process of meltwater percolation and redistribution.\nThe paper is quite challenging to comprehend as it contains several elements of work that are not easily linked and are each quite obscure. Despite several attempts to read and review, I could not form a clear narrative from the study and, more importantly, I was not persuaded by the arguments presented. Either the narrative needs to be strengthened, or the arguments need to be made more compelling with further work.\nMajor\nI have identified five major concerns that I hope will act as a guide for the authors to reformulate their paper.\n1. Loose language. There is lots of loose language and (for me) too much hyperbole in this paper. The term “runoff” in particular is central to this work but is not well or consistently defined; it is interchangeably used to describe meltwater, mobile water, and water that has left the ice sheet. In my book it is the latter only. Elsewhere the term “rivers” is used but this is informal and not helpful; it would be better to make an assessment of the intensity of meltwater flow in supraglacial channels and contrast the magnitude. And it is stated in the abstract and repeated elsewhere that the amount of runoff from the percolation zone is not known when this is simply not true – the authors used for example regional climate models that make an estimate of this and there have been many studies on meltwater retention including the identification of ice slabs on which this work is predicated. I would recommend a thorough parse of loose and hyperbolic language.\n2. Superimposed ice formation. The initial part of this work is based on in situ measurements that purport to document superimposed ice formation. I am not persuaded that the authors have demonstrated this effectively. The main approach is to measure the depth to the slab surface two years apart at a variety of locations. This is not an actual measurement of superimposed ice formation, it’s a crude approximation. The slabs could be mobile, because the site is not at the equilibrium line. More importantly, the authors are not able to resurvey the same site and so their measurement is in fact a combination of spatial and temporal variation in slab surface depth. Without knowing the former, it is impossible to deduce the latter. This calculation unfortunately underpins the entire study, and so it is quite fundamental to convince the reader that the measurements are robust. Having drilled many firn cores in this region I know that measuring the depth to a common horizon is a very challenging thing to do with precision, and that there are typically large natural variations in the depth of common horizons in any case. I wonder if repeat GPR would have been a better way to do this entire experiment. Or drilling through the entire slab. Or analysing the core chemistry to isolate annual layers. Certainly the spatial variation in slab surface depth needs to be reliable quantified somehow before the temporal variation can be deduced.\n3. Spatial variation in SIF. A main conclusion of this work is that the degree of superimposed ice formation is higher closer to surface meltwater channels than it is distal from them. This is I think the most original aspect of the study. However, the conclusion itself depends on how reliable the estimates of SIF (see point 2), but also on how strong the spatial correlation is. Looking at the data in Figure 1 I am not convinced that the spatial variation in SIF is significantly correlated with the distribution of meltwater channels; there seem to be just as many sites that buck the trend as there are that fit the hypothesis. The sites sampled are also quite randomly distributed, which doesn’t help. Wouldn’t it make sense to survey a profile tangential to a channel, for example? I think more work and quite possibly more data are needed to convince the reader that this conclusion is robust.\n4. Timing of SIF. There is a section of work on the timing of superimposed ice formation deduced from measurements of the firn dielectric properties. I found this section very difficult to navigate and not persuasive; the only thing that seems to be consistent about the probe measurements is how inconsistent they are in time and space. It seems that very few measurements were actually made and yet apparently firm conclusions are drawn from them. This is very similar to point 3; it suggests to me that the field experiment was again not up to the task and in hindsight other measurements would have been made. It’s also not clear to me how important it is to know the timing of superimposed ice formation – is that really material to this study ? My feeling is this whole is a distraction and could easily be dropped.\n5. Modelling. The study next leaps into a series of “model” experiments designed to help scale up the field observations. However, these experiments are tangentially linked to the field measurements at best. The first “model” is a physical 1D heat transfer model, the second is a linear regression between satellite radar backscatter and measured superimposed ice formation, and the third is another linear regression between the modelled superimposed ice formation and cumulative negative degree days derived (I think) from the 1D heat transfer model. The heat flow model is very loosely related to the field observations as it is contrasted to temperature profiles that are not directly related to the SIF. Moreover, we can’t really judge how well the model and observations agree from the data as presented. The satellite radar backscatter is very loosely related to in situ superimposed ice formation; many other factors affect backscatter and in any case the correlation deduced is largely owing to the clustering of the measurements into two distinct populations, without which there would be no correlation. The degree of correlation between satellite derived SIF and negative degree days is relatively weak in any case, half of the variation is not explained. I find this approach to be quite chaotic, highly empirical, and not persuasive. I would recommend the authors rethink their approach here entirely and find a better way to add context and value to their sparse field measurements. The authors do note that regional climate models with coupled firn process models are another way to study meltwater retention and I agree this would be a better way to approach this problem.\nMinor\nLine 1: I’m not sure the word rivers is helpful here as they are generally associated with land and rain fed. can you describe in glaciological terminology instead?\nLine 4: I don't think its unknown; there are estimates. need to be a bit clearer here what you mean\nLine5: should probably be clear about what you mean by runoff; is it meltwater until it is mobile, or until it is known that it will leave the ice sheet locus, or until it has left the ice sheet locus?\nLine 6: I think you mean that up to 3/4 of the annual melt can be retained by SIF (i.e. the converse of what is stated), or that the process of SIF can lead to meltwater being retained\nLine 7: Is the heat flow modelling informed somehow by the field observations? or how else are they linked?\nLine 7: How are the SAR data used? would be helpful to state\nLine 8: Quantify \"majority\"\nLine 8: Quantify \"remainder\"\nLine 9: This probably answers my earlier point about the definition of runoff\nLine 10: Using your model?\nLine 11: Does this mean total or local SMB? not sure SMB is a good metric since it is the sum of accumulation and ablation\nLine 12: Is there an uncertainty on 32 Gt?\nLine 12: To call them hotspots this is relative to other areas, but it doesn’t seem that you have assessed those?\nLine 19: Decreased not densified\nLine 25: Is this really true? i think there are estimates in the models you are using- what I guess you are saying is that they don't account for potential retention associated with the ice slabs. Thats very different to unknown\nLine 37: Running water not runoff\nLine 40: Radar returns is not a jargon. Measurements?\nLine 42: Technically the slabs themselves are superimposed ice and so yours are not the first observations\nLine 43: What is the runoff limit?\nLine 44: What is local runoff?\nFig 1: Why is fs1 high on both years?\nLine 69: You are showing that the SIF varies by 10s of cm w e over short distances, so perhaps the smb does too, let alone surface features? In which case you experiment is not well constrained\nLine 88: Why use e>3 when water is e=80?\nLine 94: TDR probe\nFig.3 : It’s very difficult in this figure to contrast the observed and modelled temperature profiles, which seem to be the only validation of the model. I suggest replotting them as e.g. a direct scatter plot with a colour coding to illustrate depth so that we can see how well they compare (or otherwise)\nLine 139: Isn't this obviously the case? What else would be causing the SIF. I don't think this is a grand conclusion\nLine 141: Just because your model agrees with SIF in some places doesn't mean you can conclude that something else must be happening elsewhere – that’s not a hypothesis you have tested by experiment. Is your model able to simulate lateral movement of meltwater? What other processes might lead to local variations in melting that are not simulated by your model? More work required here\nLine 144: It does not indicate this; it just indicates you have a match within 1sigma in some locations\nLine 156: Correlation, not relationship\nFig. 4: These data are clustered in 2 populations which is leading to most of the correlation. I don't think they are distributed well enough to develop an empirical parameterisation. why are they clustered in this way? I suggest this is explored in more detail before jumping to conclusions\nLine 174: This is all a bit contrived; your result is entirely dependent upon the relationship that you force to exist.\nFig. 6: This isn't \"wetness\", its sigma nought. I think you are jumping to too many conclusions here\nLine 211: So only half of the variance is explained\nLine 219: I'm lost here; if you have a model of the firn pack why are you counting cNDD to model SIF? why not just compute it from your actual model?\nLine 334: Such as?\nLine 336: Why are they minimum estimates?\nLine 343: This is a poorly constrained experiment as you don't know the spatial variation in firn stratigraphy. how can you be sure you are measuring superimposed ice?\nLine 358: So this is an underestimtae of the SIF?\nLine 459: Why choose the 1st percentile? is this based on theory?\nLine 461: Why likely?\nLine 466: Extrapolate or scale up, not distribute\nLine 470: This isn't determining the existence of a relationship, it is assuming the existence of a relationship\nLine 486: Detect not retrieve"}, "author_response": "Response to reviews\nConcurrent superimposed ice formation and meltwater\nrunoff on Greenland’s ice slabs, version 1\nSummary\nWe thank the three reviewers for their engagement with our first submission. We believe that we have\nsuccessfully addressed their major concerns, leading to major improvements which warrant re-submission\nthrough the Nature Communications appeals process.\nReviewer 1 found the manuscript to be well-written, thorough, and provide important observational evidence.\nTheir concerns were essentially minor. The bulk of our changes concern comments by reviewers 2 and 3.\nReviewer 2 said that the manuscript addresses an important caveat to the popular notion that ice slabs lead to\nrunoff. They found the field measurements novel and convincing but were concerned that the remote sensing\nanalysis had poorly quantified but likely large uncertainties. Reviewer 3 had some major concerns about both\nthe field measurements and the upscaling.\nOur most major and important changes are:\n1. We have developed a new upscaling methodology. We no longer rely on statistical links between\nfield-measured SIF, negative degree days and SAR. Instead, we extend our physical heat flux\nmodelling approach to estimate SIF through (a) ‘bottom-up’ (i.e. slab cold content) refreezing during\nboth summer and autumn, and (b) ‘top-down’ refreezing (i.e. negative air temperatures) during\nautumn. We still use SAR sigma0HV, but only to determine wetness duration both across the field\narea and ice-slabs-wide through application of a binary threshold. We undertake these calculations in\nboth our field area and ice-slabs-wide. It also allows a much more rigorous exploration of\nuncertainties. In our view this is a more elegant and physical approach, particularly because it allows\nus to fully differentiate between summer and autumn SIF. Please see our headline response to\nreviewer 2 for more details.\n2. Our ice-slabs-wide analysis now incorporates the new ice slab extents mapped by Culberg et al.\n(2024); as these extents are much larger than Jullien et al. (2023), our estimate of SIF onto ice slabs\nhas increased accordingly.\n3. We have recomputed our SAR wetness threshold and we have examined the sensitivity of our\ncalculations to the choice of the SAR sigma0HV threshold for mapping wetness; our findings are\nrobust across the specified sensitivity range of the Sentinel-1 SAR sensor.\n4. We disagree with the view of reviewer 3 that our in-situ stake measurements of SI are insufficient to\ndemonstrate our purported findings. In response we have extended our description of our SIF\nmeasurement protocol, including the addition of a sketch figure to the supplementary materials. Please\nalso see the rebuttal in this response document. We do however agree that our measurements to\ncapture SIF timing are not fully convincing, so have removed these accordingly.\n5. We have made many improvements to the main text which we hope improve the readability.\nOur changes in bold\nIn our longer responses, the key parts are underlined.\nThis manuscript presents estimates of superimposed ice formation on top of ice slabs along the K-Transect in\nSW Greenland. In this study, the authors present field measurements that suggest local ice formation on ice\nslabs can retain ¾ of local annual surface melt. They estimate that ice formation along ice slabs facilitated the\nretention of ~32 gigatonnes of melt between 2017 and 2022, roughly ~4% of the Greenland SMB.\nOverall, the manuscript is well-written, thorough, and presents important observational evidence for surface\nmelt refreeze onto ice slabs, which plays a role on the overall SMB of the GrIS. I have several comments,\nmost of which are relatively minor; I believe that their consideration will overall improve the manuscript.\nI have organized the comments below, grouped by section within the manuscript.\nResults:\nPg 3 line 77: Because local melt is near uniform across the study area (presumably as a result of their\nmodeling results, but this is not made explicit), they hypothesis that measured spatial variation in SIF is\nindicative of inputs not only from local melt production, but also lateral water flow. This seems entirely\nreasonable given the observed increases in SIF for sites proximal to major surface rivers. However, the\ninclusion of external references to reinforce this claim would be beneficial.\nAs the reviewer infers, we take melt as a smoothly varying field of decreasing melt with increasing elevation,\nas calculated by a surface energy balance model forced by in-situ observations. Fig 5d now illustrates how\nmelt decreases with increasing elevation. We also now refer to findings from Noel et al. 2016 in the first\nparagraph of our Results. They showed that melt rates show very little spatial variability in this region.\nIt is possible that, once visible slush fields have developed, then local melt rates might increase further as a\nfunction of their reduced albedo. However, these slush fields can only initially develop by the transport and\nponding of melt that is essentially uniform in space considering the smooth surface topography over 100s\nmetres scales and homogeneous snow surface cover (e.g. Holmes 1955; Noel et al., 2016; Clerx et al., 2024/in\npress). In other words, their origin does not lie in albedo. Nonetheless, we now note this possibility of\nsubsequent local melt enhancement in the first paragraph of the Discussion.\nPage 4 line 79: “Finally, beyond 2–4 km away from visible surface runoff, we found that SIF began to\nincrease again, as illustrated by 0.25 m w.e. at site 3u3 (Suppl. Fig. 1).” – With what degree of confidence can\nthis claim be made, given that it is reliant upon a single point wise measurement of SIF, and the measured\nvalue is within 2 standard deviations of the reported average SIF from 2022?\nOn reflection, we agree, so we have removed this statement from the main text.\nPage 7 line 139: “In contrast, larger amounts of SIF can only be explained by lateral transport of meltwater\nwhich provides additional water for SIF to continue in specific locations beyond the end of the melt season.” –\nSimilar to first comment; Isn’t it possible that this elevated SIF can also be explained by spatial heterogeneity\nin melt that is not captured in the SMB model? I tend to agree with the authors that the more plausible\nexplanation is lateral transport of meltwater, but to what degree can spatial heterogeneity of melt be ruled out?\nPlease see our response to the previous comment about spatial variability in melt.\nAn important point to mention in the refreezing during Autumn section is the role that volume scattering due\nto surface and subsurface heterogeneities (e.g. surface roughness, sastrugi, ice lenses, ice pipes, percolation\npipes) plays. Implicitly, the bias-removal approach that is employed here *should* account for this process,\nprovided the geometry of these features is assumed temporally non-varying, such that changes in the cross-\npolarization backscatter can be attributed solely to changes in dielectric permittivity rather than surface\nroughness and volume scattering. I would suggest that this point be elaborated upon in the methods section,\n2\nperhaps with reference to some relevant literature (e.g., Rignot 1995, Backscatter Model for the Unusual\nRadar Properties of the Greenland Ice Sheet).\nThanks for this suggestion. We have elaborated in the Methods (‘Detecting residual water presence’)\naccordingly.\nMethods:\nPage 16, line 396: Regarding the model producing 0.25 m w.e. melt; is this merely a consequence of the\nsubsurface model discretization being 0.25 m? If not, what is the uncertainty of this model-derived estimate of\nw.e. melt?\nNo, this is not a consequence of the sub-surface model discretisation. In each time step, the model’s grid\nmoves up/down after the occurrence of accumulation/ablation; van As et al. (2017) contains full details.\nConcerning uncertainty: the model does not come with uncertainty bounds because of the complex and non-\nlinear nature of the modelling. This makes validation a crucial aspect of SEB modelling, as is done in\nnumerous previous studies. Concerning this model specifically, see e.g. van As et al. (2017) cited in\nmanuscript. In this manuscript we therefore evaluate how sensible its melt estimates are by reference to in-situ\nsurface lowering measurements, which show excellent agreement (see Methods).\nPage 18 line 469: Section on the linear regression between measured SIF and mean crosspol backscatter:\ngiven that an uncertainty in the in situ SIF values are provided, it would be worthwhile to propagate this\nuncertainty through the linear regression for a nominal uncertainty estimate of the linearly-regressed area-\nwide SIF estimates, particularly since the R2 in the linear regression isn’t particularly strong.\nThe linear regression approach has been removed in this re-submission.\nThis manuscript describes field measurements of superimposed ice formation on top of ice slabs on\nGreenland’s K-transect in the summers of 2021 and 2022. Using repeated firn cores, time domain\nreflectometry, and thermal modeling, the authors show that a significant portion of annual melt refreezes on\ntop of the ice slabs due to ice slab cold content and cold surface temperatures during autumn. In particular, the\nauthors show that refreezing occurs even where water flows laterally over the surface, potentially for kms,\nsuggesting that surface mass loss may not scale directly with visible surface runoff during the melt season.\nThe authors then use Sentinel-1 SAR measurements to establish relationships between their field\nmeasurements of SIF and radar-inferred surface wetness measurements. They use these relationships to scale\ntheir SIF measurements first to the area surrounding their field site and then to the entire ice sheet. The\nauthors conclude that SIF accounts for about 4% of the contemporary GrIS SMB, making it a small but non-\ntrivial contributor to the overall ice sheet mass balance.\nOverall, this paper addresses an important caveat to the popular notion that ice slabs lead to complete surface\nrunoff. In particular, the field measurements are novel and convincing and suggest that superimposed ice\nformation (SIF) needs to be accounted for in order to understand how quickly surface melting becomes mass\nloss as melting advances inland on the ice sheet. However, the remote sensing analyses used to scale these\nfield measurements to the entire ice sheet are subject to large uncertainties, many of which are not well\nquantified in the current manuscript. I am particularly concerned that the approach to estimating pan-\nGreenland SIF relies on a highly uncertain relationship between SIF and nCDD in a single ~30 km x 30 km\nregion which cannot be validated anywhere else on the ice sheet, as well a background SIF rate that is\nassumed to not vary at all across the ice sheet. This opens the potential for large, uncharacterizable regional\nbiases to affect the pan-Greenland SIF estimates.\n‘[the authors use a] background SIF rate that is assumed to not vary at all across the ice sheet’. In fact, the\nbackground SIF rate used for the ice-slabs-wide analysis varies in two ways:\n(1) As in our first submission, our treatment of the conditions which permit SIF across the ice sheet is\nexplicitly spatial. At a given MAR model grid node, SIF can only occur during the melt season and\nwhen the slab is covered in (porous) snow which provides the basis for a slush matrix overlying the\nslab.\n(2) Previously, we used three different fixed initial temperature scenarios to model bottom-up summer\nSIF. In this re-submission, we instead use the state-of-the-art ice-sheet-wide 10 m depth\ntemperatures data product published by Vandecrux et al. (2024, TC) to initialise our bottom-up\nSIF modelling. We also run our model at +/- the stated uncertainty of the Vandecrux product.\nWe consider this to be a major improvement compared to the previous version of our manuscript.\nConcerning the remote sensing analyses: we concur with the reviewer’s concerns that these were subject to\nlarge uncertainties and had the potential for large, uncharacterizable regional biases to affect the pan-GrIS SIF\nestimates. In conjunction with similar concerns raised by reviewer 3, we have therefore made major\nrevisions which replace the previous statistical approach completely with physical heat flow modelling\nto make our upscaling calculations. In summary:\n(1) In areas which remain wet after the end of the melt season according to SAR, we now allow our\nsummertime “bottom-up” calculations to continue into autumn until SAR indicates that the area is\ndry, yielding the additional bottom-up SIF which occurs in autumn. We partition the contributions of\nsummer vs autumn bottom-up SIF using the last day of melting from passive microwave radiometry.\n(2) In these same areas, we now also explicitly address “top-down” refreezing as a result of negative\nsurface temperatures, beginning once the melt season has finished. These calculations are very\nsensitive to the thickness of the snowpack overlying the slush matrix; we discuss this further in the\nmanuscript.\n(3) Any pixels suspected of being lakes (i.e. still ‘residually wet’ at the end of the year) are now\ncompletely excluded from our calculations.\n(4) We have recalculated the SAR sensitivity threshold using data acquired in all areas labelled as\n‘minimum likely’ slabs by Culberg et al. (2024, TC), with little resulting change (see below).\n4\n(5) Over our field area, we now show the uncertainty in our SIF calculations associated with the stated\nsensitivity of the Sentinel-1 SAR sensor used to determine wetness duration. We show that this source\nof uncertainty is second order, i.e. wetness duration does not change appreciably within this\nsensitivity range. We also checked this on an ice-slabs-wide basis where again it is a second-order\nsource of uncertainty (see Supplementary Figure 4).\nMajor Comments:\n[1] I am a bit confused about the primary purpose of modelling of the rate of superimposed ice formation. The\nauthors state that the reason they conduct the modelling is to determine the degree to which ice slab warming\nis due to latent heating. If this is indeed the main goal, then the argument would be better supported by\ncomparing the model run presented in the paper to a model run forced by the measured surface temperature\nand which only considers heat conduction (no latent heating). The difference between those two scenarios\nwould provide a stronger quantification of the role of latent heating in warming the ice slabs.\nHowever, the primary conclusion that seems to come out of this section is that refreezing due\nto ice slab cold content can only account for measured SIF rates in regions without additional\nadvection of meltwater. If the goal is instead to accurately model local SIF, then it seems like the authors\nmight be better served by using their fully coupled SEB-firn model instead of this simplified model that\nimposes a fixed meltwater supply that may or may not be consistent with the actual melt supply and timing\nduring the season when measurements were made.\nWe agree that we should have been clearer here. Our main goal was neither to super-accurately model local\nSIF nor to determine the degree to which ice slab warming is due to latent heating. Rather, it was to provide a\nconstraint on how much SIF can occur during the melt season under the plausible assumption of a continuous\nmeltwater supply at the slab surface. Then, SIF in excess of this logically occurs during autumn, i.e. after\nactive surface melt has finished. In other words, the primary purpose of this modelling approach is to partition\nthe SIF that we measured into the melt and post-melt seasons, thereby setting the scene for a slabs-wide first-\norder estimate of SIF later in the study. With our revised modelling approach we believe that this is\nmotivation is much clearer.\nAn ideal solution would indeed be to use a fully coupled SEB-firn model. However, no such model that is\nsuitable for our application exists. We are dealing with the problem of ponded meltwater in a snow matrix\natop an impermeable and cold barrier, which is an extreme end-member of what can be considered ‘firn’. We\nunderstand from firn modelling colleagues that even when ponding subroutines are available then they often\nremain disabled for want of validation. We now describe these issues in the Methods section ‘Modelling\nSIF’. Thus, it is more informative to simplify the problem as far as possible, i.e. by considering a relatively\nthin snowpack with a continuously wet slush matrix on top of an ice slab (conditions which were met for\nmuch of the melt season according to our in-situ measurements and further bolstered by observations in Clerx\net al., 2022) and then using classical heat flux theory to determine the rate of refreezing at first order.\nNew in this re-submission is an alternative configuration of the model to estimate SIF caused by\natmospheric cooling (‘top-down’). As we describe in the paper, an issue here is that we lack good bounds on\nthe thickness of the snowpack overlying the slush matrix. We wanted to use end-of-season snowpack\nthickness simulated by MAR but found that the outputs were highly spatially incoherent, which is a further\ncaution against 'blind’ use of RCMs assuming that they always deliver the best results. The reviewer may be\ninterested to see a paper with several of the same co-authors as this submission currently in Discussion, which\ndescribes some of the difficulties with RCM behaviour at the runoff limits (https://doi.org/10.5194/egusphere-\n2024-2750).\nOverall, we contend that our chosen approach is much simpler to interpret than the outputs of a fully coupled\nSEB-firn model. In turn, we hope that our findings can be used to advance ponding schemes in firn models.\nWe have made narrative changes to the relevant parts of the manuscript accordingly.\n5\n[2] Uncertainty in the regression between measured SIF and the multi-temporal mean 𝜎̂ 𝐻𝑉 can be better\nquantified and mapped into uncertainty in the study site SIF map (Fig. 4b). One option would be to use a\nBayesian regression that can quantify the range of linear models that are consistent with the data. This would\nmake it possible to produce an uncertainty map showing the range of the plausible SIF for each pixel in Fig.\n4b given the range of plausible models. But the key point is that there should be corresponding uncertainty\nmap for Fig 4b that account for uncertainty in the SIF-𝜎̂ 𝐻𝑉 relationship and the uncertainty in the background\nSIF rate.\nSimilarly, better uncertainty quantification is needed for the SIF-nCDD regression. Particularly given the\nrelatively poor R2 value for this regression, I strongly feel that it is important to quantify the uncertainty in the\nmodel fit and clearly demonstrate how that uncertainty maps onto uncertainty in the final estimates of regional\nautumn SIF that are presented in Figure 7.\nThis analysis has been removed from our re-submission.\n[3] For setting the SAR wetness backscatter threshold, I would like to see some sensitivity tests from different\nregions of the ice sheet showing how the 1st percentile of 𝜎̂ 𝐻𝑉 prior to the melt season varies regionally. For\na reliable pan-Greenland mapping, it is very important to establish whether conditions in SW Greenland are\nconsistent across the ice sheet. Other studies that have used Sentinel-1 data to quantify either melting or the\nextent of wet snow facies over Greenland have found that performance with a single threshold can vary\nsignificantly across the ice sheet (Hu et al., 2022; Li et al., 2023). Regional bias due to a fixed wetness\nthreshold should be quantified and incorporated into assessments of overall uncertainty.\nWe thank the reviewer for raising these concerns. First, we comment on the literature citations:\n(1) The findings of Hu et al. 2022 are not as relevant to our study as they first appear. Hu et al. did not\nnormalise summertime retrievals with respect to winter, unlike in our current approach, and so\ntherefore a considerable source of variability remains.\n(2) Li et al. 2024’s methodology is more similar to ours. They emphasise that SAR-based mapping of\nwetness at a set threshold is only reliable above the snowline, but that there is no appreciable regional\nvariability above the snowline. As we consider only snow-covered surfaces above the ELA, the\ndifficulties that Li et al. experienced with selecting a threshold applicable to melt detection over bare\nice is not relevant here.\nWe recap that our threshold of -2.1 dB was quite similar to Liang et al. 2021, who rather than the 1st percentile\ntake the less conservative 5th percentile of dry snow as their threshold: -2.66 dB. Both Liang et al., 2021 and\nLi et al. 2024 employ a more conservative value of -3 dB, +/- 0.3 dB, but they had the express aim of\ncapturing surface melting conditions which they validated by reference to in-situ air temperature\nmeasurements. This is not a relevant validation metric in our case as we are seeking to identify residual\nwetness after the end of active melt, which also indicates that a 3 dB threshold would be too conservative for\nour purposes here.\nNotwithstanding this comparison, we have undertaken another examination of our threshold choice. Now, we\ncompute histograms by sampling all areas within the ‘minimum likely ice slab’ classification of Culberg\net al. (2024), rather than just the western flank of the ice sheet as previously. The 1st percentile value\nhas decreased slightly to -1.17 dB (i.e. 0.07 dB more conservative than in our previous submission).\nIn the revised manuscript we adopt the new threshold of -1.17 dB, +- 0.3 dB indicative of the\nradiometric uncertainty of Sentinel-1’s sensor as per Geudtner and Torres (2012). We use the +/- to\npropagate an uncertainty in wetness duration through to our final calculations in our field area. On an\nice-slabs-wide basis, we provide a new Supplementary Figure which shows that the residual wetness\nduration/area is essentially insensitive to the radiometric uncertainty of the sensor.\n[4] I am concerned about how well the nCDD model can be generalized to map autumn SIF rates across the\nentire ice sheet. This is conceptually quite similar to the positive degree day model used for melting, which\nhas a spatially and temporally variable rate factor that typically needs to be tuned for individual sites. I have\nsimilar concerns for the nCDD model. Essentially, the authors’ scaling assumes that air temperature is the\nonly driver of refreezing and neglects all the other components of the subsurface energy balance, to include\nfirn and snow conditions, which might vary significantly across the ice sheet. It also assumes that the amount\nof melt that is refrozen scales directly with the length of time required to fully refreeze an area, which may be\ngenerally reasonable, but would absolutely be influenced by the transient accumulation rate, subsurface\ntemperatures, etc on the local scale. Without secondary validation sites, it seems basically impossible to know\nif the regression model coefficients will apply anywhere outside of southwest Greenland.\nThis analysis has been removed from our re-submission.\n[5] How are buried lakes accounted for when calculating autumn SIF? If I have understood the methods\ncorrectly, it seems like there has been no attempt to separate that mode of water\nimpoundment from superimposed ice formation. It seems problematic to apply a relationship\nbetween nCDD and refreezing in porous snow or a slush matrix to the refreezing of large open water bodies.\nSince buried lakes will have the highest and most persistent wetness values (see the authors’ own Figure 5,\nwhich is dominated by lakes), I am concerned that including lakes may significantly bias the results.\nThanks for raising this concern. Now, we screen out likely lakes by removing pixels which are still wet in\nthe last SAR acquisition each calendar year.\nMinor Comments:\nOverall, the section titles and introductory sentences could do a better job of explaining to the reader\nthe end goal of the analyses being described in the section. For example, the section “Modelling the rate of\nsuperimposed ice formation” would be easier to follow if it started with a sentence explaining why the authors\nwant to model this quantity. The section title could also better reflect the findings instead of the methods, for\nexample, something like “Ice slab cold content a key factor in summer SIF”.\nThanks for this suggestion. We have made several changes accordingly.\nLine 34 – I am not sure I understand what the second half of this sentence is trying to express.\nWe have rephrased as we agree that this was very unclear previously.\nLine 92 – Why was there no spike in 𝜖 when this water reached and fully flooded the probe before refreezing?\nWe cannot know for sure, but as the section ‘Timing of SIF’ has been removed in response to reviewer 3\nconcerns, we have not investigated further.\nLine 185/296 – there are now estimates of ice slab extent from both modeling (Brils et al., 2024) and Sentinel-\n1 observations (Culberg et al., 2024) that provide a more up-to-date and comprehensive outline of ice slab\nextent on the GrIS.\nIndeed - though we find it a little concerning that literature published in 2023 can already be considered less\n‘up-to-date’! Nonetheless, given that we have made major changes to our upscaling approach, we have taken\nthe opportunity to replace the Jullien et al. (2023) outlines with those based on Sentinel-1 observations\nby Culberg et al. (2024). We now explore the impact that the uncertainty in ice slab extent makes on our final\nestimates of SIF.\nWe have moved the statements of ice slab extent previously found at this location into the Introduction.\nwould at least provide a constraint on the shape and range of temperature profiles in ice slab regions.\nThanks for this suggestion. Rather than using modelled temperature profiles for initialisation with all the\nuncertainties in RCM sub-surface schemes, we take advantage of another study that has been published this\nyear: Vandecrux et al. (2024, TC), which mapped observed 10 m depth temperatures ice-sheet-wide using\nmachine learning techniques, presented along with an associated uncertainty product. We use the Vandecrux\ndata to set the bottom temperature of our initial temperature profiles and then scale our observed\nprofile from FS2 accordingly. We also run our calculations +/- the stated 10 m temperature uncertainty.\nPlease see our revised Methods section for more information.\nHow is the long-term equilibrium line in Figure 5 defined?\nThis was based on RACMO data 1981-2000. But in this re-submission we have modified this figure to only\nshow wetness in the Culberg et al. (2024) maximum ice slabs extent area.\nLine 297 – I found the context for the discussion future warming and peak refreezing a bit confusing. Is the\npoint simply that we should not project current SIF rates as a function of SMB into the future due to the\nuncertainty in climate forcing and hypsometry?\nWe have made efforts to improve the narrative of this paragraph – please see full text.\nLine 476 – where did n=4 come from when only 3 PROMICE weather stations were noted in the study area?\nWhich stations were used and how are they distributed in elevation across the study area? (e.g. is it fully\nrepresentative of the elevation distribution, or does this scaling do significant extrapolation outside the data\nbounds?)\nThe middle and upper weather stations largely encompass our study area, with the upper station being situated\nat 1840 m ASL. Above this elevation SEB model input variables are extrapolated, not interpolated, and\ntherefore calculations do not exceed 2000 m ASL so that extrapolation does not become unrealistic. The\nmodel uses these observations to compute in discrete 100 m elevation bins. We computed our regression\nacross four elevation bins (1,600, 1700, 1800 and 1,900 m asl). This was not clear previously, so thanks for\ndrawing this to our attention and we have updated the text accordingly.\nLine 506 – What determines if a cell becomes snow free? Is this from MAR, an albedo remote sensing\nproduct, or something else?\nAccording to modelled MAR snowpack. This is now noted in the text.\nLine 526 – wouldn’t a brief rise above 0C potentially lead to melting that could supply more water for SIF?\nWe have added a note that, given several weeks of prior atmospheric cooling of the snowpack, it is probable\nthat this melt would simply have been refrozen in the uppermost centimetres below the surface.\nSummary\nThis paper uses a combination of field observations, satellite observations, and simple/conceptual modelling\nto address the important topic of meltwater retention in the percolation zone of the Greenland ice sheet –\nabout which there is some uncertainty. The paper builds on the finding that meltwater retention has led to the\nformation of ice slabs at depth the firn pack, and that the presence of these slabs disrupts the normal process of\nmeltwater percolation and redistribution.\nThe paper is quite challenging to comprehend as it contains several elements of work that are not easily linked\nand are each quite obscure. Despite several attempts to read and review, I could not form a clear narrative\nfrom the study and, more importantly, I was not persuaded by the arguments presented. Either the narrative\nneeds to be strengthened, or the arguments need to be made more compelling with further work.\nMajor\nI have identified five major concerns that I hope will act as a guide for the authors to reformulate their paper.\n1. Loose language. There is lots of loose language and (for me) too much hyperbole in this paper. The term\n“runoff” in particular is central to this work but is not well or consistently defined; it is interchangeably used\nto describe meltwater, mobile water, and water that has left the ice sheet. In my book it is the latter only.\nElsewhere the term “rivers” is used but this is informal and not helpful; it would be better to make an\nassessment of the intensity of meltwater flow in supraglacial channels and contrast the magnitude.\nAnd it is stated in the abstract and repeated elsewhere that the amount of runoff from the percolation zone is\nnot known when this is simply not true – the authors used for example regional climate models that make an\nestimate of this and there have been many studies on meltwater retention including the identification of ice\nslabs on which this work is predicated. I would recommend a thorough parse of loose and hyperbolic\nlanguage.\nWe have endeavoured to tighten up our use of language throughout the manuscript, particularly around\nthe term ‘runoff’ which we agree was not used as precisely as it should have been. We hope that these changes\nmeet with approval.\nThe reviewer said that they find ‘too much hyperbole’. The only example they provided is: “it is stated in the\nabstract and repeated elsewhere that the amount of runoff from the percolation zone is not known when this is\nsimply not true – the authors used for example regional climate models that make an estimate of this and there\nhave been many studies on meltwater retention including the identification of ice slabs on which this work is\npredicated.”. Concerning this specific example, we disagree: regional climate models make estimates which\nso far remain unvalidated for the percolation zone, especially where ice slabs are present. They are potentially\nsubject to major issues associated with invalid parameterisations. We currently have a manuscript co-written\nwith the MAR and RACMO teams under discussion which covers many of these themes in detail:\nhttps://doi.org/10.5194/egusphere-2024-2750. Nonetheless, to tighten the language we are happy to change\nfrom ‘unknown’ to ‘unmeasured’, ‘unvalidated’. We also note that we have added further justification\nfor our heat flux modelling approach (over RCMs) in the Methods section ‘Modelling SIF’.\nThe reviewer states that there have been ‘many studies on meltwater retention including the identification of\nice slabs:\n• The reviewer does not provide examples of these ‘many studies’.\n• Taking this comment literally, we think the reviewer is suggesting that there are already several\nstudies of meltwater retention on top of ice slabs. We are not aware that such studies exist, with the\nexceptions of already-cited studies of SIF atop glacier ice in the Canadian Arctic and Svalbard. In\nJullien et al. (2023), Culberg et al. (2024), which do not simultaneously assess meltwater retention\nprocesses.\n• There is arguably perhaps one study, Machguth et al., (2016), which focuses on meltwater retention\non a slab. In this study, the term ‘superimposed ice formation’ appears only once, in the 3rd-last\nparagraph of the Methods, only in reference to core alignment. This study undertook ‘barcode\nalignment’ on a couple of firn cores at KAN_U, which as this reviewer points out subsequently can be\nsubject to substantial uncertainties. But having re-read this study in detail for this response we have\nconcluded that the evidence of SIF they present is so limited that we now do not cite it in the\nIntroduction as an observation of Greenland SIF.\n2. Superimposed ice formation. The initial part of this work is based on in situ measurements that purport to\ndocument superimposed ice formation. I am not persuaded that the authors have demonstrated this effectively.\nThe main approach is to measure the depth to the slab surface two years apart at a variety of locations. This is\nnot an actual measurement of superimposed ice formation, it’s a crude approximation. The slabs could be\nmobile, because the site is not at the equilibrium line. More importantly, the authors are not able to resurvey\nthe same site and so their measurement is in fact a combination of spatial and temporal variation in slab\nsurface depth. Without knowing the former, it is impossible to deduce the latter. This calculation\nunfortunately underpins the entire study, and so it is quite fundamental to convince the reader that the\nmeasurements are robust. Having drilled many firn cores in this region I know that measuring the depth to a\ncommon horizon is a very challenging thing to do with precision, and that there are typically large natural\nvariations in the depth of common horizons in any case. I wonder if repeat GPR would have been a better way\nto do this entire experiment. Or drilling through the entire slab. Or analysing the core chemistry to isolate\nannual layers. Certainly the spatial variation in slab surface depth needs to be reliable quantified somehow\nbefore the temporal variation can be deduced.\nWe strongly disagree with the thrust of this comment that our in-situ observations are invalid. We highlight\nthat neither of the other two reviewers raised concerns of this nature. From our perspective there seems to\nhave been a misunderstanding of our field measurement protocol. At root, this seems to relate to the locations\nat which we make our measurements. We do resurvey the same site between years: these locations are exactly\nthe same each year with respect to the underlying stratigraphy, as we returned to a location marked with a\nbamboo stake. Each location had of course moved by several 10s of metres due to ice flow, but this process\ndoes not affect the near surface stratigraphy. To clarify this further, we have improved our explanation,\nincluding adding a supplementary figure which graphically illustrates the measurement process, and\nsupplementary discussion to discuss local variability in SIF more extensively.\nIn our judgement the suggestion of using GPR is a non-starter. If we were to make a GPR survey, we would\nmeasure the thickness of the snow layer, albeit with larger uncertainty than by drilling. We ask ourselves: how\ndoes GPR measure the thickness of the slab and in particular its annual growth? We could potentially detect\nthe bottom depth of the slab in subsequent years, but would be hampered by all the GPR-related uncertainties\nthat would reduce the precision of the outcome in such a study. And crucially, even if a change in thickness\ncould be measured beyond uncertainty, we still would not know whether the slab grew from accretion below\nor above. In contrast, this is precisely what our approach of returning to a marked surface location in year 2\nenables us to do.\nThe reviewer states “Having drilled so many firn cores in this region I know that measuring the depth to a\ncommon horizon is very challenging…”. We have no way of assessing what expertise this reviewer has, nor\nthe scope of the region they are considering. The reviewer may have drilled firn cores at high elevations, but\nperhaps not into ice slabs. Our cores are from the slabs, so are basically solid ice. In these cores, to say where\nsnow ends and firn starts is not always clear, but the horizon where the slab begins is always abundantly clear.\nFor more details, the core logs can be examined in our data repository.\n3. Spatial variation in SIF. A main conclusion of this work is that the degree of superimposed ice formation is\nhigher closer to surface meltwater channels than it is distal from them. This is I think the most original aspect\nof the study. However, the conclusion itself depends on how reliable the estimates of SIF (see point 2), but\nalso on how strong the spatial correlation is. Looking at the data in Figure 1 I am not convinced that the\nspatial variation in SIF is significantly correlated with the distribution of meltwater channels; there seem to be\njust as many sites that buck the trend as there are that fit the hypothesis. The sites sampled are also quite\nrandomly distributed, which doesn’t help. Wouldn’t it make sense to survey a profile tangential to a channel,\nfor example? I think more work and quite possibly more data are needed to convince the reader that this\nconclusion is robust.\nAs the reviewer notes, “the conclusion itself depends on how reliable the estimates of SIF (see point 2)” are.\nAs we respond above, we are confident in our SIF estimates and believe that there has been a\nmisunderstanding here.\nConcerning spatial variation, Fig. 1 provides a high-level overview of the SIF magnitudes. Further evidence is\npresented in Supplementary Figure 1, which explicitly compares the quantity of SIF with the Euclidean\ndistance to surface drainage visible on optical satellite imagery. Our measurements were concentrated on two\ntransects, with a spacing that reflects the proximity to surface drainage. Sites FS1, 122, camp, 231, 3u3\nconstitute one transect. The sites shown in Fig. 1d constitute another transect, which is quasi-tangential to the\nmain surface channel. We now note these transect configurations in the main text. Concerning a profile\ntruly tangential to a channel: the location of the channel varies somewhat between years, so this is less feasible\nto deploy than it might first appear, and also what a transect equidistant of a stream is one year, would not be\nthe next. Concerning the measurement stations, these were installed to acquire data in three distinctively\ndifferent surface slope regimes and with contrasting surface hydrology in any given melt summer.\nPlease see also our response to reviewer 1 concerning possible spatial variability in melt.\n4. Timing of SIF. There is a section of work on the timing of superimposed ice formation deduced from measurements of the\nfirn dielectric properties. I found this section very difficult to navigate and not persuasive; the only thing that\nseems to be consistent about the probe measurements is how inconsistent they are in time and space. It seems\nthat very few measurements were actually made and yet apparently firm conclusions are drawn from them.\nThis is very similar to point 3; it suggests to me that the field experiment was again not up to the task and\nin hindsight other measurements would have been made. It’s also not clear to me how important it is to know\nthe timing of superimposed ice formation – is that really material to this study ? My feeling is this whole is a distraction and could easily be dropped.\nOn reflection, we agree that knowing the sub-seasonal timing of SIF is not important to our key messages and\nis challenging to interpret from our data. We have therefore removed this whole section as suggested.\n5. Modelling. The study next leaps into a series of “model” experiments designed to help scale up the field\nobservations. However, these experiments are tangentially linked to the field measurements at best. The first\n“model” is a physical 1D heat transfer model, the second is a linear regression between satellite radar\nbackscatter and measured superimposed ice formation, and the third is another linear regression between the\nmodelled superimposed ice formation and cumulative negative degree days derived (I think) from the 1D heat\ntransfer model. The heat flow model is very loosely related to the field observations as it is contrasted to\ntemperature profiles that are not directly related to the SIF. Moreover, we can’t really judge how well the\nmodel and observations agree from the data as presented.\nWe are grateful for these comments, which together with similar comments by reviewer 2 caused us to\nfundamentally rethink how we conduct our upscaling. We believe that we have now successfully\naddressed this concern and that our response to reviewer 2 answers this major comment.\nThe satellite radar backscatter is very loosely related to in situ superimposed ice formation; many other factors\naffect backscatter and in any case the correlation deduced is largely owing to the clustering of the\nmeasurements into two distinct populations, without which there would be no correlation.\nWe recognise that in our previous submission there was the potential to conflate radar backscatter with in-situ\nSIF, primarily because of the statistical way in which we attempted to conduct our upscaling; however, this\nwas not our intention. To be clear, we consider backscatter variations in autumn to tell us only about residual\nwetness in the near-surface. It is the disappearance of this backscatter that we attribute to SIF, which we now\nmodel physically instead of our previous problematic statistical approach. In the revised manuscript this is\nnow clearer because we talk only of residual wetness [duration] and not also in terms of sigmaHV0\nbackscatter strength like previously.\nThe degree of correlation between satellite derived SIF and negative degree days is relatively weak in any\ncase, half of the variation is not explained. I find this approach to be quite chaotic, highly empirical, and not\npersuasive. I would recommend the authors rethink their approach here entirely and find a better way to add\ncontext and value to their sparse field measurements. The authors do note that regional climate models with\ncoupled firn process models are another way to study meltwater retention and I agree this would be a better\nway to approach this problem.\nAs explained in our response summary, we have taken heed of these concerns about the upscaling\napproach and, we believe, now have a much better methodology with which to deliver our upscaled\ncalculations.\nConcerning the suggestion here that RCMs with coupled firn process models would be a ‘better way’ to\napproach the problem: we disagree, and also note that we do not state that they constitute a ‘better way’. As\nwe have tried to clarify in our re-submission, there are substantial uncertainties/problems remaining with how\nRCMs and even observation-based SEB models represent the sub-surface, which makes it desirable to instead\nadopt a simplified approach to modelling SIF. Please see in particular the Methods section ‘Modelling\nSIF’.\nMinor comments\nLine 1: I’m not sure the word rivers is helpful here as they are generally associated with land and rain fed. can\nyou describe in glaciological terminology instead?\nThe term ‘rivers’ is used elsewhere in the literature. Non exhaustive examples include line 2 of Yang and\nSmith (2016, JGR); Smith et al. (2015, PNAS). We have decided to stick with this intuitive terminology rather\nthan ‘supraglacial streams and rivers’ at every mention.\nLine 4: I don't think its unknown; there are estimates. need to be a bit clearer here what you mean\nWe have changed to ‘unmeasured’.\nLine5: should probably be clear about what you mean by runoff; is it meltwater until it is mobile, or until it is\nknown that it will leave the ice sheet locus, or until it has left the ice sheet locus?\nHere, in-line with the definition of runoff, water that is assumed to leave the ice sheet locus – updated to\nrefer to ‘visible runoff’ which is then linked to previous studies in the Introduction.\nLine 6: I think you mean that up to 3/4 of the annual melt can be retained by SIF (i.e. the converse of what is\nstated), or that the process of SIF can lead to meltwater being retained\nWe are a confused by this comment – the text said that ‘SIF….can retain up to roughly three quarters of\nlocally-produced annual melt’, which doesn’t seem the be converse. Nonetheless, we have tried to clarify in\nour revised Abstract.\nLine 7: Is the heat flow modelling informed somehow by the field observations? or how else are they linked?\nThis comment refers to the abstract. We have added a note that the modelling is initialised by in-situ field\nobservations. The full information is available in the main text as before.\nWe have added a note that they are used to derive residual wetness.\nLine 8: Quantify \"majority\"\nWe prefer not to do this in the abstract, as while the numbers are always a majority (i.e. greater than half),\nthe\nprecise value depends on parameter choices. The reader should read the main text for these details.\nLine 8: Quantify \"remainder\"\nSee previous comment.\nLine 9: This probably answers my earlier point about the definition of runoff\nYes, we believe so; see also modified abstract.\nLine 10: Using your model?\nWe are not sure what this comment refers to, sorry.\nLine 11: Does this mean total or local SMB? not sure SMB is a good metric since it is the sum of\naccumulation and ablation\nThis is a good point – we concur that SMB was a poor choice of metric. We now instead express this in\nterms of runoff estimates made using Cryosat-2 altimetry measurements. We estimate that recent SIF\nhas decreased recent runoff by 2-3 %.\nLine 12: Is there an uncertainty on 32 Gt?\nThis value has changed with our new upscaling approach; we now provide a ‘most-likely’ value (60 Gt)\ntogether with a range (27-75 Gt). The range is primarily a function of which ice slab extent we examine,\nbut\nalso varies depending on 10 m firn temperatures and autumn snowpack thickness.\nLine 12: To call them hotspots this is relative to other areas, but it doesn’t seem that you have assessed those?\nConsidering refreezing in absolute terms (metres water equivalent), they are indeed hotspots compared to the\nrest of the ice sheet. Nowhere else is this combination of (a) meltwater supply and (b) a snow/firn matrix in\nwhich to refreeze onto the underlying ice met.\nLine 19: Decreased not densified\nChanged\nLine 25: Is this really true? i think there are estimates in the models you are using- what I guess you are saying\nis that they don't account for potential retention associated with the ice slabs. Thats very different to unknown\nWe have replaced with ‘unmeasured’.\nLine 37: Running water not runoff\nChanged\nLine 40: Radar returns is not a jargon. Measurements?\npreviously (see for example MacFerrin et al. 2019 Supplementary).\nLine 42: Technically the slabs themselves are superimposed ice and so yours are not the first observations\nWe are observing superimposed ice formation, i.e. the process of forming superimposed ice, as we state in the\ntext. Slab mapping has indeed already mapped superimposed ice.\nLine 43: What is the runoff limit?\nIn lines 23-24 we now define the visible runoff limit alongside the visible runoff area.\nLine 44: What is local runoff?\nChanged to ‘net water discharge from the hydrological catchment’.\nFig 1: Why is fs1 high on both years?\nThat there was no slush field at the surface doesn’t mean that there was no slush beneath the surface. Visible\nhydrological features do not tell the whole story, so we also need to look at the sub-surface. C-band SAR\nclearly shows that this area remained wet for a long time after the end of the melt season (Fig. 5a, see wetter\nfeature between FS2 and downstream later, between 1730 and 1760 m asl).\nLine 69: You are showing that the SIF varies by 10s of cm w e over short distances, so perhaps the smb does\ntoo, let alone surface features? In which case you experiment is not well constrained\nPlease see our responses to reviewer 1 concerning surface melt rates.\nLine 88: Why use e>3 when water is e=80?\nThe TDR measurements have been removed from this re-submission.\nLine 94: TDR probe\nThe TDR measurements have been removed from this re-submission.\nFig.3 : It’s very difficult in this figure to contrast the observed and modelled temperature profiles, which seem\nto be the only validation of the model. I suggest replotting them as e.g. a direct scatter plot with a colour\ncoding to illustrate depth so that we can see how well they compare (or otherwise)\nWe have redrawn this figure panel to make the comparisons between observed vs modelled temperature\nprofiles clearer. We have also added some text to the Results which addresses the offset between\nobserved and modelled profiles at the end of the melt season.\nLine 139: Isn't this obviously the case? What else would be causing the SIF. I don't think this is a grand\nconclusion\nWe are not exactly sure what `this’ refers to, as one sentence finishes and another starts at this line number.\nAssuming that ‘this’ refers to the sentence in L137-139, we do not make out that this is a grand conclusion.\nRather, this statement supported the sentence which followed, i.e. ‘In contrast, larger amount of SIF can only\nbe explained by lateral transport of meltwater….’. We have made various minor changes to better explain\nour model findings.\nLine 141: Just because your model agrees with SIF in some places doesn't mean you can conclude that\n14\nsomething else must be happening elsewhere – that’s not a hypothesis you have tested by experiment. Is your\nmodel able to simulate lateral movement of meltwater? What other processes might lead to local variations in\nmelting that are not simulated by your model? More work required here\nThe majority of cores show rather similar SIF. We clearly show that simulated background refreezing agrees\nto SIF in several locations. It clearly does not agree in a few locations. We therefore looked for the reasons for\nthe disagreement. It is obvious that refreezing from only the melting directly above a given location cannot\nvary so enormously in space given how homogeneous the snow-covered surface is and the rather small\nelevation envelope that we worked within (see also our response to reviewer 1). We also observe the largest\nrefreezing rates where there is most slush. To us, it is logical that these locations should therefore have less\ncold content in the slab. Finally, it is unclear to us what the benefit of including lateral movement in our 1-D\nmodel would be.\nLine 144: It does not indicate this; it just indicates you have a match within 1sigma in some locations\nThis line says ‘Since our model indicates that the localised larger quantities of superimposed ice cannot be\naccommodated by slab cold content during the melt season alone…’.\nOur model does indicate this as the modelled summertime SIF is valid for most of our measurement sites, not\njust the one where we compare. Of course, we do not have more firn temperature data with which to undertake\na more thorough comparison. But in the context of how much larger some of SI measurements are (by\ndecimetres), then we consider a match within 1 sigma (0.06 m w.e.) sufficient for our assertions.\nLine 156: Correlation, not relationship\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nFig. 4: These data are clustered in 2 populations which is leading to most of the correlation. I don't think they\nare distributed well enough to develop an empirical parameterisation. why are they clustered in this way? I\nsuggest this is explored in more detail before jumping to conclusions\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nLine 174: This is all a bit contrived; your result is entirely dependent upon the relationship that you force to\nexist.\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nFig. 6: This isn't \"wetness\", its sigma nought. I think you are jumping to too many conclusions here\nAs we describe in the Methods, we derive a binary wetness indication by applying a threshold to detrended\nsigma0HV. Thus, Figure 6 did present wetness (duration), not sigma0. This approach is based on\npreviously published work, most notably Liang et al. (2021). Note that Fig. 6 is now Suppl. Fig. 5.\nLine 211: So only half of the variance is explained\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 219: I'm lost here; if you have a model of the firn pack why are you counting cNDD to model SIF? why\nnot just compute it from your actual model?\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nWe now summarise the specific issues (difficulty locating the slab surface at FS3; uncertainty related to\ndrilling into SIF until metal chips appeared).\nLine 336: Why are they minimum estimates?\nPlease see response to previous comment.\nLine 343: This is a poorly constrained experiment as you don't know the spatial variation in firn stratigraphy.\nhow can you be sure you are measuring superimposed ice?\nWe believe that this concern is addressed by our more detailed explanation of our field measurement\nprotocol and Supplementary Discussion; specifically, during our repeat measurements, we found very little\nspatial variability at metres scales.\nLine 358: So this is an underestimtae of the SIF?\nThis paragraph is about the cSMB, not about SIF. By definition (Cogley et al. 2011), SIF refers only to\nrefreezing on top of impermeable horizons, in this case the ice slab.\nLine 459: Why choose the 1st percentile? is this based on theory?\nPlease see our response to Reviewer 2.\nLine 461: Why likely?\nThanks for pointing this out – on reflection, we do not think we can quantify this statement that the threshold\nis ‘likely conservative’, so have removed it in this re-submission.\nLine 466: Extrapolate or scale up, not distribute\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 470: This isn't determining the existence of a relationship, it is assuming the existence of a relationship\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 486: Detect not retrieve\nChanged.\n16"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Overall, I find that the authors have adequately addressed my comments.\nThe inclusion of supplementary Figure 1 is appreciated, as this partially clarifies some of the concerns of reviewer #3.\nI also appreciate the incorporation of the recent Culberg et al. 2024 results of ice slab maps into the resubmitted manuscript.\nI appreciate the removal of the statistical upscaling of remote sensing observations approach, as the large uncertainties,\nlimited correlation between SIF and crosspol backscatter, and lack of validation outside of the study area made these results\nsuspect. While I appreciate the new ‘bottom-up’ ‘top-down’ approach for upscaling SIF estimates, I do have some minor\nquestions about this new method.\n-Line 470: Could the authors briefly comment on to what degree the assumption of a constant slab thickness of 12 m has on\nthe bottom-up SIF calculation via the 1-D heat flow model now adopted in the revised manuscript? I appreciate that there is\nno existing estimate of slab thickness, but if we take as end members thicknesses of 5 m and 15 m (overestimate of 7 m vs.\nunderestimate of 3 m), how sensitive are the bottom-up SIF estimates to changes in the slab thickness assumed?\n-line 570: An analogous comparison can be drawn to the ‘top-down’ case, where the sensitivity in snowpack thickness is\nnow assessed. It appears an attempt at this was done with thicknesses of 0.5, 1.0, and 1.5, but this is not discussed any\nfurther. What were the results from this sensitivity analysis, and what snowpack thickness was ultimately assumed? The\n‘Localised SIF continues into autumn’ section suggests that a thickness of 0.6 m was ultimately used, and a comparison was\nmade to 0.3 m and 0.9 m, but this is seemingly contrary to what is described in the methods section, unless I misunderstand.\nClarification here would be helpful.", "reviewer_2": "Response to reviews\nConcurrent superimposed ice formation and meltwater\nrunoff on Greenland’s ice slabs, version 1\nSummary\nWe thank the three reviewers for their engagement with our first submission. We believe that we have\nsuccessfully addressed their major concerns, leading to major improvements which warrant re-submission\nthrough the Nature Communications appeals process.\nReviewer 1 found the manuscript to be well-written, thorough, and provide important observational evidence.\nTheir concerns were essentially minor. The bulk of our changes concern comments by reviewers 2 and 3.\nReviewer 2 said that the manuscript addresses an important caveat to the popular notion that ice slabs lead to\nrunoff. They found the field measurements novel and convincing but were concerned that the remote sensing\nanalysis had poorly quantified but likely large uncertainties. Reviewer 3 had some major concerns about both\nthe field measurements and the upscaling.\nOur most major and important changes are:\n1. We have developed a new upscaling methodology. We no longer rely on statistical links between\nfield-measured SIF, negative degree days and SAR. Instead, we extend our physical heat flux\nmodelling approach to estimate SIF through (a) ‘bottom-up’ (i.e. slab cold content) refreezing during\nboth summer and autumn, and (b) ‘top-down’ refreezing (i.e. negative air temperatures) during\nautumn. We still use SAR sigma0HV, but only to determine wetness duration both across the field\narea and ice-slabs-wide through application of a binary threshold. We undertake these calculations in\nboth our field area and ice-slabs-wide. It also allows a much more rigorous exploration of\nuncertainties. In our view this is a more elegant and physical approach, particularly because it allows\nus to fully differentiate between summer and autumn SIF. Please see our headline response to\nreviewer 2 for more details.\n2. Our ice-slabs-wide analysis now incorporates the new ice slab extents mapped by Culberg et al.\n(2024); as these extents are much larger than Jullien et al. (2023), our estimate of SIF onto ice slabs\nhas increased accordingly.\n3. We have recomputed our SAR wetness threshold and we have examined the sensitivity of our\ncalculations to the choice of the SAR sigma0HV threshold for mapping wetness; our findings are\nrobust across the specified sensitivity range of the Sentinel-1 SAR sensor.\n4. We disagree with the view of reviewer 3 that our in-situ stake measurements of SI are insufficient to\ndemonstrate our purported findings. In response we have extended our description of our SIF\nmeasurement protocol, including the addition of a sketch figure to the supplementary materials. Please\nalso see the rebuttal in this response document. We do however agree that our measurements to\ncapture SIF timing are not fully convincing, so have removed these accordingly.\n5. We have made many improvements to the main text which we hope improve the readability.\nOur changes in bold\nIn our longer responses, the key parts are underlined.\n1\nReviewer #1 (Remarks to the Author):\nThis manuscript presents estimates of superimposed ice formation on top of ice slabs along the K-Transect in\nSW Greenland. In this study, the authors present field measurements that suggest local ice formation on ice\nslabs can retain ¾ of local annual surface melt. They estimate that ice formation along ice slabs facilitated the\nretention of ~32 gigatonnes of melt between 2017 and 2022, roughly ~4% of the Greenland SMB.\nOverall, the manuscript is well-written, thorough, and presents important observational evidence for surface\nmelt refreeze onto ice slabs, which plays a role on the overall SMB of the GrIS. I have several comments,\nmost of which are relatively minor; I believe that their consideration will overall improve the manuscript.\nI have organized the comments below, grouped by section within the manuscript.\nResults:\nPg 3 line 77: Because local melt is near uniform across the study area (presumably as a result of their\nmodeling results, but this is not made explicit), they hypothesis that measured spatial variation in SIF is\nindicative of inputs not only from local melt production, but also lateral water flow. This seems entirely\nreasonable given the observed increases in SIF for sites proximal to major surface rivers. However, the\ninclusion of external references to reinforce this claim would be beneficial.\nAs the reviewer infers, we take melt as a smoothly varying field of decreasing melt with increasing elevation,\nas calculated by a surface energy balance model forced by in-situ observations. Fig 5d now illustrates how\nmelt decreases with increasing elevation. We also now refer to findings from Noel et al. 2016 in the first\nparagraph of our Results. They showed that melt rates show very little spatial variability in this region.\nIt is possible that, once visible slush fields have developed, then local melt rates might increase further as a\nfunction of their reduced albedo. However, these slush fields can only initially develop by the transport and\nponding of melt that is essentially uniform in space considering the smooth surface topography over 100s\nmetres scales and homogeneous snow surface cover (e.g. Holmes 1955; Noel et al., 2016; Clerx et al., 2024/in\npress). In other words, their origin does not lie in albedo. Nonetheless, we now note this possibility of\nsubsequent local melt enhancement in the first paragraph of the Discussion.\nPage 4 line 79: “Finally, beyond 2–4 km away from visible surface runoff, we found that SIF began to\nincrease again, as illustrated by 0.25 m w.e. at site 3u3 (Suppl. Fig. 1).” – With what degree of confidence can\nthis claim be made, given that it is reliant upon a single point wise measurement of SIF, and the measured\nvalue is within 2 standard deviations of the reported average SIF from 2022?\nOn reflection, we agree, so we have removed this statement from the main text.\nPage 7 line 139: “In contrast, larger amounts of SIF can only be explained by lateral transport of meltwater\nwhich provides additional water for SIF to continue in specific locations beyond the end of the melt season.” –\nSimilar to first comment; Isn’t it possible that this elevated SIF can also be explained by spatial heterogeneity\nin melt that is not captured in the SMB model? I tend to agree with the authors that the more plausible\nexplanation is lateral transport of meltwater, but to what degree can spatial heterogeneity of melt be ruled out?\nPlease see our response to the previous comment about spatial variability in melt.\nAn important point to mention in the refreezing during Autumn section is the role that volume scattering due\nto surface and subsurface heterogeneities (e.g. surface roughness, sastrugi, ice lenses, ice pipes, percolation\npipes) plays. Implicitly, the bias-removal approach that is employed here *should* account for this process,\nprovided the geometry of these features is assumed temporally non-varying, such that changes in the cross-\npolarization backscatter can be attributed solely to changes in dielectric permittivity rather than surface\nroughness and volume scattering. I would suggest that this point be elaborated upon in the methods section,\n2\nperhaps with reference to some relevant literature (e.g., Rignot 1995, Backscatter Model for the Unusual\nRadar Properties of the Greenland Ice Sheet).\nThanks for this suggestion. We have elaborated in the Methods (‘Detecting residual water presence’)\naccordingly.\nMethods:\nPage 16, line 396: Regarding the model producing 0.25 m w.e. melt; is this merely a consequence of the\nsubsurface model discretization being 0.25 m? If not, what is the uncertainty of this model-derived estimate of\nw.e. melt?\nNo, this is not a consequence of the sub-surface model discretisation. In each time step, the model’s grid\nmoves up/down after the occurrence of accumulation/ablation; van As et al. (2017) contains full details.\nConcerning uncertainty: the model does not come with uncertainty bounds because of the complex and non-\nlinear nature of the modelling. This makes validation a crucial aspect of SEB modelling, as is done in\nnumerous previous studies. Concerning this model specifically, see e.g. van As et al. (2017) cited in\nmanuscript. In this manuscript we therefore evaluate how sensible its melt estimates are by reference to in-situ\nsurface lowering measurements, which show excellent agreement (see Methods).\nPage 18 line 469: Section on the linear regression between measured SIF and mean crosspol backscatter:\ngiven that an uncertainty in the in situ SIF values are provided, it would be worthwhile to propagate this\nuncertainty through the linear regression for a nominal uncertainty estimate of the linearly-regressed area-\nwide SIF estimates, particularly since the R2 in the linear regression isn’t particularly strong.\nThe linear regression approach has been removed in this re-submission.\n3\nReviewer #2\nThis manuscript describes field measurements of superimposed ice formation on top of ice slabs on\nGreenland’s K-transect in the summers of 2021 and 2022. Using repeated firn cores, time domain\nreflectometry, and thermal modeling, the authors show that a significant portion of annual melt refreezes on\ntop of the ice slabs due to ice slab cold content and cold surface temperatures during autumn. In particular, the\nauthors show that refreezing occurs even where water flows laterally over the surface, potentially for kms,\nsuggesting that surface mass loss may not scale directly with visible surface runoff during the melt season.\nThe authors then use Sentinel-1 SAR measurements to establish relationships between their field\nmeasurements of SIF and radar-inferred surface wetness measurements. They use these relationships to scale\ntheir SIF measurements first to the area surrounding their field site and then to the entire ice sheet. The\nauthors conclude that SIF accounts for about 4% of the contemporary GrIS SMB, making it a small but non-\ntrivial contributor to the overall ice sheet mass balance.\nOverall, this paper addresses an important caveat to the popular notion that ice slabs lead to complete surface\nrunoff. In particular, the field measurements are novel and convincing and suggest that superimposed ice\nformation (SIF) needs to be accounted for in order to understand how quickly surface melting becomes mass\nloss as melting advances inland on the ice sheet. However, the remote sensing analyses used to scale these\nfield measurements to the entire ice sheet are subject to large uncertainties, many of which are not well\nquantified in the current manuscript. I am particularly concerned that the approach to estimating pan-\nGreenland SIF relies on a highly uncertain relationship between SIF and nCDD in a single ~30 km x 30 km\nregion which cannot be validated anywhere else on the ice sheet, as well a background SIF rate that is\nassumed to not vary at all across the ice sheet. This opens the potential for large, uncharacterizable regional\nbiases to affect the pan-Greenland SIF estimates.\n‘[the authors use a] background SIF rate that is assumed to not vary at all across the ice sheet’. In fact, the\nbackground SIF rate used for the ice-slabs-wide analysis varies in two ways:\n(1) As in our first submission, our treatment of the conditions which permit SIF across the ice sheet is\nexplicitly spatial. At a given MAR model grid node, SIF can only occur during the melt season and\nwhen the slab is covered in (porous) snow which provides the basis for a slush matrix overlying the\nslab.\n(2) Previously, we used three different fixed initial temperature scenarios to model bottom-up summer\nSIF. In this re-submission, we instead use the state-of-the-art ice-sheet-wide 10 m depth\ntemperatures data product published by Vandecrux et al. (2024, TC) to initialise our bottom-up\nSIF modelling. We also run our model at +/- the stated uncertainty of the Vandecrux product.\nWe consider this to be a major improvement compared to the previous version of our manuscript.\nConcerning the remote sensing analyses: we concur with the reviewer’s concerns that these were subject to\nlarge uncertainties and had the potential for large, uncharacterizable regional biases to affect the pan-GrIS SIF\nestimates. In conjunction with similar concerns raised by reviewer 3, we have therefore made major\nrevisions which replace the previous statistical approach completely with physical heat flow modelling\nto make our upscaling calculations. In summary:\n(1) In areas which remain wet after the end of the melt season according to SAR, we now allow our\nsummertime “bottom-up” calculations to continue into autumn until SAR indicates that the area is\ndry, yielding the additional bottom-up SIF which occurs in autumn. We partition the contributions of\nsummer vs autumn bottom-up SIF using the last day of melting from passive microwave radiometry.\n(2) In these same areas, we now also explicitly address “top-down” refreezing as a result of negative\nsurface temperatures, beginning once the melt season has finished. These calculations are very\nsensitive to the thickness of the snowpack overlying the slush matrix; we discuss this further in the\nmanuscript.\n(3) Any pixels suspected of being lakes (i.e. still ‘residually wet’ at the end of the year) are now\ncompletely excluded from our calculations.\n(4) We have recalculated the SAR sensitivity threshold using data acquired in all areas labelled as\n‘minimum likely’ slabs by Culberg et al. (2024, TC), with little resulting change (see below).\n4\n(5) Over our field area, we now show the uncertainty in our SIF calculations associated with the stated\nsensitivity of the Sentinel-1 SAR sensor used to determine wetness duration. We show that this source\nof uncertainty is second order, i.e. wetness duration does not change appreciably within this\nsensitivity range. We also checked this on an ice-slabs-wide basis where again it is a second-order\nsource of uncertainty (see Supplementary Figure 4).\nMajor Comments:\n[1] I am a bit confused about the primary purpose of modelling of the rate of superimposed ice formation. The\nauthors state that the reason they conduct the modelling is to determine the degree to which ice slab warming\nis due to latent heating. If this is indeed the main goal, then the argument would be better supported by\ncomparing the model run presented in the paper to a model run forced by the measured surface temperature\nand which only considers heat conduction (no latent heating). The difference between those two scenarios\nwould provide a stronger quantification of the role of latent heating in warming the ice slabs.\nHowever, the primary conclusion that seems to come out of this section is that refreezing due\nto ice slab cold content can only account for measured SIF rates in regions without additional\nadvection of meltwater. If the goal is instead to accurately model local SIF, then it seems like the authors\nmight be better served by using their fully coupled SEB-firn model instead of this simplified model that\nimposes a fixed meltwater supply that may or may not be consistent with the actual melt supply and timing\nduring the season when measurements were made.\nWe agree that we should have been clearer here. Our main goal was neither to super-accurately model local\nSIF nor to determine the degree to which ice slab warming is due to latent heating. Rather, it was to provide a\nconstraint on how much SIF can occur during the melt season under the plausible assumption of a continuous\nmeltwater supply at the slab surface. Then, SIF in excess of this logically occurs during autumn, i.e. after\nactive surface melt has finished. In other words, the primary purpose of this modelling approach is to partition\nthe SIF that we measured into the melt and post-melt seasons, thereby setting the scene for a slabs-wide first-\norder estimate of SIF later in the study. With our revised modelling approach we believe that this is\nmotivation is much clearer.\nAn ideal solution would indeed be to use a fully coupled SEB-firn model. However, no such model that is\nsuitable for our application exists. We are dealing with the problem of ponded meltwater in a snow matrix\natop an impermeable and cold barrier, which is an extreme end-member of what can be considered ‘firn’. We\nunderstand from firn modelling colleagues that even when ponding subroutines are available then they often\nremain disabled for want of validation. We now describe these issues in the Methods section ‘Modelling\nSIF’. Thus, it is more informative to simplify the problem as far as possible, i.e. by considering a relatively\nthin snowpack with a continuously wet slush matrix on top of an ice slab (conditions which were met for\nmuch of the melt season according to our in-situ measurements and further bolstered by observations in Clerx\net al., 2022) and then using classical heat flux theory to determine the rate of refreezing at first order.\nNew in this re-submission is an alternative configuration of the model to estimate SIF caused by\natmospheric cooling (‘top-down’). As we describe in the paper, an issue here is that we lack good bounds on\nthe thickness of the snowpack overlying the slush matrix. We wanted to use end-of-season snowpack\nthickness simulated by MAR but found that the outputs were highly spatially incoherent, which is a further\ncaution against 'blind’ use of RCMs assuming that they always deliver the best results. The reviewer may be\ninterested to see a paper with several of the same co-authors as this submission currently in Discussion, which\ndescribes some of the difficulties with RCM behaviour at the runoff limits (https://doi.org/10.5194/egusphere-\n2024-2750).\nOverall, we contend that our chosen approach is much simpler to interpret than the outputs of a fully coupled\nSEB-firn model. In turn, we hope that our findings can be used to advance ponding schemes in firn models.\nWe have made narrative changes to the relevant parts of the manuscript accordingly.\n5\n[2] Uncertainty in the regression between measured SIF and the multi-temporal mean 𝜎̂ 𝐻𝑉 can be better\nquantified and mapped into uncertainty in the study site SIF map (Fig. 4b). One option would be to use a\nBayesian regression that can quantify the range of linear models that are consistent with the data. This would\nmake it possible to produce an uncertainty map showing the range of the plausible SIF for each pixel in Fig.\n4b given the range of plausible models. But the key point is that there should be corresponding uncertainty\nmap for Fig 4b that account for uncertainty in the SIF-𝜎̂ 𝐻𝑉 relationship and the uncertainty in the background\nSIF rate.\nSimilarly, better uncertainty quantification is needed for the SIF-nCDD regression. Particularly given the\nrelatively poor R2 value for this regression, I strongly feel that it is important to quantify the uncertainty in the\nmodel fit and clearly demonstrate how that uncertainty maps onto uncertainty in the final estimates of regional\nautumn SIF that are presented in Figure 7.\nThis analysis has been removed from our re-submission.\n[3] For setting the SAR wetness backscatter threshold, I would like to see some sensitivity tests from different\nregions of the ice sheet showing how the 1st percentile of 𝜎̂ 𝐻𝑉 prior to the melt season varies regionally. For\na reliable pan-Greenland mapping, it is very important to establish whether conditions in SW Greenland are\nconsistent across the ice sheet. Other studies that have used Sentinel-1 data to quantify either melting or the\nextent of wet snow facies over Greenland have found that performance with a single threshold can vary\nsignificantly across the ice sheet (Hu et al., 2022; Li et al., 2023). Regional bias due to a fixed wetness\nthreshold should be quantified and incorporated into assessments of overall uncertainty.\nWe thank the reviewer for raising these concerns. First, we comment on the literature citations:\n(1) The findings of Hu et al. 2022 are not as relevant to our study as they first appear. Hu et al. did not\nnormalise summertime retrievals with respect to winter, unlike in our current approach, and so\ntherefore a considerable source of variability remains.\n(2) Li et al. 2024’s methodology is more similar to ours. They emphasise that SAR-based mapping of\nwetness at a set threshold is only reliable above the snowline, but that there is no appreciable regional\nvariability above the snowline. As we consider only snow-covered surfaces above the ELA, the\ndifficulties that Li et al. experienced with selecting a threshold applicable to melt detection over bare\nice is not relevant here.\nWe recap that our threshold of -2.1 dB was quite similar to Liang et al. 2021, who rather than the 1st percentile\ntake the less conservative 5th percentile of dry snow as their threshold: -2.66 dB. Both Liang et al., 2021 and\nLi et al. 2024 employ a more conservative value of -3 dB, +/- 0.3 dB, but they had the express aim of\ncapturing surface melting conditions which they validated by reference to in-situ air temperature\nmeasurements. This is not a relevant validation metric in our case as we are seeking to identify residual\nwetness after the end of active melt, which also indicates that a 3 dB threshold would be too conservative for\nour purposes here.\nNotwithstanding this comparison, we have undertaken another examination of our threshold choice. Now, we\ncompute histograms by sampling all areas within the ‘minimum likely ice slab’ classification of Culberg\net al. (2024), rather than just the western flank of the ice sheet as previously. The 1st percentile value\nhas decreased slightly to -1.17 dB (i.e. 0.07 dB more conservative than in our previous submission).\nIn the revised manuscript we adopt the new threshold of -1.17 dB, +- 0.3 dB indicative of the\nradiometric uncertainty of Sentinel-1’s sensor as per Geudtner and Torres (2012). We use the +/- to\npropagate an uncertainty in wetness duration through to our final calculations in our field area. On an\nice-slabs-wide basis, we provide a new Supplementary Figure which shows that the residual wetness\nduration/area is essentially insensitive to the radiometric uncertainty of the sensor.\n[4] I am concerned about how well the nCDD model can be generalized to map autumn SIF rates across the\nentire ice sheet. This is conceptually quite similar to the positive degree day model used for melting, which\n6\nhas a spatially and temporally variable rate factor that typically needs to be tuned for individual sites. I have\nsimilar concerns for the nCDD model. Essentially, the authors’ scaling assumes that air temperature is the\nonly driver of refreezing and neglects all the other components of the subsurface energy balance, to include\nfirn and snow conditions, which might vary significantly across the ice sheet. It also assumes that the amount\nof melt that is refrozen scales directly with the length of time required to fully refreeze an area, which may be\ngenerally reasonable, but would absolutely be influenced by the transient accumulation rate, subsurface\ntemperatures, etc on the local scale. Without secondary validation sites, it seems basically impossible to know\nif the regression model coefficients will apply anywhere outside of southwest Greenland.\nThis analysis has been removed from our re-submission.\n[5] How are buried lakes accounted for when calculating autumn SIF? If I have understood the methods\ncorrectly, it seems like there has been no attempt to separate that mode of water\nimpoundment from superimposed ice formation. It seems problematic to apply a relationship\nbetween nCDD and refreezing in porous snow or a slush matrix to the refreezing of large open water bodies.\nSince buried lakes will have the highest and most persistent wetness values (see the authors’ own Figure 5,\nwhich is dominated by lakes), I am concerned that including lakes may significantly bias the results.\nThanks for raising this concern. Now, we screen out likely lakes by removing pixels which are still wet in\nthe last SAR acquisition each calendar year.\nMinor Comments:\nOverall, the section titles and introductory sentences could do a better job of explaining to the reader the end\ngoal of the analyses being described in the section. For example, the section “Modelling the rate of\nsuperimposed ice formation” would be easier to follow if it started with a sentence explaining why the authors\nwant to model this quantity. The section title could also better reflect the findings instead of the methods,\nfor example, something like “Ice slab cold content a key factor in summer SIF”.\nThanks for this suggestion. We have made several changes accordingly.\nLine 34 – I am not sure I understand what the second half of this sentence is trying to express.\nWe have rephrased as we agree that this was very unclear previously.\nLine 92 – Why was there no spike in 𝜖 when this water reached and fully flooded the probe before refreezing?\nWe cannot know for sure, but as the section ‘Timing of SIF’ has been removed in response to reviewer 3\nconcerns, we have not investigated further.\nLine 185/296 – there are now estimates of ice slab extent from both modeling (Brils et al., 2024) and Sentinel-\n1 observations (Culberg et al., 2024) that provide a more up-to-date and comprehensive outline of ice slab\nextent on the GrIS.\nIndeed - though we find it a little concerning that literature published in 2023 can already be considered less\n‘up-to-date’! Nonetheless, given that we have made major changes to our upscaling approach, we have taken\nthe opportunity to replace the Jullien et al. (2023) outlines with those based on Sentinel-1 observations\nby Culberg et al. (2024). We now explore the impact that the uncertainty in ice slab extent makes on our final\nestimates of SIF.\nWe have moved the statements of ice slab extent previously found at this location into the Introduction.\n7\nLine 223 – modeled temperature profile, for example, from IMAU FDM with ice slabs (Brils et al., 2024)\nwould at least provide a constraint on the shape and range of temperature profiles in ice slab regions.\nThanks for this suggestion. Rather than using modelled temperature profiles for initialisation with all the\nuncertainties in RCM sub-surface schemes, we take advantage of another study that has been published this\nyear: Vandecrux et al. (2024, TC), which mapped observed 10 m depth temperatures ice-sheet-wide using\nmachine learning techniques, presented along with an associated uncertainty product. We use the Vandecrux\ndata to set the bottom temperature of our initial temperature profiles and then scale our observed\nprofile from FS2 accordingly. We also run our calculations +/- the stated 10 m temperature uncertainty.\nPlease see our revised Methods section for more information.\nHow is the long-term equilibrium line in Figure 5 defined?\nThis was based on RACMO data 1981-2000. But in this re-submission we have modified this figure to only\nshow wetness in the Culberg et al. (2024) maximum ice slabs extent area.\nLine 297 – I found the context for the discussion future warming and peak refreezing a bit confusing. Is the\npoint simply that we should not project current SIF rates as a function of SMB into the future due to the\nuncertainty in climate forcing and hypsometry?\nWe have made efforts to improve the narrative of this paragraph – please see full text.\nLine 476 – where did n=4 come from when only 3 PROMICE weather stations were noted in the study area?\nWhich stations were used and how are they distributed in elevation across the study area? (e.g. is it fully\nrepresentative of the elevation distribution, or does this scaling do significant extrapolation outside the data\nbounds?)\nThe middle and upper weather stations largely encompass our study area, with the upper station being situated\nat 1840 m ASL. Above this elevation SEB model input variables are extrapolated, not interpolated, and\ntherefore calculations do not exceed 2000 m ASL so that extrapolation does not become unrealistic. The\nmodel uses these observations to compute in discrete 100 m elevation bins. We computed our regression\nacross four elevation bins (1,600, 1700, 1800 and 1,900 m asl). This was not clear previously, so thanks for\ndrawing this to our attention and we have updated the text accordingly.\nLine 506 – What determines if a cell becomes snow free? Is this from MAR, an albedo remote sensing\nproduct, or something else?\nAccording to modelled MAR snowpack. This is now noted in the text.\nLine 526 – wouldn’t a brief rise above 0C potentially lead to melting that could supply more water for SIF?\nWe have added a note that, given several weeks of prior atmospheric cooling of the snowpack, it is probable\nthat this melt would simply have been refrozen in the uppermost centimetres below the surface.\n8\nReviewer #3 (Remarks to the Author):\nSummary\nThis paper uses a combination of field observations, satellite observations, and simple/conceptual modelling\nto address the important topic of meltwater retention in the percolation zone of the Greenland ice sheet –\nabout which there is some uncertainty. The paper builds on the finding that meltwater retention has led to the\nformation of ice slabs at depth the firn pack, and that the presence of these slabs disrupts the normal process of\nmeltwater percolation and redistribution.\nThe paper is quite challenging to comprehend as it contains several elements of work that are not easily linked\nand are each quite obscure. Despite several attempts to read and review, I could not form a clear narrative\nfrom the study and, more importantly, I was not persuaded by the arguments presented. Either the narrative\nneeds to be strengthened, or the arguments need to be made more compelling with further work.\nMajor\nI have identified five major concerns that I hope will act as a guide for the authors to reformulate their paper.\n1. Loose language. There is lots of loose language and (for me) too much hyperbole in this paper. The term\n“runoff” in particular is central to this work but is not well or consistently defined; it is interchangeably used\nto describe meltwater, mobile water, and water that has left the ice sheet. In my book it is the latter only.\nElsewhere the term “rivers” is used but this is informal and not helpful; it would be better to make an\nassessment of the intensity of meltwater flow in supraglacial channels and contrast the magnitude.\nAnd it is stated in the abstract and repeated elsewhere that the amount of runoff from the percolation zone is\nnot known when this is simply not true – the authors used for example regional climate models that make an\nestimate of this and there have been many studies on meltwater retention including the identification of ice\nslabs on which this work is predicated. I would recommend a thorough parse of loose and hyperbolic\nlanguage.\nWe have endeavoured to tighten up our use of language throughout the manuscript, particularly around\nthe term ‘runoff’ which we agree was not used as precisely as it should have been. We hope that these changes\nmeet with approval.\nThe reviewer said that they find ‘too much hyperbole’. The only example they provided is: “it is stated in the\nabstract and repeated elsewhere that the amount of runoff from the percolation zone is not known when this is\nsimply not true – the authors used for example regional climate models that make an estimate of this and there\nhave been many studies on meltwater retention including the identification of ice slabs on which this work is\npredicated.”. Concerning this specific example, we disagree: regional climate models make estimates which\nso far remain unvalidated for the percolation zone, especially where ice slabs are present. They are potentially\nsubject to major issues associated with invalid parameterisations. We currently have a manuscript co-written\nwith the MAR and RACMO teams under discussion which covers many of these themes in detail:\nhttps://doi.org/10.5194/egusphere-2024-2750. Nonetheless, to tighten the language we are happy to change\nfrom ‘unknown’ to ‘unmeasured’, ‘unvalidated’. We also note that we have added further justification\nfor our heat flux modelling approach (over RCMs) in the Methods section ‘Modelling SIF’.\nThe reviewer states that there have been ‘many studies on meltwater retention including the identification of\nice slabs:\n• The reviewer does not provide examples of these ‘many studies’.\n• Taking this comment literally, we think the reviewer is suggesting that there are already several\nstudies of meltwater retention on top of ice slabs. We are not aware that such studies exist, with the\nexceptions of already-cited studies of SIF atop glacier ice in the Canadian Arctic and Svalbard. In\n9\nGreenland, we are only aware of slab mapping/identification exercises, e.g. MacFerrin et al. (2019),\nJullien et al. (2023), Culberg et al. (2024), which do not simultaneously assess meltwater retention\nprocesses.\n• There is arguably perhaps one study, Machguth et al., (2016), which focuses on meltwater retention\non a slab. In this study, the term ‘superimposed ice formation’ appears only once, in the 3rd-last\nparagraph of the Methods, only in reference to core alignment. This study undertook ‘barcode\nalignment’ on a couple of firn cores at KAN_U, which as this reviewer points out subsequently can be\nsubject to substantial uncertainties. But having re-read this study in detail for this response we have\nconcluded that the evidence of SIF they present is so limited that we now do not cite it in the\nIntroduction as an observation of Greenland SIF.\n2. Superimposed ice formation. The initial part of this work is based on in situ measurements that purport to\ndocument superimposed ice formation. I am not persuaded that the authors have demonstrated this effectively.\nThe main approach is to measure the depth to the slab surface two years apart at a variety of locations. This is\nnot an actual measurement of superimposed ice formation, it’s a crude approximation. The slabs could be\nmobile, because the site is not at the equilibrium line. More importantly, the authors are not able to resurvey\nthe same site and so their measurement is in fact a combination of spatial and temporal variation in slab\nsurface depth. Without knowing the former, it is impossible to deduce the latter. This calculation\nunfortunately underpins the entire study, and so it is quite fundamental to convince the reader that the\nmeasurements are robust. Having drilled many firn cores in this region I know that measuring the depth to a\ncommon horizon is a very challenging thing to do with precision, and that there are typically large natural\nvariations in the depth of common horizons in any case. I wonder if repeat GPR would have been a better way\nto do this entire experiment. Or drilling through the entire slab. Or analysing the core chemistry to isolate\nannual layers. Certainly the spatial variation in slab surface depth needs to be reliable quantified somehow\nbefore the temporal variation can be deduced.\nWe strongly disagree with the thrust of this comment that our in-situ observations are invalid. We highlight\nthat neither of the other two reviewers raised concerns of this nature. From our perspective there seems to\nhave been a misunderstanding of our field measurement protocol. At root, this seems to relate to the locations\nat which we make our measurements. We do resurvey the same site between years: these locations are exactly\nthe same each year with respect to the underlying stratigraphy, as we returned to a location marked with a\nbamboo stake. Each location had of course moved by several 10s of metres due to ice flow, but this process\ndoes not affect the near surface stratigraphy. To clarify this further, we have improved our explanation,\nincluding adding a supplementary figure which graphically illustrates the measurement process, and\nsupplementary discussion to discuss local variability in SIF more extensively.\nIn our judgement the suggestion of using GPR is a non-starter. If we were to make a GPR survey, we would\nmeasure the thickness of the snow layer, albeit with larger uncertainty than by drilling. We ask ourselves: how\ndoes GPR measure the thickness of the slab and in particular its annual growth? We could potentially detect\nthe bottom depth of the slab in subsequent years, but would be hampered by all the GPR-related uncertainties\nthat would reduce the precision of the outcome in such a study. And crucially, even if a change in thickness\ncould be measured beyond uncertainty, we still would not know whether the slab grew from accretion below\nor above. In contrast, this is precisely what our approach of returning to a marked surface location in year 2\nenables us to do.\nThe reviewer states “Having drilled so many firn cores in this region I know that measuring the depth to a\ncommon horizon is very challenging…”. We have no way of assessing what expertise this reviewer has, nor\nthe scope of the region they are considering. The reviewer may have drilled firn cores at high elevations, but\nperhaps not into ice slabs. Our cores are from the slabs, so are basically solid ice. In these cores, to say where\nsnow ends and firn starts is not always clear, but the horizon where the slab begins is always abundantly clear.\nFor more details, the core logs can be examined in our data repository.\n3. Spatial variation in SIF. A main conclusion of this work is that the degree of superimposed ice formation is\nhigher closer to surface meltwater channels than it is distal from them. This is I think the most original aspect\nof the study. However, the conclusion itself depends on how reliable the estimates of SIF (see point 2),\nbut\n10\nalso on how strong the spatial correlation is. Looking at the data in Figure 1 I am not convinced that the\nspatial variation in SIF is significantly correlated with the distribution of meltwater channels; there seem to be\njust as many sites that buck the trend as there are that fit the hypothesis. The sites sampled are also quite\nrandomly distributed, which doesn’t help. Wouldn’t it make sense to survey a profile tangential to a channel,\nfor example? I think more work and quite possibly more data are needed to convince the reader that this\nconclusion is robust.\nAs the reviewer notes, “the conclusion itself depends on how reliable the estimates of SIF (see point 2)” are.\nAs we respond above, we are confident in our SIF estimates and believe that there has been a\nmisunderstanding here.\nConcerning spatial variation, Fig. 1 provides a high-level overview of the SIF magnitudes. Further evidence is\npresented in Supplementary Figure 1, which explicitly compares the quantity of SIF with the Euclidean\ndistance to surface drainage visible on optical satellite imagery. Our measurements were concentrated on two\ntransects, with a spacing that reflects the proximity to surface drainage. Sites FS1, 122, camp, 231, 3u3\nconstitute one transect. The sites shown in Fig. 1d constitute another transect, which is quasi-tangential to the\nmain surface channel. We now note these transect configurations in the main text. Concerning a profile\ntruly tangential to a channel: the location of the channel varies somewhat between years, so this is less feasible\nto deploy than it might first appear, and also what a transect equidistant of a stream is one year, would not be\nthe next. Concerning the measurement stations, these were installed to acquire data in three distinctively\ndifferent surface slope regimes and with contrasting surface hydrology in any given melt summer.\nPlease see also our response to reviewer 1 concerning possible spatial variability in melt.\n4. Timing of SIF. There is a section of work on the timing of superimposed ice formation deduced from\nmeasurements of the firn dielectric properties. I found this section very difficult to navigate and not\npersuasive; the only thing that seems to be consistent about the probe measurements is how inconsistent they\nare in time and space. It seems that very few measurements were actually made and yet apparently firm\nconclusions are drawn from them. This is very similar to point 3; it suggests to me that the field experiment\nwas again not up to the task and in hindsight other measurements would have been made. It’s also not clear to\nme how important it is to know the timing of superimposed ice formation – is that really material to this study\n? My feeling is this whole is a distraction and could easily be dropped.\nOn reflection, we agree that knowing the sub-seasonal timing of SIF is not important to our key messages and\nis challenging to interpret from our data. We have therefore removed this whole section as suggested.\n5. Modelling. The study next leaps into a series of “model” experiments designed to help scale up the field\nobservations. However, these experiments are tangentially linked to the field measurements at best. The first\n“model” is a physical 1D heat transfer model, the second is a linear regression between satellite radar\nbackscatter and measured superimposed ice formation, and the third is another linear regression between the\nmodelled superimposed ice formation and cumulative negative degree days derived (I think) from the 1D heat\ntransfer model. The heat flow model is very loosely related to the field observations as it is contrasted to\ntemperature profiles that are not directly related to the SIF. Moreover, we can’t really judge how well the\nmodel and observations agree from the data as presented.\nWe are grateful for these comments, which together with similar comments by reviewer 2 caused us to\nfundamentally rethink how we conduct our upscaling. We believe that we have now successfully\naddressed this concern and that our response to reviewer 2 answers this major comment.\nThe satellite radar backscatter is very loosely related to in situ superimposed ice formation; many other factors\naffect backscatter and in any case the correlation deduced is largely owing to the clustering of the\nmeasurements into two distinct populations, without which there would be no correlation.\nWe recognise that in our previous submission there was the potential to conflate radar backscatter with in-situ\nSIF, primarily because of the statistical way in which we attempted to conduct our upscaling; however, this\n11\nwas not our intention. To be clear, we consider backscatter variations in autumn to tell us only about residual\nwetness in the near-surface. It is the disappearance of this backscatter that we attribute to SIF, which we now\nmodel physically instead of our previous problematic statistical approach. In the revised manuscript this is\nnow clearer because we talk only of residual wetness [duration] and not also in terms of sigmaHV0\nbackscatter strength like previously.\nThe degree of correlation between satellite derived SIF and negative degree days is relatively weak in any\ncase, half of the variation is not explained. I find this approach to be quite chaotic, highly empirical, and not\npersuasive. I would recommend the authors rethink their approach here entirely and find a better way to add\ncontext and value to their sparse field measurements. The authors do note that regional climate models with\ncoupled firn process models are another way to study meltwater retention and I agree this would be a better\nway to approach this problem.\nAs explained in our response summary, we have taken heed of these concerns about the upscaling\napproach and, we believe, now have a much better methodology with which to deliver our upscaled\ncalculations.\nConcerning the suggestion here that RCMs with coupled firn process models would be a ‘better way’ to\napproach the problem: we disagree, and also note that we do not state that they constitute a ‘better way’.\nAs we have tried to clarify in our re-submission, there are substantial uncertainties/problems remaining with how\nRCMs and even observation-based SEB models represent the sub-surface, which makes it desirable to instead\nadopt a simplified approach to modelling SIF. Please see in particular the Methods section ‘Modelling\nsuperimposed ice formation’.\nMinor comments\nLine 1: I’m not sure the word rivers is helpful here as they are generally associated with land and rain fed. can\nyou describe in glaciological terminology instead?\nThe term ‘rivers’ is used elsewhere in the literature. Non exhaustive examples include line 2 of Yang and\nSmith (2016, JGR); Smith et al. (2015, PNAS). We have decided to stick with this intuitive terminology rather\nthan ‘supraglacial streams and rivers’ at every mention.\nLine 4: I don't think its unknown; there are estimates. need to be a bit clearer here what you mean\nWe have changed to ‘unmeasured’.\nLine5: should probably be clear about what you mean by runoff; is it meltwater until it is mobile, or until it is\nknown that it will leave the ice sheet locus, or until it has left the ice sheet locus?\nHere, in-line with the definition of runoff, water that is assumed to leave the ice sheet locus – updated to\nrefer to ‘visible runoff’ which is then linked to previous studies in the Introduction.\nLine 6: I think you mean that up to 3/4 of the annual melt can be retained by SIF (i.e. the converse of what is\nstated), or that the process of SIF can lead to meltwater being retained\nWe are a confused by this comment – the text said that ‘SIF….can retain up to roughly three quarters of\nlocally-produced annual melt’, which doesn’t seem the be converse. Nonetheless, we have tried to clarify in\nour revised Abstract.\nLine 7: Is the heat flow modelling informed somehow by the field observations? or how else are they linked?\nThis comment refers to the abstract. We have added a note that the modelling is initialised by in-situ field\nobservations. The full information is available in the main text as before.\n12\nLine 7: How are the SAR data used? would be helpful to state\nWe have added a note that they are used to derive residual wetness.\nLine 8: Quantify \"majority\"\nWe prefer not to do this in the abstract, as while the numbers are always a majority (i.e. greater than half),\nthe\nprecise value depends on parameter choices. The reader should read the main text for these details.\nLine 8: Quantify \"remainder\"\nSee previous comment.\nLine 9: This probably answers my earlier point about the definition of runoff\nYes, we believe so; see also modified abstract.\nLine 10: Using your model?\nWe are not sure what this comment refers to, sorry.\nLine 11: Does this mean total or local SMB? not sure SMB is a good metric since it is the sum of\naccumulation and ablation\nThis is a good point – we concur that SMB was a poor choice of metric. We now instead express this in\nterms of runoff estimates made using Cryosat-2 altimetry measurements. We estimate that recent SIF\nhas decreased recent runoff by 2-3 %.\nLine 12: Is there an uncertainty on 32 Gt?\nThis value has changed with our new upscaling approach; we now provide a ‘most-likely’ value (60 Gt)\ntogether with a range (27-75 Gt). The range is primarily a function of which ice slab extent we examine, but\nalso varies depending on 10 m firn temperatures and autumn snowpack thickness.\nLine 12: To call them hotspots this is relative to other areas, but it doesn’t seem that you have assessed those?\nConsidering refreezing in absolute terms (metres water equivalent), they are indeed hotspots compared to the\nrest of the ice sheet. Nowhere else is this combination of (a) meltwater supply and (b) a snow/firn\nmatrix in\nwhich to refreeze onto the underlying ice met.\nLine 19: Decreased not densified\nChanged\nLine 25: Is this really true? i think there are estimates in the models you are using- what I guess you are saying\nis that they don't account for potential retention associated with the ice slabs. Thats very different to unknown\nWe have replaced with ‘unmeasured’.\nLine 37: Running water not runoff\nChanged\nLine 40: Radar returns is not a jargon. Measurements?\n13\nWe have changed to “measurements” as suggested, but note that this term has been used in the literature\npreviously (see for example MacFerrin et al. 2019 Supplementary).\nLine 42: Technically the slabs themselves are superimposed ice and so yours are not the first observations\nWe are observing superimposed ice formation, i.e. the process of forming superimposed ice, as we state in the\ntext. Slab mapping has indeed already mapped superimposed ice.\nLine 43: What is the runoff limit?\nIn lines 23-24 we now define the visible runoff limit alongside the visible runoff area.\nLine 44: What is local runoff?\nChanged to ‘net water discharge from the hydrological catchment’.\nFig 1: Why is fs1 high on both years?\nThat there was no slush field at the surface doesn’t mean that there was no slush beneath the surface. Visible\nhydrological features do not tell the whole story, so we also need to look at the sub-surface. C-band SAR\nclearly shows that this area remained wet for a long time after the end of the melt season (Fig. 5a, see wetter\nfeature between FS2 and downstream later, between 1730 and 1760 m asl).\nLine 69: You are showing that the SIF varies by 10s of cm w e over short distances, so perhaps the smb does\ntoo, let alone surface features? In which case you experiment is not well constrained\nPlease see our responses to reviewer 1 concerning surface melt rates.\nLine 88: Why use e>3 when water is e=80?\nThe TDR measurements have been removed from this re-submission.\nLine 94: TDR probe\nThe TDR measurements have been removed from this re-submission.\nFig.3 : It’s very difficult in this figure to contrast the observed and modelled temperature profiles, which seem\nto be the only validation of the model. I suggest replotting them as e.g. a direct scatter plot with a colour\ncoding to illustrate depth so that we can see how well they compare (or otherwise)\nWe have redrawn this figure panel to make the comparisons between observed vs modelled temperature\nprofiles clearer. We have also added some text to the Results which addresses the offset between\nobserved and modelled profiles at the end of the melt season.\nLine 139: Isn't this obviously the case? What else would be causing the SIF. I don't think this is a grand\nconclusion\nWe are not exactly sure what `this’ refers to, as one sentence finishes and another starts at this line number.\nAssuming that ‘this’ refers to the sentence in L137-139, we do not make out that this is a grand conclusion.\nRather, this statement supported the sentence which followed, i.e. ‘In contrast, larger amount of SIF can only\nbe explained by lateral transport of meltwater….’. We have made various minor changes to better explain\nour model findings.\nLine 141: Just because your model agrees with SIF in some places doesn't mean you can conclude that\n14\nsomething else must be happening elsewhere – that’s not a hypothesis you have tested by experiment. Is your\nmodel able to simulate lateral movement of meltwater? What other processes might lead to local variations in\nmelting that are not simulated by your model? More work required here\nThe majority of cores show rather similar SIF. We clearly show that simulated background refreezing agrees\nto SIF in several locations. It clearly does not agree in a few locations. We therefore looked for the reasons for\nthe disagreement. It is obvious that refreezing from only the melting directly above a given location cannot\nvary so enormously in space given how homogeneous the snow-covered surface is and the rather small\nelevation envelope that we worked within (see also our response to reviewer 1). We also observe the largest\nrefreezing rates where there is most slush. To us, it is logical that these locations should therefore have less\ncold content in the slab. Finally, it is unclear to us what the benefit of including lateral movement in our 1-D\nmodel would be.\nLine 144: It does not indicate this; it just indicates you have a match within 1sigma in some locations\nThis line says ‘Since our model indicates that the localised larger quantities of superimposed ice cannot be\naccommodated by slab cold content during the melt season alone…’.\nOur model does indicate this as the modelled summertime SIF is valid for most of our measurement sites, not\njust the one where we compare. Of course, we do not have more firn temperature data with which to undertake\na more thorough comparison. But in the context of how much larger some of SI measurements are (by\ndecimetres), then we consider a match within 1 sigma (0.06 m w.e.) sufficient for our assertions.\nLine 156: Correlation, not relationship\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nFig. 4: These data are clustered in 2 populations which is leading to most of the correlation. I don't think they\nare distributed well enough to develop an empirical parameterisation. why are they clustered in this way? I\nsuggest this is explored in more detail before jumping to conclusions\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nLine 174: This is all a bit contrived; your result is entirely dependent upon the relationship that you force to\nexist.\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nFig. 6: This isn't \"wetness\", its sigma nought. I think you are jumping to too many conclusions here\nAs we describe in the Methods, we derive a binary wetness indication by applying a threshold to detrended\nsigma0HV. Thus, Figure 6 did present wetness (duration), not sigma0. This approach is based on\npreviously published work, most notably Liang et al. (2021). Note that Fig. 6 is now Suppl. Fig. 5.\nLine 211: So only half of the variance is explained\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 219: I'm lost here; if you have a model of the firn pack why are you counting cNDD to model SIF? why\nnot just compute it from your actual model?\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\n15\nLine 334: Such as?\nWe now summarise the specific issues (difficulty locating the slab surface at FS3; uncertainty related to\ndrilling into SIF until metal chips appeared).\nLine 336: Why are they minimum estimates?\nPlease see response to previous comment.\nLine 343: This is a poorly constrained experiment as you don't know the spatial variation in firn stratigraphy.\nhow can you be sure you are measuring superimposed ice?\nWe believe that this concern is addressed by our more detailed explanation of our field measurement\nprotocol and Supplementary Discussion; specifically, during our repeat measurements, we found very little\nspatial variability at metres scales.\nLine 358: So this is an underestimtae of the SIF?\nThis paragraph is about the cSMB, not about SIF. By definition (Cogley et al. 2011), SIF refers only to\nrefreezing on top of impermeable horizons, in this case the ice slab.\nLine 459: Why choose the 1st percentile? is this based on theory?\nPlease see our response to Reviewer 2.\nLine 461: Why likely?\nThanks for pointing this out – on reflection, we do not think we can quantify this statement that the threshold\nis ‘likely conservative’, so have removed it in this re-submission.\nLine 466: Extrapolate or scale up, not distribute\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 470: This isn't determining the existence of a relationship, it is assuming the existence of a relationship\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 486: Detect not retrieve\nChanged.\n16\nResponse to reviews\nConcurrent superimposed ice formation and meltwater runoff on\nGreenland’s ice slabs, version 2\nWe are grateful to the reviewers for their time engaging with our revised submission. In response we\nhave made further substantial changes which we believe have improved the manuscript.\nMain changes:\n• We have established that bottom-up SIF is not sensitive to slab thickness within observation-\nbased bounds (R#1, R#2).\n• We have added a sensitivity analysis to the supplementary materials about the potential\nimportance of summer top-down SIF – it is not large in comparison to other SIF sources and so\nwe can neglect it (R#2).\n• While considering summer top-down SIF, we realised that we should also incorporate refreezing\nof irreducible water content in the overlying snowpack in autumn. We have re-run our\ncalculations for both the field study area and slabs-wide with this process included. This\nreduces autumn top-down SIF and in turn slightly reduces total SIF compared to our previous\nsubmission. See our response to R#2.1 for more information.\n• We have improved several Figures concerning the field study (R#3, R#2). In particular, Fig. 1 is\nclarified based on R#3’s suggestions and now incorporates distance to drainage. Fig 4 has\nimproved labelling and now focuses only on autumn SIF. We provide a new Fig. 5 showing total\nSIF across the study area and including observed vs. modelled SIF.\n• When examining SIF at each stake, we now consider slope angle in addition to distance to\nsurface drainage (R#2, R#3). See Fig. 1d and the rewritten Results section ‘Superimposed ice\nformation depends on near-surface water flow’.\n• For our field study area we now dedicate “Localised SIF continues into autumn’ to presenting\nFig. 4, while a new section ‘Total SIF at the runoff limit’ presents Fig. 5.\n• Other small changes made with the intent of improving comprehension can be identified in the\ntracked-changes manuscript.\nLine numbers below refer to the tracked-changes manuscript.\nSlightly revised bottom-up SIF model outputs for field area\nDuring our revisions we noticed that the bottom-up SIF model outputs in our field area used a model\nrun initialised with non-isothermal snowpack due to the wrong start date (2022-07-05T18:00). It should\nhave started slightly later (2022-07-06T14:15) to match when the snowpack first became isothermal.\nThis caused very large initial heat fluxes which are unphysical within our model formulation.\nWe have re-run the model with the correct start date. This increases our summer-time bottom-up\nSIF at FS2 very slightly, from 0.140 to 0.146 m w.e (rounding to 0.15 m w.e). This is slightly closer to\nour measured average SIF of 0.16 m w.e. away from surface water features. The difference is\ncaused by refreezing over the first several days of the model run.\nReviewer #1 (Remarks to the Author):\nOverall, I find that the authors have adequately addressed my comments.\nThe inclusion of supplementary Figure 1 is appreciated, as this partially clarifies some of the concerns\nof reviewer #3.\nI also appreciate the incorporation of the recent Culberg et al. 2024 results of ice slab maps into the\nresubmitted manuscript.\nWe’re very pleased to receive these positive comments.\nI appreciate the removal of the statistical upscaling of remote sensing observations approach, as the\nlarge uncertainties, limited correlation between SIF and crosspol backscatter, and lack of validation\noutside of the study area made these results suspect. While I appreciate the new ‘bottom-up’ ‘top-\ndown’ approach for upscaling SIF estimates, I do have some minor questions about this new method.\n-Line 470: Could the authors briefly comment on to what degree the assumption of a constant slab\nthickness of 12 m has on the bottom-up SIF calculation via the 1-D heat flow model now adopted in the\nrevised manuscript? I appreciate that there is no existing estimate of slab thickness, but if we take as\nend members thicknesses of 5 m and 15 m (overestimate of 7 m vs. underestimate of 3 m), how\nsensitive are the bottom-up SIF estimates to changes in the slab thickness assumed?\nWe have made a new sensitivity analysis in which we run the model at our field site FS2 for slab\nthicknesses of 5 m and 15 m. This shows that the importance of the slab thickness to SIF is\nnegligible and can be disregarded. (L538-541).\n-line 570: An analogous comparison can be drawn to the ‘top-down’ case, where the sensitivity in\nsnowpack thickness is now assessed. It appears an attempt at this was done with thicknesses of 0.5,\n1.0, and 1.5, but this is not discussed any further. What were the results from this sensitivity analysis,\nand what snowpack thickness was ultimately assumed? The ‘Localised SIF continues into autumn’\nsection suggests that a thickness of 0.6 m was ultimately used, and a comparison was made to 0.3 m\nand 0.9 m, but this is seemingly contrary to what is described in the methods section, unless I\nmisunderstand. Clarification here would be helpful.\nLine 570 refers to the reference list, so we had to guess here. We think this comment concerns the\nchoices of snowpack thickness for top-down SIF simulations between (a) the field and (b) the ice-sheet\ncases. As a recap:\n• In the field area, we used 0.3, 0.6 and 0.9 m overlying snowpack thicknesses, as introduced in\n‘Localised SIF continues into Autumn’ and mentioned later in the Methods section ‘Modelling\nsuperimposed ice formation’. The sensitivities are presented in Fig. 4c. We use 0.6 m as our\ncentral estimate based on analysis of the thermistor string measurements at FS2.\n• Ice-slabs-wide, we used 0.5, 1.0 and 1.5 m overlying snowpack thicknesses, as introduced in\n‘SIF across Greenland’s ice slabs’ and then in more detail in the Methods section ‘Ice-slabs-\nwide wetness and SIF -> Bottom-up SIF’. The sensitivities are shown in Fig. 6d (shading around\neach line represents the range in SIF associated with each choice of snowpack).\nWe tried to identify changes to make these details clearer. We have added a line to the Methods\nsection ‘Modelling superimposed ice formation’ (L428-9):\nHere we describe the model formulation and its application to our field study area. Its ice-\nslabs-wide application is described separately in `Ice-slabs-wide wetness and\nsuperimposed ice formation', below.\nReviewer 2\nThe authors have put in a great deal of work to revise this manuscript and have resolved the majority of\nmy concerns. The new physical approach to upscaling their observations is a significant improvement\non the previous empirical scaling, and although it remains subject to large uncertainties, these are now\nmuch better quantified and explained to the reader. I also appreciated the attention to the written\nnarrative, which is now much clearer than in the original submission. I have one major comments on\nthe heat conduction models and a few minor comments for clarity that I describe below.\nWe are grateful to the reviewer for these positive comments.\nMajor Comments:\n[1] Overall, I think the authors make a reasonable argument for using a simple heat conduction model\nto estimate bottom-up and top-down SIF, given the many challenges in modeling slush and lateral\nwater transport in firn. However, I have a couple of questions/thoughts on this approach.\n• First, it seems that top-down refreezing may still be important in the summer, since the\nthermistor measurements in Figure 2b show a period of almost 15 days at the end of August\nwhere temperatures at and above the ice slab drop below 0C and the ice slab itself begins to\ncool again slightly. Currently, the modeling approach only considers bottom-up refreezing\nduring the summer. I understand that there is more uncertainty in top-down refreezing in the\nsummer due to uncertainties about the overlying snow thickness and its change with time.\nHowever, I would encourage the authors to think about whether there is a way that they can\nreasonably incorporate summer top-down refreezing (or quantify some uncertainty associated\nwith not including it) for their pan-Greenland estimates. I advocate this because the study\nencompasses some cold years where it is highly likely that the snowpack did not remain at 0C\nfor the entire melt season.\nWe agree that this issue is worthy of further consideration. We have undertaken a new sensitivity\nanalysis using a combination of PMR melt extent duration and our top-down SIF model running\nwith a fixed -10c surface temperature and 5% irreducible water content in the overlying snowpack.\nThis analysis indicates that while thin (0.5 m) snowpacks overlying the slush may result in some\nappreciable summer top-down SIF, potentially as much as 6-7 Gt over the 2017-2022 study period,\nemploying a thicker (1.0 m snowpack) reduces this value dramatically to less than 1 Gt. In relative\nterms, these findings indicate summer top-down SIF generates a maximum of ~10% of total SIF but\nlikely much less.\nGiven that top-down summer SIF is not a major contributor and in the interests of keeping the\nmain text a manageable length, we make a small note in the main text (L240-243). We have added\nthe full details of the analysis to the Supplementary Discussion (‘Top-down refreezing during\nsummer melt hiatuses’) along with a new Supplementary Figure 6.\nWe also note that here we considered for the first time the impact of irreducible water content in the\noverlying snowpack pore space, which does not contribute to SIF but must be refrozen before the\nsurface-originating cold wave can begin to form SIF in the slush matrix. This process also applies at\nthe start of autumn, so we modified our top-down autumn SIF runs to start with 5% irreducible\nwater content (in-between the classical 7% value of Colbeck (1974) versus the depth-dependent\nscheme of models such as MAR which use a range of 7 to 2% with increasing depth), rather than\n0% as previously. This reformulation has the impact of very slightly reducing top-down SIF\n(because it cannot commence until the overlying snowpack has frozen); this effect is visible in Fig.\n4c. We have updated the autumn model runs in both our field area and ice-slabs-wide, which\nresults in reduced autumn top-down SIF and therefore a slight reduction in slabs-wide total SIF.\n• Second, I think it would be more accurate to describe the results of this model as the “SIF\ncapacity” or “SIF potential” since this is a purely thermal model that assumes infinite and\ninstantaneous melt availability. Therefore, the mass of refrozen ice calculated by the model is\nreally an upper bound – the maximum that could be refrozen over the time period of interest\ngiven the initial thermal state at depth and the evolving air temperature. Melt availability would\nmodify the true amount of SI that forms. I think this is a relatively minor reframing that can help\navoid any overstating of the results.\nConcerning the model basis as thermal, assuming infinite/instant melt supply: We have added this\nsentiment using the suggested wording where we first introduce the model (L119-20).\n‘Melt availability would modify the true amount of SI that forms.’. We highlight that our slabs-wide\nresults already take this into account. We exclude any SIF which forms in excess of total seasonal melt\n(albeit according to MAR melt which as we highlight in the Discussion has its own limits). But in this\nsense our calculations are more realistic rather than simply maximum bounds, especially when also\nconsidering our use of van de Crux et al. (2024) T10m uncertainties to bracket our calculations.\nMore broadly, we are of the view that there remain substantial uncertainties and unconstrained\nparameters, especially in treatment of top-down SIF. These prevent us from being able to draw the\nconclusion that our central ice-slabs-wide SIF estimate is truly an indication of maximum SIF potential\n– it may be an over-estimate but it may alternatively be an under-estimate.\nMinor Comments:\nLine 81 – might note here that the cSMB is calculated from repeat firn cores so that the reader does not\nhave to go digging in the methods to figure out if you are using a climate model, SEB model,\nmeasurements, etc.\nAdded.\nFigure 3 – I am unsure what the question marks on the lateral boundaries of the snow, firn, and slush\nare intended to convey.\nAs noted in the caption, these question marks indicate unconstrained snowpack and slush matrix\nthicknesses. We have made some small adjustments to the figure (all question marks same colour,\nand added horizontal ticks between the question marks) which we hope make this clearer.\nFigure 4\n• What is the y axis on the histogram? I did not know how to interpret [-].\nThis was unitless on a 0-1 scale. We have multiplied to percentages and labelled accordingly.\n• Panel c, consider adding a legend or labels for the different green lines and the yellow line to the\nfigure itself. In the caption, consider mentioning the source of the temperature timeseries.\nWe have added a legend; the caption now indicates source of temperatures.\n• Either label or explain in the caption what A↑, A↓, and S↑ stand for. I figured it out", "reviewer_3": "Overall\nI have spent another day reviewing this paper which for me seems to be too much at this late stage.\nThe authors have chosen to drop key aspects of their original study, including the SAR based retrieval for superimposed ice\nformation (SIF) and also the timing work. They have instead added more modelling, and they make use of new external\ndatasets on ice slab extent and firn temperature. These are all good decisions, and the paper is clearer and more plausible for\nthe changes.\nI don’t think the authors have addressed the main concerns raised by the reviewers. In several cases they are simply\nrejected with arguments in place of evidence and I don’t think this is acceptable. I would have expected the rebuttal letter to\nsummarise the main concerns and how they have been addressed. Instead, the authors summarise what they believe to be\nthe most important issues instead of the major issues the reviewers identified.\nI still have major concerns with the study:\n1. Failure to adequately address reviewer concerns\nReviewer 2 notes that the small study area may not be representative of the ice sheet; this is not properly addressed.\nReviewer 2 notes that new methods should be validated; this is still not the case with the wetness derived from Sentinel-1\ndata.\nReviewer 3 notes that runoff is not unknown and yet the authors continue to insist it is.\nReviewer 3 notes that the authors experiment does not measure SIF directly, it measures the depth to the ice layer surface.\nThe authors do not address this.\nReviewer 3 notes that the correlation between melting and proximity of surface water is not clear. The authors do not\naddress this.\n2. There is no error estimate, meaning the findings are not confident.\nIn response to reviewer 1 the authors state that “the model does not come with uncertainty bounds because of the complex\nand nonlinear nature of the modelling.” This is not acceptable; an error estimate should be mandatory for publication.\nThe new data in fact lead to a revision of the retained SIF from 32 gigtaons to 60 gigatons between 2017 and 2022. No\nuncertainty estimate is given either in the original or revised manuscript, so it is not possible to say whether these two\nestimates are in agreement within experimental error. But this is a key take home message – the authors main conclusion\ncomes without any assessment of its confidence. It is also very obvious to this reviewer, but not the eventual readers, that\ntheir estimate changed by more than a factor 2 during the review process without any assessment of its reliability. This is for\nme a major concern, and the authors should I believe include an estimated error on this primary conclusion.\n3. The main finding, that SIF is only 2-3% of the runoff, is not stated.\nThe main finding – that 60 Gt of ice is retained – amounts to a 2 to 3% adjustment to the ice sheet runoff during the same\nperiod, well within the estimated 16% certainty of that measurement. I think it would be helpful for the authors to reflect this in\ntheir abstract.\n4. The paper is too verbose, poorly structured, and still has plenty of loose language.\nThe abstract contains more loose language. For example, it is stated that “However, the amount of runoff versus retention\nfrom these areas has not been measured.” but the runoff has certainly been measured, so why is this noted? If “these areas”\nmeans the ice slabs then say so, otherwise I read it to mean “the lower accumulation area” as the previous sentence states.\nThe phrases “at first order” and “most likely” are in there; what do these mean ? Why say “the majority y of superimposed ice\nforms during the melt season” instead of stating the estimated value? This is really quite elementary stuff and should be\npicked up before submission by authors.\nThe paper structure doesn’t distinguish adequately between the observational and modelling aspects of the work, which are\nquite distinct and should be clearly separated to the reader. Instead there are sub-headings related to processes and\nconclusions which is not helpful at all. I suggest laying the results out in a more tradiational manner; observations then\nmodels etc. There is also a tendency to draw firm conclusions from the model (also in the sub-headings) when I am not so\nsure that this is truly diagnostic as the models are primed to deliver results related to the proximity of water, for example.\nOverall I think too much of the paper is given over to the modelling when this is the least reliable aspect of the work,\nmore emphasis could be given to the observations.\n5. The in situ SIF data are not as strong as the authors would like us to believe.\nAlthough 3 sites close to rivers show high (0.5-0.6 mwe) SIF, one site also very close shows low (0.2 mwe) SIF. All are\nwithin 100 m of a feature and so there is really not enough information to draw a firm conclusion other than that the\nregional rate of SIF is about 0.2 to 0.3 mwe and that rates can be 2-3 times higher close to surface rivers. This is not the same\nas saying that “proximity to visible surface hydrological pathways strongly determines SIF magnitude”. I actually think the\nlateral water flow is the more interesting conclusion that could be drawn and explored from this dataset; one could set up a\nsimple model of a river surrounded by firn and assess the percolation. The authors state that “at least the majority of SIF can\nbe explained by direct refreezing of local surface melt” but the actual fraction can be estimated from the data; why not do this?\n6. The figures don’t showcase the main results.\nFig.1 should be redrawn. The location map (panel a) is the largest component while the field measurements of SIF (panels\nb, c, and d) are the smallest. This should be the other way around; the field measurements are the main data collected for\nthis study and should be prominent. Supp Fig 2 should be in Fig 1 also and the chart annotated with the site names shown in\nFig. 1 panels c,d,e; it’s probably the most valuable new information provided.\nSupp Fig 3 shows a comparison between the model and observed SIF. This figure should also be in the main text (again\nwith field sites labelled) as it is the main evidence to support use of the model for regional and ice sheet wide estimates.\nIts not great though, is it (r=0.79)? You could use the rms departure as a very rough estimate of the model certainty,\nthough the in situ measurements themselves are subject to uncertainty too. I guess its about 0.1 to 0.2 mwe/yr certainty and\nthat’s comparable to the background SIF of 0.16 mwe which would mean a ~100% error.\nFig 4 shows maps of the wetness duration (b) and the % of melt that refreezes (e). These are interesting but a map of the SIF\nwith the field observations from Fig 1bc,d would be far more valuable thing to show as it would show clearly how the model\ncompares to reality. The other panels in Fig. 4 could be moved to the supplementary materials.\nFig. 6. The components of SIF (panels b,c,d) can be moved to the supplementary materials; this would declutter the image\nand refocus it to just to the total SIF and its spatial / temporal variability.", "reviewer_author_response": "Response to reviews\nConcurrent superimposed ice formation and meltwater\nrunoff on Greenland’s ice slabs, version 1\nSummary\nWe thank the three reviewers for their engagement with our first submission. We believe that we have\nsuccessfully addressed their major concerns, leading to major improvements which warrant re-submission\nthrough the Nature Communications appeals process.\nReviewer 1 found the manuscript to be well-written, thorough, and provide important observational evidence.\nTheir concerns were essentially minor. The bulk of our changes concern comments by reviewers 2 and 3.\nReviewer 2 said that the manuscript addresses an important caveat to the popular notion that ice slabs lead to\nrunoff. They found the field measurements novel and convincing but were concerned that the remote sensing\nanalysis had poorly quantified but likely large uncertainties. Reviewer 3 had some major concerns about both\nthe field measurements and the upscaling.\nOur most major and important changes are:\n1. We have developed a new upscaling methodology. We no longer rely on statistical links between\nfield-measured SIF, negative degree days and SAR. Instead, we extend our physical heat flux\nmodelling approach to estimate SIF through (a) ‘bottom-up’ (i.e. slab cold content) refreezing during\nboth summer and autumn, and (b) ‘top-down’ refreezing (i.e. negative air temperatures) during\nautumn. We still use SAR sigma0HV, but only to determine wetness duration both across the field\narea and ice-slabs-wide through application of a binary threshold. We undertake these calculations in\nboth our field area and ice-slabs-wide. It also allows a much more rigorous exploration of\nuncertainties. In our view this is a more elegant and physical approach, particularly because it allows\nus to fully differentiate between summer and autumn SIF. Please see our headline response to\nreviewer 2 for more details.\n2. Our ice-slabs-wide analysis now incorporates the new ice slab extents mapped by Culberg et al.\n(2024); as these extents are much larger than Jullien et al. (2023), our estimate of SIF onto ice slabs\nhas increased accordingly.\n3. We have recomputed our SAR wetness threshold and we have examined the sensitivity of our\ncalculations to the choice of the SAR sigma0HV threshold for mapping wetness; our findings are\nrobust across the specified sensitivity range of the Sentinel-1 SAR sensor.\n4. We disagree with the view of reviewer 3 that our in-situ stake measurements of SI are insufficient to\ndemonstrate our purported findings. In response we have extended our description of our SIF\nmeasurement protocol, including the addition of a sketch figure to the supplementary materials. Please\nalso see the rebuttal in this response document. We do however agree that our measurements to\ncapture SIF timing are not fully convincing, so have removed these accordingly.\n5. We have made many improvements to the main text which we hope improve the readability.\nOur changes in bold\nIn our longer responses, the key parts are underlined.\nThis manuscript presents estimates of superimposed ice formation on top of ice slabs along the K-Transect in\nSW Greenland. In this study, the authors present field measurements that suggest local ice formation on ice\nslabs can retain ¾ of local annual surface melt. They estimate that ice formation along ice slabs facilitated the\nretention of ~32 gigatonnes of melt between 2017 and 2022, roughly ~4% of the Greenland SMB.\nOverall, the manuscript is well-written, thorough, and presents important observational evidence for surface\nmelt refreeze onto ice slabs, which plays a role on the overall SMB of the GrIS. I have several comments,\nmost of which are relatively minor; I believe that their consideration will overall improve the manuscript.\nI have organized the comments below, grouped by section within the manuscript.\nResults:\nPg 3 line 77: Because local melt is near uniform across the study area (presumably as a result of their\nmodeling results, but this is not made explicit), they hypothesis that measured spatial variation in SIF is\nindicative of inputs not only from local melt production, but also lateral water flow. This seems entirely\nreasonable given the observed increases in SIF for sites proximal to major surface rivers. However, the\ninclusion of external references to reinforce this claim would be beneficial.\nAs the reviewer infers, we take melt as a smoothly varying field of decreasing melt with increasing elevation,\nas calculated by a surface energy balance model forced by in-situ observations. Fig 5d now illustrates how\nmelt decreases with increasing elevation. We also now refer to findings from Noel et al. 2016 in the first\nparagraph of our Results. They showed that melt rates show very little spatial variability in this region.\nIt is possible that, once visible slush fields have developed, then local melt rates might increase further as a\nfunction of their reduced albedo. However, these slush fields can only initially develop by the transport and\nponding of melt that is essentially uniform in space considering the smooth surface topography over 100s\nmetres scales and homogeneous snow surface cover (e.g. Holmes 1955; Noel et al., 2016; Clerx et al., 2024/in\npress). In other words, their origin does not lie in albedo. Nonetheless, we now note this possibility of\nsubsequent local melt enhancement in the first paragraph of the Discussion.\nPage 4 line 79: “Finally, beyond 2–4 km away from visible surface runoff, we found that SIF began to\nincrease again, as illustrated by 0.25 m w.e. at site 3u3 (Suppl. Fig. 1).” – With what degree of confidence can\nthis claim be made, given that it is reliant upon a single point wise measurement of SIF, and the measured\nvalue is within 2 standard deviations of the reported average SIF from 2022?\nOn reflection, we agree, so we have removed this statement from the main text.\nPage 7 line 139: “In contrast, larger amounts of SIF can only be explained by lateral transport of meltwater\nwhich provides additional water for SIF to continue in specific locations beyond the end of the melt season.” –\nSimilar to first comment; Isn’t it possible that this elevated SIF can also be explained by spatial heterogeneity\nin melt that is not captured in the SMB model? I tend to agree with the authors that the more plausible\nexplanation is lateral transport of meltwater, but to what degree can spatial heterogeneity of melt be ruled out?\nPlease see our response to the previous comment about spatial variability in melt.\nAn important point to mention in the refreezing during Autumn section is the role that volume scattering due\nto surface and subsurface heterogeneities (e.g. surface roughness, sastrugi, ice lenses, ice pipes, percolation\npipes) plays. Implicitly, the bias-removal approach that is employed here *should* account for this process,\nprovided the geometry of these features is assumed temporally non-varying, such that changes in the cross-\npolarization backscatter can be attributed solely to changes in dielectric permittivity rather than surface\nroughness and volume scattering. I would suggest that this point be elaborated upon in the methods section,\n2\nperhaps with reference to some relevant literature (e.g., Rignot 1995, Backscatter Model for the Unusual\nRadar Properties of the Greenland Ice Sheet).\nThanks for this suggestion. We have elaborated in the Methods (‘Detecting residual water presence’)\naccordingly.\nMethods:\nPage 16, line 396: Regarding the model producing 0.25 m w.e. melt; is this merely a consequence of the\nsubsurface model discretization being 0.25 m? If not, what is the uncertainty of this model-derived estimate of\nw.e. melt?\nNo, this is not a consequence of the sub-surface model discretisation. In each time step, the model’s grid\nmoves up/down after the occurrence of accumulation/ablation; van As et al. (2017) contains full details.\nConcerning uncertainty: the model does not come with uncertainty bounds because of the complex and non-\nlinear nature of the modelling. This makes validation a crucial aspect of SEB modelling, as is done in\nnumerous previous studies. Concerning this model specifically, see e.g. van As et al. (2017) cited in\nmanuscript. In this manuscript we therefore evaluate how sensible its melt estimates are by reference to in-situ\nsurface lowering measurements, which show excellent agreement (see Methods).\nPage 18 line 469: Section on the linear regression between measured SIF and mean crosspol backscatter:\ngiven that an uncertainty in the in situ SIF values are provided, it would be worthwhile to propagate this\nuncertainty through the linear regression for a nominal uncertainty estimate of the linearly-regressed area-\nwide SIF estimates, particularly since the R2 in the linear regression isn’t particularly strong.\nThe linear regression approach has been removed in this re-submission.\nThis manuscript describes field measurements of superimposed ice formation on top of ice slabs on\nGreenland’s K-transect in the summers of 2021 and 2022. Using repeated firn cores, time domain\nreflectometry, and thermal modeling, the authors show that a significant portion of annual melt refreezes on\ntop of the ice slabs due to ice slab cold content and cold surface temperatures during autumn. In particular, the\nauthors show that refreezing occurs even where water flows laterally over the surface, potentially for kms,\nsuggesting that surface mass loss may not scale directly with visible surface runoff during the melt season.\nThe authors then use Sentinel-1 SAR measurements to establish relationships between their field\nmeasurements of SIF and radar-inferred surface wetness measurements. They use these relationships to scale\ntheir SIF measurements first to the area surrounding their field site and then to the entire ice sheet. The\nauthors conclude that SIF accounts for about 4% of the contemporary GrIS SMB, making it a small but non-\ntrivial contributor to the overall ice sheet mass balance.\nOverall, this paper addresses an important caveat to the popular notion that ice slabs lead to complete surface\nrunoff. In particular, the field measurements are novel and convincing and suggest that superimposed ice\nformation (SIF) needs to be accounted for in order to understand how quickly surface melting becomes mass\nloss as melting advances inland on the ice sheet. However, the remote sensing analyses used to scale these\nfield measurements to the entire ice sheet are subject to large uncertainties, many of which are not well\nquantified in the current manuscript. I am particularly concerned that the approach to estimating pan-\nGreenland SIF relies on a highly uncertain relationship between SIF and nCDD in a single ~30 km x 30 km\nregion which cannot be validated anywhere else on the ice sheet, as well a background SIF rate that is\nassumed to not vary at all across the ice sheet. This opens the potential for large, uncharacterizable regional\nbiases to affect the pan-Greenland SIF estimates.\n‘[the authors use a] background SIF rate that is assumed to not vary at all across the ice sheet’. In fact, the\nbackground SIF rate used for the ice-slabs-wide analysis varies in two ways:\n(1) As in our first submission, our treatment of the conditions which permit SIF across the ice sheet is\nexplicitly spatial. At a given MAR model grid node, SIF can only occur during the melt season and\nwhen the slab is covered in (porous) snow which provides the basis for a slush matrix overlying the\nslab.\n(2) Previously, we used three different fixed initial temperature scenarios to model bottom-up summer\nSIF. In this re-submission, we instead use the state-of-the-art ice-sheet-wide 10 m depth\ntemperatures data product published by Vandecrux et al. (2024, TC) to initialise our bottom-up\nSIF modelling. We also run our model at +/- the stated uncertainty of the Vandecrux product.\nWe consider this to be a major improvement compared to the previous version of our manuscript.\nConcerning the remote sensing analyses: we concur with the reviewer’s concerns that these were subject to\nlarge uncertainties and had the potential for large, uncharacterizable regional biases to affect the pan-GrIS SIF\nestimates. In conjunction with similar concerns raised by reviewer 3, we have therefore made major\nrevisions which replace the previous statistical approach completely with physical heat flow modelling\nto make our upscaling calculations. In summary:\n(1) In areas which remain wet after the end of the melt season according to SAR, we now allow our\nsummertime “bottom-up” calculations to continue into autumn until SAR indicates that the area is\ndry, yielding the additional bottom-up SIF which occurs in autumn. We partition the contributions of\nsummer vs autumn bottom-up SIF using the last day of melting from passive microwave radiometry.\n(2) In these same areas, we now also explicitly address “top-down” refreezing as a result of negative\nsurface temperatures, beginning once the melt season has finished. These calculations are very\nsensitive to the thickness of the snowpack overlying the slush matrix; we discuss this further in the\nmanuscript.\n(3) Any pixels suspected of being lakes (i.e. still ‘residually wet’ at the end of the year) are now\ncompletely excluded from our calculations.\n(4) We have recalculated the SAR sensitivity threshold using data acquired in all areas labelled as\n‘minimum likely’ slabs by Culberg et al. (2024, TC), with little resulting change (see below).\n4\n(5) Over our field area, we now show the uncertainty in our SIF calculations associated with the stated\nsensitivity of the Sentinel-1 SAR sensor used to determine wetness duration. We show that this source\nof uncertainty is second order, i.e. wetness duration does not change appreciably within this\nsensitivity range. We also checked this on an ice-slabs-wide basis where again it is a second-order\nsource of uncertainty (see Supplementary Figure 4).\nMajor Comments:\n[1] I am a bit confused about the primary purpose of modelling of the rate of superimposed ice formation. The\nauthors state that the reason they conduct the modelling is to determine the degree to which ice slab warming\nis due to latent heating. If this is indeed the main goal, then the argument would be better supported by\ncomparing the model run presented in the paper to a model run forced by the measured surface temperature\nand which only considers heat conduction (no latent heating). The difference between those two scenarios\nwould provide a stronger quantification of the role of latent heating in warming the ice slabs.\nHowever, the primary conclusion that seems to come out of this section is that refreezing due\nto ice slab cold content can only account for measured SIF rates in regions without additional\nadvection of meltwater. If the goal is instead to accurately model local SIF, then it seems like the authors\nmight be better served by using their fully coupled SEB-firn model instead of this simplified model that\nimposes a fixed meltwater supply that may or may not be consistent with the actual melt supply and timing\nduring the season when measurements were made.\nWe agree that we should have been clearer here. Our main goal was neither to super-accurately model local\nSIF nor to determine the degree to which ice slab warming is due to latent heating. Rather, it was to provide a\nconstraint on how much SIF can occur during the melt season under the plausible assumption of a continuous\nmeltwater supply at the slab surface. Then, SIF in excess of this logically occurs during autumn, i.e. after\nactive surface melt has finished. In other words, the primary purpose of this modelling approach is to partition\nthe SIF that we measured into the melt and post-melt seasons, thereby setting the scene for a slabs-wide first-\norder estimate of SIF later in the study. With our revised modelling approach we believe that this is\nmotivation is much clearer.\nAn ideal solution would indeed be to use a fully coupled SEB-firn model. However, no such model that is\nsuitable for our application exists. We are dealing with the problem of ponded meltwater in a snow matrix\natop an impermeable and cold barrier, which is an extreme end-member of what can be considered ‘firn’. We\nunderstand from firn modelling colleagues that even when ponding subroutines are available then they often\nremain disabled for want of validation. We now describe these issues in the Methods section ‘Modelling\nSIF’. Thus, it is more informative to simplify the problem as far as possible, i.e. by considering a relatively\nthin snowpack with a continuously wet slush matrix on top of an ice slab (conditions which were met for\nmuch of the melt season according to our in-situ measurements and further bolstered by observations in Clerx\net al., 2022) and then using classical heat flux theory to determine the rate of refreezing at first order.\nNew in this re-submission is an alternative configuration of the model to estimate SIF caused by\natmospheric cooling (‘top-down’). As we describe in the paper, an issue here is that we lack good bounds on\nthe thickness of the snowpack overlying the slush matrix. We wanted to use end-of-season snowpack\nthickness simulated by MAR but found that the outputs were highly spatially incoherent, which is a further\ncaution against 'blind’ use of RCMs assuming that they always deliver the best results. The reviewer may be\ninterested to see a paper with several of the same co-authors as this submission currently in Discussion, which\ndescribes some of the difficulties with RCM behaviour at the runoff limits (https://doi.org/10.5194/egusphere-\n2024-2750).\nOverall, we contend that our chosen approach is much simpler to interpret than the outputs of a fully coupled\nSEB-firn model. In turn, we hope that our findings can be used to advance ponding schemes in firn models.\nWe have made narrative changes to the relevant parts of the manuscript accordingly.\n5\n[2] Uncertainty in the regression between measured SIF and the multi-temporal mean 𝜎̂ 𝐻𝑉 can be better\nquantified and mapped into uncertainty in the study site SIF map (Fig. 4b). One option would be to use a\nBayesian regression that can quantify the range of linear models that are consistent with the data. This would\nmake it possible to produce an uncertainty map showing the range of the plausible SIF for each pixel in Fig.\n4b given the range of plausible models. But the key point is that there should be corresponding uncertainty\nmap for Fig 4b that account for uncertainty in the SIF-𝜎̂ 𝐻𝑉 relationship and the uncertainty in the background\nSIF rate.\nSimilarly, better uncertainty quantification is needed for the SIF-nCDD regression. Particularly given the\nrelatively poor R2 value for this regression, I strongly feel that it is important to quantify the uncertainty in the\nmodel fit and clearly demonstrate how that uncertainty maps onto uncertainty in the final estimates of regional\nautumn SIF that are presented in Figure 7.\nThis analysis has been removed from our re-submission.\n[3] For setting the SAR wetness backscatter threshold, I would like to see some sensitivity tests from different\nregions of the ice sheet showing how the 1st percentile of 𝜎̂ 𝐻𝑉 prior to the melt season varies regionally. For\na reliable pan-Greenland mapping, it is very important to establish whether conditions in SW Greenland are\nconsistent across the ice sheet. Other studies that have used Sentinel-1 data to quantify either melting or the\nextent of wet snow facies over Greenland have found that performance with a single threshold can vary\nsignificantly across the ice sheet (Hu et al., 2022; Li et al., 2023). Regional bias due to a fixed wetness\nthreshold should be quantified and incorporated into assessments of overall uncertainty.\nWe thank the reviewer for raising these concerns. First, we comment on the literature citations:\n(1) The findings of Hu et al. 2022 are not as relevant to our study as they first appear. Hu et al. did not\nnormalise summertime retrievals with respect to winter, unlike in our current approach, and so\ntherefore a considerable source of variability remains.\n(2) Li et al. 2024’s methodology is more similar to ours. They emphasise that SAR-based mapping of\nwetness at a set threshold is only reliable above the snowline, but that there is no appreciable regional\nvariability above the snowline. As we consider only snow-covered surfaces above the ELA, the\ndifficulties that Li et al. experienced with selecting a threshold applicable to melt detection over bare\nice is not relevant here.\nWe recap that our threshold of -2.1 dB was quite similar to Liang et al. 2021, who rather than the 1st percentile\ntake the less conservative 5th percentile of dry snow as their threshold: -2.66 dB. Both Liang et al., 2021 and\nLi et al. 2024 employ a more conservative value of -3 dB, +/- 0.3 dB, but they had the express aim of\ncapturing surface melting conditions which they validated by reference to in-situ air temperature\nmeasurements. This is not a relevant validation metric in our case as we are seeking to identify residual\nwetness after the end of active melt, which also indicates that a 3 dB threshold would be too conservative for\nour purposes here.\nNotwithstanding this comparison, we have undertaken another examination of our threshold choice. Now, we\ncompute histograms by sampling all areas within the ‘minimum likely ice slab’ classification of Culberg\net al. (2024), rather than just the western flank of the ice sheet as previously. The 1st percentile value\nhas decreased slightly to -1.17 dB (i.e. 0.07 dB more conservative than in our previous submission).\nIn the revised manuscript we adopt the new threshold of -1.17 dB, +- 0.3 dB indicative of the\nradiometric uncertainty of Sentinel-1’s sensor as per Geudtner and Torres (2012). We use the +/- to\npropagate an uncertainty in wetness duration through to our final calculations in our field area. On an\nice-slabs-wide basis, we provide a new Supplementary Figure which shows that the residual wetness\nduration/area is essentially insensitive to the radiometric uncertainty of the sensor.\n[4] I am concerned about how well the nCDD model can be generalized to map autumn SIF rates across the\nentire ice sheet. This is conceptually quite similar to the positive degree day model used for melting, which\nhas a spatially and temporally variable rate factor that typically needs to be tuned for individual sites. I have\nsimilar concerns for the nCDD model. Essentially, the authors’ scaling assumes that air temperature is the\nonly driver of refreezing and neglects all the other components of the subsurface energy balance, to include\nfirn and snow conditions, which might vary significantly across the ice sheet. It also assumes that the amount\nof melt that is refrozen scales directly with the length of time required to fully refreeze an area, which may be\ngenerally reasonable, but would absolutely be influenced by the transient accumulation rate, subsurface\ntemperatures, etc on the local scale. Without secondary validation sites, it seems basically impossible to know\nif the regression model coefficients will apply anywhere outside of southwest Greenland.\nThis analysis has been removed from our re-submission.\n[5] How are buried lakes accounted for when calculating autumn SIF? If I have understood the methods\ncorrectly, it seems like there has been no attempt to separate that mode of water\nimpoundment from superimposed ice formation. It seems problematic to apply a relationship\nbetween nCDD and refreezing in porous snow or a slush matrix to the refreezing of large open water bodies.\nSince buried lakes will have the highest and most persistent wetness values (see the authors’ own Figure 5,\nwhich is dominated by lakes), I am concerned that including lakes may significantly bias the results.\nThanks for raising this concern. Now, we screen out likely lakes by removing pixels which are still wet in\nthe last SAR acquisition each calendar year.\nMinor Comments:\nOverall, the section titles and introductory sentences could do a better job of explaining to the reader\nthe end goal of the analyses being described in the section. For example, the section “Modelling the rate of\nsuperimposed ice formation” would be easier to follow if it started with a sentence explaining why the authors\nwant to model this quantity. The section title could also better reflect the findings instead of the methods, for\nexample, something like “Ice slab cold content a key factor in summer SIF”.\nThanks for this suggestion. We have made several changes accordingly.\nLine 34 – I am not sure I understand what the second half of this sentence is trying to express.\nWe have rephrased as we agree that this was very unclear previously.\nLine 92 – Why was there no spike in 𝜖 when this water reached and fully flooded the probe before refreezing?\nWe cannot know for sure, but as the section ‘Timing of SIF’ has been removed in response to reviewer 3\nconcerns, we have not investigated further.\nLine 185/296 – there are now estimates of ice slab extent from both modeling (Brils et al., 2024) and Sentinel-\n1 observations (Culberg et al., 2024) that provide a more up-to-date and comprehensive outline of ice slab\nextent on the GrIS.\nIndeed - though we find it a little concerning that literature published in 2023 can already be considered less\n‘up-to-date’! Nonetheless, given that we have made major changes to our upscaling approach, we have taken\nthe opportunity to replace the Jullien et al. (2023) outlines with those based on Sentinel-1 observations\nby Culberg et al. (2024). We now explore the impact that the uncertainty in ice slab extent makes on our final\nestimates of SIF.\nWe have moved the statements of ice slab extent previously found at this location into the Introduction.\nwould at least provide a constraint on the shape and range of temperature profiles in ice slab regions.\nThanks for this suggestion. Rather than using modelled temperature profiles for initialisation with all the\nuncertainties in RCM sub-surface schemes, we take advantage of another study that has been published this\nyear: Vandecrux et al. (2024, TC), which mapped observed 10 m depth temperatures ice-sheet-wide using\nmachine learning techniques, presented along with an associated uncertainty product. We use the Vandecrux\ndata to set the bottom temperature of our initial temperature profiles and then scale our observed\nprofile from FS2 accordingly. We also run our calculations +/- the stated 10 m temperature uncertainty.\nPlease see our revised Methods section for more information.\nHow is the long-term equilibrium line in Figure 5 defined?\nThis was based on RACMO data 1981-2000. But in this re-submission we have modified this figure to only\nshow wetness in the Culberg et al. (2024) maximum ice slabs extent area.\nLine 297 – I found the context for the discussion future warming and peak refreezing a bit confusing. Is the\npoint simply that we should not project current SIF rates as a function of SMB into the future due to the\nuncertainty in climate forcing and hypsometry?\nWe have made efforts to improve the narrative of this paragraph – please see full text.\nLine 476 – where did n=4 come from when only 3 PROMICE weather stations were noted in the study area?\nWhich stations were used and how are they distributed in elevation across the study area? (e.g. is it fully\nrepresentative of the elevation distribution, or does this scaling do significant extrapolation outside the data\nbounds?)\nThe middle and upper weather stations largely encompass our study area, with the upper station being situated\nat 1840 m ASL. Above this elevation SEB model input variables are extrapolated, not interpolated, and\ntherefore calculations do not exceed 2000 m ASL so that extrapolation does not become unrealistic. The\nmodel uses these observations to compute in discrete 100 m elevation bins. We computed our regression\nacross four elevation bins (1,600, 1700, 1800 and 1,900 m asl). This was not clear previously, so thanks for\ndrawing this to our attention and we have updated the text accordingly.\nLine 506 – What determines if a cell becomes snow free? Is this from MAR, an albedo remote sensing\nproduct, or something else?\nAccording to modelled MAR snowpack. This is now noted in the text.\nLine 526 – wouldn’t a brief rise above 0C potentially lead to melting that could supply more water for SIF?\nWe have added a note that, given several weeks of prior atmospheric cooling of the snowpack, it is probable\nthat this melt would simply have been refrozen in the uppermost centimetres below the surface.\nSummary\nThis paper uses a combination of field observations, satellite observations, and simple/conceptual modelling\nto address the important topic of meltwater retention in the percolation zone of the Greenland ice sheet –\nabout which there is some uncertainty. The paper builds on the finding that meltwater retention has led to the\nformation of ice slabs at depth the firn pack, and that the presence of these slabs disrupts the normal process of\nmeltwater percolation and redistribution.\nThe paper is quite challenging to comprehend as it contains several elements of work that are not easily linked\nand are each quite obscure. Despite several attempts to read and review, I could not form a clear narrative\nfrom the study and, more importantly, I was not persuaded by the arguments presented. Either the narrative\nneeds to be strengthened, or the arguments need to be made more compelling with further work.\nMajor\nI have identified five major concerns that I hope will act as a guide for the authors to reformulate their paper.\n1. Loose language. There is lots of loose language and (for me) too much hyperbole in this paper. The term\n“runoff” in particular is central to this work but is not well or consistently defined; it is interchangeably used\nto describe meltwater, mobile water, and water that has left the ice sheet. In my book it is the latter only.\nElsewhere the term “rivers” is used but this is informal and not helpful; it would be better to make an\nassessment of the intensity of meltwater flow in supraglacial channels and contrast the magnitude.\nAnd it is stated in the abstract and repeated elsewhere that the amount of runoff from the percolation zone is\nnot known when this is simply not true – the authors used for example regional climate models that make an\nestimate of this and there have been many studies on meltwater retention including the identification of ice\nslabs on which this work is predicated. I would recommend a thorough parse of loose and hyperbolic\nlanguage.\nWe have endeavoured to tighten up our use of language throughout the manuscript, particularly around\nthe term ‘runoff’ which we agree was not used as precisely as it should have been. We hope that these changes\nmeet with approval.\nThe reviewer said that they find ‘too much hyperbole’. The only example they provided is: “it is stated in the\nabstract and repeated elsewhere that the amount of runoff from the percolation zone is not known when this is\nsimply not true – the authors used for example regional climate models that make an estimate of this and there\nhave been many studies on meltwater retention including the identification of ice slabs on which this work is\npredicated.”. Concerning this specific example, we disagree: regional climate models make estimates which\nso far remain unvalidated for the percolation zone, especially where ice slabs are present. They are potentially\nsubject to major issues associated with invalid parameterisations. We currently have a manuscript co-written\nwith the MAR and RACMO teams under discussion which covers many of these themes in detail:\nhttps://doi.org/10.5194/egusphere-2024-2750. Nonetheless, to tighten the language we are happy to change\nfrom ‘unknown’ to ‘unmeasured’, ‘unvalidated’. We also note that we have added further justification\nfor our heat flux modelling approach (over RCMs) in the Methods section ‘Modelling SIF’.\nThe reviewer states that there have been ‘many studies on meltwater retention including the identification of\nice slabs:\n• The reviewer does not provide examples of these ‘many studies’.\n• Taking this comment literally, we think the reviewer is suggesting that there are already several\nstudies of meltwater retention on top of ice slabs. We are not aware that such studies exist, with the\nexceptions of already-cited studies of SIF atop glacier ice in the Canadian Arctic and Svalbard. In\nJullien et al. (2023), Culberg et al. (2024), which do not simultaneously assess meltwater retention\nprocesses.\n• There is arguably perhaps one study, Machguth et al., (2016), which focuses on meltwater retention\non a slab. In this study, the term ‘superimposed ice formation’ appears only once, in the 3rd-last\nparagraph of the Methods, only in reference to core alignment. This study undertook ‘barcode\nalignment’ on a couple of firn cores at KAN_U, which as this reviewer points out subsequently can be\nsubject to substantial uncertainties. But having re-read this study in detail for this response we have\nconcluded that the evidence of SIF they present is so limited that we now do not cite it in the\nIntroduction as an observation of Greenland SIF.\n2. Superimposed ice formation. The initial part of this work is based on in situ measurements that purport to\ndocument superimposed ice formation. I am not persuaded that the authors have demonstrated this effectively.\nThe main approach is to measure the depth to the slab surface two years apart at a variety of locations. This is\nnot an actual measurement of superimposed ice formation, it’s a crude approximation. The slabs could be\nmobile, because the site is not at the equilibrium line. More importantly, the authors are not able to resurvey\nthe same site and so their measurement is in fact a combination of spatial and temporal variation in slab\nsurface depth. Without knowing the former, it is impossible to deduce the latter. This calculation\nunfortunately underpins the entire study, and so it is quite fundamental to convince the reader that the\nmeasurements are robust. Having drilled many firn cores in this region I know that measuring the depth to a\ncommon horizon is a very challenging thing to do with precision, and that there are typically large natural\nvariations in the depth of common horizons in any case. I wonder if repeat GPR would have been a better way\nto do this entire experiment. Or drilling through the entire slab. Or analysing the core chemistry to isolate\nannual layers. Certainly the spatial variation in slab surface depth needs to be reliable quantified somehow\nbefore the temporal variation can be deduced.\nWe strongly disagree with the thrust of this comment that our in-situ observations are invalid. We highlight\nthat neither of the other two reviewers raised concerns of this nature. From our perspective there seems to\nhave been a misunderstanding of our field measurement protocol. At root, this seems to relate to the locations\nat which we make our measurements. We do resurvey the same site between years: these locations are exactly\nthe same each year with respect to the underlying stratigraphy, as we returned to a location marked with a\nbamboo stake. Each location had of course moved by several 10s of metres due to ice flow, but this process\ndoes not affect the near surface stratigraphy. To clarify this further, we have improved our explanation,\nincluding adding a supplementary figure which graphically illustrates the measurement process, and\nsupplementary discussion to discuss local variability in SIF more extensively.\nIn our judgement the suggestion of using GPR is a non-starter. If we were to make a GPR survey, we would\nmeasure the thickness of the snow layer, albeit with larger uncertainty than by drilling. We ask ourselves: how\ndoes GPR measure the thickness of the slab and in particular its annual growth? We could potentially detect\nthe bottom depth of the slab in subsequent years, but would be hampered by all the GPR-related uncertainties\nthat would reduce the precision of the outcome in such a study. And crucially, even if a change in thickness\ncould be measured beyond uncertainty, we still would not know whether the slab grew from accretion below\nor above. In contrast, this is precisely what our approach of returning to a marked surface location in year 2\nenables us to do.\nThe reviewer states “Having drilled so many firn cores in this region I know that measuring the depth to a\ncommon horizon is very challenging…”. We have no way of assessing what expertise this reviewer has, nor\nthe scope of the region they are considering. The reviewer may have drilled firn cores at high elevations, but\nperhaps not into ice slabs. Our cores are from the slabs, so are basically solid ice. In these cores, to say where\nsnow ends and firn starts is not always clear, but the horizon where the slab begins is always abundantly clear.\nFor more details, the core logs can be examined in our data repository.\n3. Spatial variation in SIF. A main conclusion of this work is that the degree of superimposed ice formation is\nhigher closer to surface meltwater channels than it is distal from them. This is I think the most original aspect\nof the study. However, the conclusion itself depends on how reliable the estimates of SIF (see point 2), but\nalso on how strong the spatial correlation is. Looking at the data in Figure 1 I am not convinced that the\nspatial variation in SIF is significantly correlated with the distribution of meltwater channels; there seem to be\njust as many sites that buck the trend as there are that fit the hypothesis. The sites sampled are also quite\nrandomly distributed, which doesn’t help. Wouldn’t it make sense to survey a profile tangential to a channel,\nfor example? I think more work and quite possibly more data are needed to convince the reader that this\nconclusion is robust.\nAs the reviewer notes, “the conclusion itself depends on how reliable the estimates of SIF (see point 2)” are.\nAs we respond above, we are confident in our SIF estimates and believe that there has been a\nmisunderstanding here.\nConcerning spatial variation, Fig. 1 provides a high-level overview of the SIF magnitudes. Further evidence is\npresented in Supplementary Figure 1, which explicitly compares the quantity of SIF with the Euclidean\ndistance to surface drainage visible on optical satellite imagery. Our measurements were concentrated on two\ntransects, with a spacing that reflects the proximity to surface drainage. Sites FS1, 122, camp, 231, 3u3\nconstitute one transect. The sites shown in Fig. 1d constitute another transect, which is quasi-tangential to the\nmain surface channel. We now note these transect configurations in the main text. Concerning a profile\ntruly tangential to a channel: the location of the channel varies somewhat between years, so this is less feasible\nto deploy than it might first appear, and also what a transect equidistant of a stream is one year, would not be\nthe next. Concerning the measurement stations, these were installed to acquire data in three distinctively\ndifferent surface slope regimes and with contrasting surface hydrology in any given melt summer.\nPlease see also our response to reviewer 1 concerning possible spatial variability in melt.\n4. Timing of SIF. There is a section of work on the timing of superimposed ice formation deduced from measurements of the\nfirn dielectric properties. I found this section very difficult to navigate and not persuasive; the only thing that\nseems to be consistent about the probe measurements is how inconsistent they are in time and space. It seems\nthat very few measurements were actually made and yet apparently firm conclusions are drawn from them.\nThis is very similar to point 3; it suggests to me that the field experiment was again not up to the task and\nin hindsight other measurements would have been made. It’s also not clear to me how important it is to know\nthe timing of superimposed ice formation – is that really material to this study ? My feeling is this whole is a distraction and could easily be dropped.\nOn reflection, we agree that knowing the sub-seasonal timing of SIF is not important to our key messages and\nis challenging to interpret from our data. We have therefore removed this whole section as suggested.\n5. Modelling. The study next leaps into a series of “model” experiments designed to help scale up the field\nobservations. However, these experiments are tangentially linked to the field measurements at best. The first\n“model” is a physical 1D heat transfer model, the second is a linear regression between satellite radar\nbackscatter and measured superimposed ice formation, and the third is another linear regression between the\nmodelled superimposed ice formation and cumulative negative degree days derived (I think) from the 1D heat\ntransfer model. The heat flow model is very loosely related to the field observations as it is contrasted to\ntemperature profiles that are not directly related to the SIF. Moreover, we can’t really judge how well the\nmodel and observations agree from the data as presented.\nWe are grateful for these comments, which together with similar comments by reviewer 2 caused us to\nfundamentally rethink how we conduct our upscaling. We believe that we have now successfully\naddressed this concern and that our response to reviewer 2 answers this major comment.\nThe satellite radar backscatter is very loosely related to in situ superimposed ice formation; many other factors\naffect backscatter and in any case the correlation deduced is largely owing to the clustering of the\nmeasurements into two distinct populations, without which there would be no correlation.\nWe recognise that in our previous submission there was the potential to conflate radar backscatter with in-situ\nSIF, primarily because of the statistical way in which we attempted to conduct our upscaling; however, this\nwas not our intention. To be clear, we consider backscatter variations in autumn to tell us only about residual\nwetness in the near-surface. It is the disappearance of this backscatter that we attribute to SIF, which we now\nmodel physically instead of our previous problematic statistical approach. In the revised manuscript this is\nnow clearer because we talk only of residual wetness [duration] and not also in terms of sigmaHV0\nbackscatter strength like previously.\nThe degree of correlation between satellite derived SIF and negative degree days is relatively weak in any\ncase, half of the variation is not explained. I find this approach to be quite chaotic, highly empirical, and not\npersuasive. I would recommend the authors rethink their approach here entirely and find a better way to add\ncontext and value to their sparse field measurements. The authors do note that regional climate models with\ncoupled firn process models are another way to study meltwater retention and I agree this would be a better\nway to approach this problem.\nAs explained in our response summary, we have taken heed of these concerns about the upscaling\napproach and, we believe, now have a much better methodology with which to deliver our upscaled\ncalculations.\nConcerning the suggestion here that RCMs with coupled firn process models would be a ‘better way’ to\napproach the problem: we disagree, and also note that we do not state that they constitute a ‘better way’. As\nwe have tried to clarify in our re-submission, there are substantial uncertainties/problems remaining with how\nRCMs and even observation-based SEB models represent the sub-surface, which makes it desirable to instead\nadopt a simplified approach to modelling SIF. Please see in particular the Methods section ‘Modelling\nSIF’.\nMinor comments\nLine 1: I’m not sure the word rivers is helpful here as they are generally associated with land and rain fed. can\nyou describe in glaciological terminology instead?\nThe term ‘rivers’ is used elsewhere in the literature. Non exhaustive examples include line 2 of Yang and\nSmith (2016, JGR); Smith et al. (2015, PNAS). We have decided to stick with this intuitive terminology rather\nthan ‘supraglacial streams and rivers’ at every mention.\nLine 4: I don't think its unknown; there are estimates. need to be a bit clearer here what you mean\nWe have changed to ‘unmeasured’.\nLine5: should probably be clear about what you mean by runoff; is it meltwater until it is mobile, or until it is\nknown that it will leave the ice sheet locus, or until it has left the ice sheet locus?\nHere, in-line with the definition of runoff, water that is assumed to leave the ice sheet locus – updated to\nrefer to ‘visible runoff’ which is then linked to previous studies in the Introduction.\nLine 6: I think you mean that up to 3/4 of the annual melt can be retained by SIF (i.e. the converse of what is\nstated), or that the process of SIF can lead to meltwater being retained\nWe are a confused by this comment – the text said that ‘SIF….can retain up to roughly three quarters of\nlocally-produced annual melt’, which doesn’t seem the be converse. Nonetheless, we have tried to clarify in\nour revised Abstract.\nLine 7: Is the heat flow modelling informed somehow by the field observations? or how else are they linked?\nThis comment refers to the abstract. We have added a note that the modelling is initialised by in-situ field\nobservations. The full information is available in the main text as before.\nWe have added a note that they are used to derive residual wetness.\nLine 8: Quantify \"majority\"\nWe prefer not to do this in the abstract, as while the numbers are always a majority (i.e. greater than half),\nthe\nprecise value depends on parameter choices. The reader should read the main text for these details.\nLine 8: Quantify \"remainder\"\nSee previous comment.\nLine 9: This probably answers my earlier point about the definition of runoff\nYes, we believe so; see also modified abstract.\nLine 10: Using your model?\nWe are not sure what this comment refers to, sorry.\nLine 11: Does this mean total or local SMB? not sure SMB is a good metric since it is the sum of\naccumulation and ablation\nThis is a good point – we concur that SMB was a poor choice of metric. We now instead express this in\nterms of runoff estimates made using Cryosat-2 altimetry measurements. We estimate that recent SIF\nhas decreased recent runoff by 2-3 %.\nLine 12: Is there an uncertainty on 32 Gt?\nThis value has changed with our new upscaling approach; we now provide a ‘most-likely’ value (60 Gt)\ntogether with a range (27-75 Gt). The range is primarily a function of which ice slab extent we examine,\nbut\nalso varies depending on 10 m firn temperatures and autumn snowpack thickness.\nLine 12: To call them hotspots this is relative to other areas, but it doesn’t seem that you have assessed those?\nConsidering refreezing in absolute terms (metres water equivalent), they are indeed hotspots compared to the\nrest of the ice sheet. Nowhere else is this combination of (a) meltwater supply and (b) a snow/firn matrix in\nwhich to refreeze onto the underlying ice met.\nLine 19: Decreased not densified\nChanged\nLine 25: Is this really true? i think there are estimates in the models you are using- what I guess you are saying\nis that they don't account for potential retention associated with the ice slabs. Thats very different to unknown\nWe have replaced with ‘unmeasured’.\nLine 37: Running water not runoff\nChanged\nLine 40: Radar returns is not a jargon. Measurements?\npreviously (see for example MacFerrin et al. 2019 Supplementary).\nLine 42: Technically the slabs themselves are superimposed ice and so yours are not the first observations\nWe are observing superimposed ice formation, i.e. the process of forming superimposed ice, as we state in the\ntext. Slab mapping has indeed already mapped superimposed ice.\nLine 43: What is the runoff limit?\nIn lines 23-24 we now define the visible runoff limit alongside the visible runoff area.\nLine 44: What is local runoff?\nChanged to ‘net water discharge from the hydrological catchment’.\nFig 1: Why is fs1 high on both years?\nThat there was no slush field at the surface doesn’t mean that there was no slush beneath the surface. Visible\nhydrological features do not tell the whole story, so we also need to look at the sub-surface. C-band SAR\nclearly shows that this area remained wet for a long time after the end of the melt season (Fig. 5a, see wetter\nfeature between FS2 and downstream later, between 1730 and 1760 m asl).\nLine 69: You are showing that the SIF varies by 10s of cm w e over short distances, so perhaps the smb does\ntoo, let alone surface features? In which case you experiment is not well constrained\nPlease see our responses to reviewer 1 concerning surface melt rates.\nLine 88: Why use e>3 when water is e=80?\nThe TDR measurements have been removed from this re-submission.\nLine 94: TDR probe\nThe TDR measurements have been removed from this re-submission.\nFig.3 : It’s very difficult in this figure to contrast the observed and modelled temperature profiles, which seem\nto be the only validation of the model. I suggest replotting them as e.g. a direct scatter plot with a colour\ncoding to illustrate depth so that we can see how well they compare (or otherwise)\nWe have redrawn this figure panel to make the comparisons between observed vs modelled temperature\nprofiles clearer. We have also added some text to the Results which addresses the offset between\nobserved and modelled profiles at the end of the melt season.\nLine 139: Isn't this obviously the case? What else would be causing the SIF. I don't think this is a grand\nconclusion\nWe are not exactly sure what `this’ refers to, as one sentence finishes and another starts at this line number.\nAssuming that ‘this’ refers to the sentence in L137-139, we do not make out that this is a grand conclusion.\nRather, this statement supported the sentence which followed, i.e. ‘In contrast, larger amount of SIF can only\nbe explained by lateral transport of meltwater….’. We have made various minor changes to better explain\nour model findings.\nLine 141: Just because your model agrees with SIF in some places doesn't mean you can conclude that\n14\nsomething else must be happening elsewhere – that’s not a hypothesis you have tested by experiment. Is your\nmodel able to simulate lateral movement of meltwater? What other processes might lead to local variations in\nmelting that are not simulated by your model? More work required here\nThe majority of cores show rather similar SIF. We clearly show that simulated background refreezing agrees\nto SIF in several locations. It clearly does not agree in a few locations. We therefore looked for the reasons for\nthe disagreement. It is obvious that refreezing from only the melting directly above a given location cannot\nvary so enormously in space given how homogeneous the snow-covered surface is and the rather small\nelevation envelope that we worked within (see also our response to reviewer 1). We also observe the largest\nrefreezing rates where there is most slush. To us, it is logical that these locations should therefore have less\ncold content in the slab. Finally, it is unclear to us what the benefit of including lateral movement in our 1-D\nmodel would be.\nLine 144: It does not indicate this; it just indicates you have a match within 1sigma in some locations\nThis line says ‘Since our model indicates that the localised larger quantities of superimposed ice cannot be\naccommodated by slab cold content during the melt season alone…’.\nOur model does indicate this as the modelled summertime SIF is valid for most of our measurement sites, not\njust the one where we compare. Of course, we do not have more firn temperature data with which to undertake\na more thorough comparison. But in the context of how much larger some of SI measurements are (by\ndecimetres), then we consider a match within 1 sigma (0.06 m w.e.) sufficient for our assertions.\nLine 156: Correlation, not relationship\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nFig. 4: These data are clustered in 2 populations which is leading to most of the correlation. I don't think they\nare distributed well enough to develop an empirical parameterisation. why are they clustered in this way? I\nsuggest this is explored in more detail before jumping to conclusions\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nLine 174: This is all a bit contrived; your result is entirely dependent upon the relationship that you force to\nexist.\nThis empirical parameterisation has been removed in this re-submission. Please see our response summary.\nFig. 6: This isn't \"wetness\", its sigma nought. I think you are jumping to too many conclusions here\nAs we describe in the Methods, we derive a binary wetness indication by applying a threshold to detrended\nsigma0HV. Thus, Figure 6 did present wetness (duration), not sigma0. This approach is based on\npreviously published work, most notably Liang et al. (2021). Note that Fig. 6 is now Suppl. Fig. 5.\nLine 211: So only half of the variance is explained\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 219: I'm lost here; if you have a model of the firn pack why are you counting cNDD to model SIF? why\nnot just compute it from your actual model?\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nWe now summarise the specific issues (difficulty locating the slab surface at FS3; uncertainty related to\ndrilling into SIF until metal chips appeared).\nLine 336: Why are they minimum estimates?\nPlease see response to previous comment.\nLine 343: This is a poorly constrained experiment as you don't know the spatial variation in firn stratigraphy.\nhow can you be sure you are measuring superimposed ice?\nWe believe that this concern is addressed by our more detailed explanation of our field measurement\nprotocol and Supplementary Discussion; specifically, during our repeat measurements, we found very little\nspatial variability at metres scales.\nLine 358: So this is an underestimtae of the SIF?\nThis paragraph is about the cSMB, not about SIF. By definition (Cogley et al. 2011), SIF refers only to\nrefreezing on top of impermeable horizons, in this case the ice slab.\nLine 459: Why choose the 1st percentile? is this based on theory?\nPlease see our response to Reviewer 2.\nLine 461: Why likely?\nThanks for pointing this out – on reflection, we do not think we can quantify this statement that the threshold\nis ‘likely conservative’, so have removed it in this re-submission.\nLine 466: Extrapolate or scale up, not distribute\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 470: This isn't determining the existence of a relationship, it is assuming the existence of a relationship\nThis comment relates to the empirical parameterisation that has been removed in this re-submission.\nLine 486: Detect not retrieve\nChanged.\nConcurrent superimposed ice formation and meltwater runoff on\nGreenland’s ice slabs, version 2\nWe are grateful to the reviewers for their time engaging with our revised submission. In response we\nhave made further substantial changes which we believe have improved the manuscript.\nMain changes:\n• We have established that bottom-up SIF is not sensitive to slab thickness within observation-\nbased bounds (R#1, R#2).\n• We have added a sensitivity analysis to the supplementary materials about the potential\nimportance of summer top-down SIF – it is not large in comparison to other SIF sources and so\nwe can neglect it (R#2).\n• While considering summer top-down SIF, we realised that we should also incorporate refreezing\nof irreducible water content in the overlying snowpack in autumn. We have re-run our\ncalculations for both the field study area and slabs-wide with this process included. This\nreduces autumn top-down SIF and in turn slightly reduces total SIF compared to our previous\nsubmission. See our response to R#2.1 for more information.\n• We have improved several Figures concerning the field study (R#3, R#2). In particular, Fig. 1 is\nclarified based on R#3’s suggestions and now incorporates distance to drainage. Fig 4 has\nimproved labelling and now focuses only on autumn SIF. We provide a new Fig. 5 showing total\nSIF across the study area and including observed vs. modelled SIF.\n• When examining SIF at each stake, we now consider slope angle in addition to distance to\nsurface drainage (R#2, R#3). See Fig. 1d and the rewritten Results section ‘Superimposed ice\nformation depends on near-surface water flow’.\n• For our field study area we now dedicate “Localised SIF continues into autumn’ to presenting\nFig. 4, while a new section ‘Total SIF at the runoff limit’ presents Fig. 5.\n• Other small changes made with the intent of improving comprehension can be identified in the\ntracked-changes manuscript.\nLine numbers below refer to the tracked-changes manuscript.\nSlightly revised bottom-up SIF model outputs for field area\nDuring our revisions we noticed that the bottom-up SIF model outputs in our field area used a model\nrun initialised with non-isothermal snowpack due to the wrong start date (2022-07-05T18:00). It should\nhave started slightly later (2022-07-06T14:15) to match when the snowpack first became isothermal.\nThis caused very large initial heat fluxes which are unphysical within our model formulation.\nWe have re-run the model with the correct start date. This increases our summer-time bottom-up\nSIF at FS2 very slightly, from 0.140 to 0.146 m w.e (rounding to 0.15 m w.e). This is slightly closer to\nour measured average SIF of 0.16 m w.e. away from surface water features. The difference is\ncaused by refreezing over the first several days of the model run.\nReviewer #1 (Remarks to the Author):\nOverall, I find that the authors have adequately addressed my comments.\nThe inclusion of supplementary Figure 1 is appreciated, as this partially clarifies some of the concerns\nof reviewer #3.\nI also appreciate the incorporation of the recent Culberg et al. 2024 results of ice slab maps into the\nresubmitted manuscript.\nWe’re very pleased to receive these positive comments.\nI appreciate the removal of the statistical upscaling of remote sensing observations approach, as the\nlarge uncertainties, limited correlation between SIF and crosspol backscatter, and lack of validation\noutside of the study area made these results suspect. While I appreciate the new ‘bottom-up’ ‘top-\ndown’ approach for upscaling SIF estimates, I do have some minor questions about this new method.\n-Line 470: Could the authors briefly comment on to what degree the assumption of a constant slab\nthickness of 12 m has on the bottom-up SIF calculation via the 1-D heat flow model now adopted in the\nrevised manuscript? I appreciate that there is no existing estimate of slab thickness, but if we take as\nend members thicknesses of 5 m and 15 m (overestimate of 7 m vs. underestimate of 3 m), how\nsensitive are the bottom-up SIF estimates to changes in the slab thickness assumed?\nWe have made a new sensitivity analysis in which we run the model at our field site FS2 for slab\nthicknesses of 5 m and 15 m. This shows that the importance of the slab thickness to SIF is\nnegligible and can be disregarded. (L538-541).\n-line 570: An analogous comparison can be drawn to the ‘top-down’ case, where the sensitivity in\nsnowpack thickness is now assessed. It appears an attempt at this was done with thicknesses of 0.5,\n1.0, and 1.5, but this is not discussed any further. What were the results from this sensitivity analysis,\nand what snowpack thickness was ultimately assumed? The ‘Localised SIF continues into autumn’\nsection suggests that a thickness of 0.6 m was ultimately used, and a comparison was\nmade to 0.3 m and 0.9 m, but this is seemingly contrary to what is described in the methods section,\nunless I misunderstand. Clarification here would be helpful.\nLine 570 refers to the reference list, so we had to guess here. We think this comment concerns the\nchoices of snowpack thickness for top-down SIF simulations between (a) the field and (b) the ice-sheet\ncases. As a recap:\n• In the field area, we used 0.3, 0.6 and 0.9 m overlying snowpack thicknesses, as introduced in\n‘Localised SIF continues into Autumn’ and mentioned later in the Methods section ‘Modelling\nsuperimposed ice formation’. The sensitivities are presented in Fig. 4c. We use 0.6 m as our\ncentral estimate based on analysis of the thermistor string measurements at FS2.\n• Ice-slabs-wide, we used 0.5, 1.0 and 1.5 m overlying snowpack thicknesses, as introduced in\n‘SIF across Greenland’s ice slabs’ and then in more detail in the Methods section ‘Ice-slabs-\nwide wetness and SIF -> Bottom-up SIF’. The sensitivities are shown in Fig. 6d (shading around\neach line represents the range in SIF associated with each choice of snowpack).\nWe tried to identify changes to make these details clearer. We have added a line to the Methods\nsection ‘Modelling superimposed ice formation’ (L428-9):\nHere we describe the model formulation and its application to our field study area. Its ice-\nslabs-wide application is described separately in `Ice-slabs-wide wetness and\nsuperimposed ice formation', below.\nReviewer 2\nThe authors have put in a great deal of work to revise this manuscript and have resolved the majority of\nmy concerns. The new physical approach to upscaling their observations is a significant improvement\non the previous empirical scaling, and although it remains subject to large uncertainties, these are now\nmuch better quantified and explained to the reader. I also appreciated the attention to the written\nnarrative, which is now much clearer than in the original submission. I have one major comments on\nthe heat conduction models and a few minor comments for clarity that I describe below.\nWe are grateful to the reviewer for these positive comments.\nMajor Comments:\n[1] Overall, I think the authors make a reasonable argument for using a simple heat conduction\nmodel to estimate bottom-up and top-down SIF, given the many challenges in modeling slush and\nlateral water transport in firn. However, I have a couple of questions/thoughts on this approach.\n• First, it seems that top-down refreezing may still be important in the summer, since the\nthermistor measurements in Figure 2b show a period of almost 15 days at the end of August\nwhere temperatures at and above the ice slab drop below 0C and the ice slab itself begins\nto cool again slightly. Currently, the modeling approach only considers bottom-up\nrefreezing during the summer. I understand that there is more uncertainty in top-down\nrefreezing in the summer due to uncertainties about the overlying snow thickness and its\nchange with time.\nHowever, I would encourage the authors to think about whether there is a way that they can\nreasonably incorporate summer top-down refreezing (or quantify some uncertainty associated\nwith not including it) for their pan-Greenland estimates. I advocate this because the study\nencompasses some cold years where it is highly likely that the snowpack did not remain at 0C\nfor the entire melt season.\nWe agree that this issue is worthy of further consideration. We have undertaken a new sensitivity\nanalysis using a combination of PMR melt extent duration and our top-down SIF model running\nwith a fixed -10c surface temperature and 5% irreducible water content in the overlying snowpack.\nThis analysis indicates that while thin (0.5 m) snowpacks overlying the slush may result in some\nappreciable summer top-down SIF, potentially as much as 6-7 Gt over the 2017-2022 study period,\nemploying a thicker (1.0 m snowpack) reduces this value dramatically to less than 1 Gt. In relative\nterms, these findings indicate summer top-down SIF generates a maximum of ~10% of total SIF but\nlikely much less.\nGiven that top-down summer SIF is not a major contributor and in the interests of keeping the\nmain text a manageable length, we make a small note in the main text (L240-243). We have added\nthe full details of the analysis to the Supplementary Discussion (‘Top-down refreezing during\nsummer melt hiatuses’) along with a new Supplementary Figure 6.\nWe also note that here we considered for the first time the impact of irreducible water content in the\noverlying snowpack pore space, which does not contribute to SIF but must be refrozen before the\nsurface-originating cold wave can begin to form SIF in the slush matrix. This process also applies at\nthe start of autumn, so we modified our top-down autumn SIF runs to start with 5% irreducible\nwater content (in-between the classical 7% value of Colbeck (1974) versus the depth-dependent\nscheme of models such as MAR which use a range of 7 to 2% with increasing depth), rather than\n0% as previously. This reformulation has the impact of very slightly reducing top-down SIF\n(because it cannot commence until the overlying snowpack has frozen); this effect is visible in Fig.\n4c. We have updated the autumn model runs in both our field area and ice-slabs-wide, which\nresults in reduced autumn top-down SIF and therefore a slight reduction in slabs-wide total SIF.\n• Second, I think it would be more accurate to describe the results of this model as the “SIF\ncapacity” or “SIF potential” since this is a purely thermal model that assumes infinite and\ninstantaneous melt availability. Therefore, the mass of refrozen ice calculated by the model is\nreally an upper bound – the maximum that could be refrozen over the time period of interest\ngiven the initial thermal state at depth and the evolving air temperature. Melt\navailability would modify the true amount of SI that forms. I think this is a relatively minor\nreframing that can help avoid any overstating of the results.\nConcerning the model basis as thermal, assuming infinite/instant melt supply: We have added this\nsentiment using the suggested wording where we first introduce the model (L119-20).\n‘Melt availability would modify the true amount of SI that forms.’. We highlight that our slabs-wide\nresults already take this into account. We exclude any SIF which forms in excess of total seasonal melt\n(albeit according to MAR melt which as we highlight in the Discussion has its own limits). But in this\nsense our calculations are more realistic rather than simply maximum bounds, especially when also\nconsidering our use of van de Crux et al. (2024) T10m uncertainties to bracket our calculations.\nMore broadly, we are of the view that there remain substantial uncertainties and unconstrained\nparameters, especially in treatment of top-down SIF. These prevent us from being able to draw the\nconclusion that our central ice-slabs-wide SIF estimate is truly an indication of maximum SIF potential\n– it may be an over-estimate but it may alternatively be an under-estimate.\nMinor Comments:\nLine 81 – might note here that the cSMB is calculated from repeat firn cores so that the reader does\nnot have to go digging in the methods to figure out if you are using a climate model, SEB model,\nmeasurements, etc.\nAdded.\nFigure 3 – I am unsure what the question marks on the lateral boundaries of the snow, firn, and slush\nare intended to convey.\nAs noted in the caption, these question marks indicate unconstrained snowpack and slush matrix\nthicknesses. We have made some small adjustments to the figure (all question marks same colour,\nand added horizontal ticks between the question marks) which we hope make this clearer.\nFigure 4\n• What is the y axis on the histogram? I did not know how to interpret [-].\nThis was unitless on a 0-1 scale. We have multiplied to percentages and labelled accordingly.\n• Panel c, consider adding a legend or labels for the different green lines and the yellow line to\nthe figure itself. In the caption, consider mentioning the source of the temperature timeseries.\nWe have added a legend; the caption now indicates source of temperatures.\n• Either label or explain in the caption what A↑, A↓, and S↑ stand for. I figured it out\neventually, but it would be nice to have it spelled out somewhere.\nSorry – this information should already have been in the caption, now added.\n• Panel f – I was confused about how the standard deviation of the observed SIF plays into the\ndifferent modeling scenarios, since the model is (as far as I understand) calculating\nsummer SIF independent of the observations.\nWe agree that this was a confusing analysis. We have therefore modified our approach: we now\nestimate uncertainty in bottom-up SIF (both summer and autumn) by modifying the 10 m depth\ninitialisation temperature by +/- 1 deg C, much as we already did in our previous submission for the ice-\nslabs-wide analysis. Please see the revised Fig. 5 and the section ‘Total SIF at the visible runoff\nlimit’ (L170-209).\n• In general, I found the discussion of the total SIF calculations in the 10 kms above FS2\nconfusing. It was not totally clear to me how the summer SIF was calculated, what the\nrelationship to the measured SIF uncertainty was, orwhy that region was chosen for these\ncalculations.\nWe have reworked the text in this section to take account of major changes to the previous Fig. 4,\nmainly requested by ref. 3. We hope that these details are now clearer (L170-209).\nLine 214-216 – I would guess that autumn SIF seems to dominate because these ice slab detections\nabout the firn aquifers are likely capturing proto-aquifer regions that didn’t quite make the radiometric\ncut to be masked out a full-blown firn aquifers. So they would likely have a similar tendency to freeze up\nresidual water slowly over the autumn.\nYes, we agree; we have extended our explanation accordingly (L270-271).\nLine 471 – could be good to check the sensitivity of the results to this assumption, since it will impart\nsome degree of spatial bias.\nWe have now undertaken a sensitivity analysis, the results are insensitive to this assumption over\nthe slab thickness range (5, 15 m). See also our response to R1 about the same topic. (L538-541).\nSupplementary Figure 2 – could the x axis be stretched a bit? At first glance, this is not super convincing\nbecause it looks like there are three stacked dots at ~100 m distance with very different\nSIFs, rather than a sharp decline in SIF with distance, but that may be simply an effect of the plot’s\naspect ratio.\nAt the suggestion of R3, we have moved SFig. 2 into panel d of Fig. 1. We have also taken this\nopportunity to add slope information to the plot and to revamp the Results text (L75-108). See also\nour response to R3.6.\nSupplementary Figure 3 – please include a 1:1 line since the x and y axis scales are not the same.\nWe have made this change and incorporated this plot into the revised Fig. 4.\nCode availability: The README file for the code seems to be limited to some of the remote sensing\nanalysis and does not include any notes on the heat flux modeling code. However, the commenting in\nthe scripts themselves seems to be pretty good and the code flowchart looks very helpful. I would just\nencourage the authors to make sure there is an entry for every script in the README with at least a\nquick blurb about what purpose that script serves (what it calculates or what figure it reproduces, etc)\nand ideally notes on any data dependencies to run the script.\nThanks for taking the time to look at these extra materials. We have added extra information to the\ncode README file, particularly concentrating on the main analysis scripts. Given the complexity of\nthe underlying computations (use of GEE and a HPC in addition to local processing) and the multiple\ndata sources, the flowchart remains the most appropriate overview of how each script is linked to data,\nother scripts and outputs. We have therefore also improved the flowchart by sub-dividing it into\nanalysis sections and highlighting the contributions of the SIF model in yellow.\nReviewer #3 (Remarks to the Author):\nOverall\nI have spent another day reviewing this paper which for me seems to be too much at this late stage.\nThe authors have chosen to drop key aspects of their original study, including the SAR based retrieval\nfor superimposed ice formation (SIF) and also the timing work. They have instead added more\nmodelling, and they make use of new external datasets on ice slab extent and firn temperature. These\nare all good decisions, and the paper is clearer and more plausible for the changes.\nWe thank the reviewer for reviewing our revised manuscript and we are pleased that they found many of\nour major changes welcome.\nI don’t think the authors have addressed the main concerns raised by the reviewers. In several cases\nthey are simply rejected with arguments in place of evidence and I don’t think this is acceptable.\nI would have expected the rebuttal letter to summarise the main concerns and how they have been\naddressed. Instead, the authors summarise what they believe to be the most important issues instead\nof the major issues the reviewers identified.\nIn our summary section we aimed to signpost the most major changes we had made to satisfy the main\nconcerns. We highlight that R2, who gave a critical appraisal last time, now says that “The authors have\nput in a great deal of work to revise this manuscript and have resolved the majority of my concerns”.\nI still have major concerns with the study:\n1. Failure to adequately address reviewer concerns\n1.1. Reviewer 2 notes that the small study area may not be representative of the ice sheet; this is not\nproperly addressed.\nR2 themself raised no further concerns about representativeness in this round of reviews.\nThis said, we believe that R3 is referring to paragraph 2 of the previous R2: “I am particularly concerned\nthat the approach to estimating pan-Greenland SIF relies on a highly uncertain relationship between SIF\nand nCDD in a single ~30 km x 30 km region which cannot be validated anywhere else on the ice sheet,\nas well a background SIF rate that is assumed to not vary at all across the ice sheet.”. However, the\nstatistical cNDD approach was already completely removed in the last submission, while the\nbackground SIF varies as a function of initial temperature at 10 m depth.\n1.2. Reviewer 2 notes that new methods should be validated; this is still not the case with the wetness\nderived from Sentinel-1 data.\nR2’s request for “new methods” relating to Sentinel-1 to be validated were as follows:\n• Comment 3 called for sensitivity testing of the SAR wetness backscatter threshold (not\nvalidation, strictly speaking). We undertook this work as suggested. R2 had no further\ncomments to make on this part of the re-submission.\n• Comment 4 concerned the generalisation of our defunct cNDD model in conjunction with SAR-\nderived wetness to estimate slabs-wide SIF. However, as this cNDD model is no longer part of\nour analysis, then any such validation would be without merit.\n1.3. Reviewer 3 notes that runoff is not unknown and yet the authors continue to insist it is.\nThis statement is factually incorrect. We do not continue to insist that runoff is unknown. As we stated\nin our previous response, we already changed ‘unknown’ to ‘unmeasured’ and/or ‘unvalidated’ in the\nmanuscript text accordingly.\nA search on the previous version of the manuscript for the word ‘unknown’ yields the following uses:\n• L391: in the context of a slush layer of unknown thickness.\n• Supplementary L88: “However, it is unknown how much refreezing takes place inside the\nsnow/firn once surface melt starts but before the snow/firn interval is fully water-saturated.”.\nWe are content to delete the phrase ‘Runoff from ice slab areas remains unmeasured.’ from the\nsecond paragraph of our Introduction as it is not essential to the message of the paragraph (L31).\nAs we requested in our last response, we would be happy to be directed to any literature which has\nmeasured runoff at high elevations/the runoff limit, because we still cannot find any such studies. So\nfar, this reviewer did not cite any literature in their reviews.\nIf we are misunderstanding the reviewer’s concern, then we would be happy to receive a few sentences\nof explanation so that we stand a realistic chance of resolving it.\n1.4. Reviewer 3 notes that the authors experiment does not measure SIF directly, it measures the depth\nto the ice layer surface. The authors do not address this.\nWe think that this comment refers to R3.2, particularly “The main approach is to measure the depth to\nthe slab surface two years apart at a variety of locations. This is not an actual measurement of\nsuperimposed ice formation, it’s a crude approximation.” We had been under the impression that the\nspecific concern here related to their following sentence about slab mobility and its impact on revisiting\nthe same firn facies in successive years, but it seems we mis-interpreted, sorry.\nThe statement in the present comment is partially correct: our measurement of SIF is based on the\nthickness of ice which freezes onto the ice slab each year. This is similar to, but not the same thing as,\n“the depth to the ice layer surface”. As we explain in our Methods and Supplementary Methods, we\nconvert these ice slab thickness increases to metres water equivalent, which requires consideration of\nfirn density. Thus, in this sense, SIF is indeed not measured ‘directly’, but is calculated through density\nchange in the firn column over a given interval.\nAll this said, we do not use the terminology ‘direct measurement’ anywhere in the manuscript, so we\nare still at a loss as to what the concern is here.\n1.5. Reviewer 3 notes that the correlation between melting and proximity of surface water is not clear.\nThe authors do not address this.\n(Clarification: Considering the previous review comment 3.1, we believe that the reviewer means the\ncorrelation between SIF (not melting as stated above) and proximity of surface water.)\nWe concur that a more detailed look at the relationship between SIF and surface water proximity\nis warranted and thank the reviewer for raising this point again. Please see our response to\ncomment 5.\n2. There is no error estimate, meaning the findings are not confident.\nIn response to reviewer 1 the authors state that “the model does not come with uncertainty bounds\nbecause of the complex and nonlinear nature of the modelling.” This is not acceptable; an error\nestimate should be mandatory for publication.\nThis comment concerns the K-Transect SEB model of van As et al. (2017) whose outputs of melt and\nsurface temperature we use here. In lieu of uncertainty bounds, its performance in modelling melt is\nevaluated and validated with respect to observations, both in the original study and by us here.\nvan As et al. (2017) found that at the end of a model run spanning seven melt seasons, the model\nunder-estimated cumulative SMB by 0.5 m w.e at KAN_U. van As et al. (2017) used the textbook surface\nroughness parameters of 0.02 and 1 mm for snow and ice, respectively.\nHere, we evaluated the model’s melt in our field area for the 2022 melt season by reference to our in-\nsitu surface lowering measurements; the model produced 0.25 m w.e. melt available for runoff, in\nexcellent agreement with our observation-derived values of 0.23-0.27 m w.e. over the same time period\n(see Methods → ‘Melt available for runoff or SIF in field area’, L406-411).\nThe new data in fact lead to a revision of the retained SIF from 32 gigtaons to 60 gigatons between 2017\nand 2022. No uncertainty estimate is given either in the original or revised manuscript, so it is not\npossible to say whether these two estimates are in agreement within experimental error. But this is a\nkey take home message – the authors main conclusion comes without any assessment of its\nconfidence. It is also very obvious to this reviewer, but not the eventual readers, that their estimate changed by more\nthan a factor 2 during the review process without any assessment of its reliability. This is for\nme a major concern, and the authors should I believe include an estimated error on this primary conclusion.\nRespectfully, we reiterate that our central estimate changed by roughly a factor 2 because using ice\nslab extents from Culberg et al. (2024; based on SAR) increased slab extents by factor 2 compared to\nthe previous product (Jullien et al., 2023; based on individual OIB flightlines). This was in response to\nthe suggestion by R2 to use the Culberg product.\n3. The main finding, that SIF is only 2-3% of the runoff, is not stated.\nThe main finding – that 60 Gt of ice is retained – amounts to a 2 to 3% adjustment to the ice sheet runoff\nduring the same period, well within the estimated 16% certainty of that measurement. I think it would\nbe helpful for the authors to reflect this in their abstract.\nWhile we have no problem with the reviewer forming their own personal views on the main takeaways of\nthis study, for us as the authors the main finding is not that “SIF is only 2-3% of runoff”, hence this\ninformation is found only in the middle of the Discussion. With respect, this is a rather reductive view\nand indicates why we are reluctant to include this value in the abstract or concluding remarks where it\nwould lack the detailed context of the Discussion.\nIce-sheet-wide, SIF is indeed currently a relatively small runoff reducer. But our main finding and chief\nnovelty lies in using our field measurements to establish bounds on the likely magnitude of this process\nfor the first time. Locally, SIF can be enormous as a percentage, even in areas where visible satellite\nimagery is rather more suggestive of pervasive runoff. The importance of this process may also increase\nin future as slabs expand. We shine light on a previously neglected process which will help to reduce\nuncertainties such as the 16% range in the Slater et al. (2021) runoff estimate.\n4. The paper is too verbose, poorly structured, and still has plenty of loose language.\nThe abstract contains more loose language. For example, it is stated that “However, the amount of\nrunoff versus retention from these areas has not been measured.” but the runoff has certainly been\nmeasured, so why is this noted?\nWe see two possibilities for measurements.\n(1) Concerning in-situ measurements, we are not aware of such literature, so we are rather\nconfident in our view that our in-situ, stake-based measurements in the vicinity of the\nGreenland Ice Sheet’s runoff limits are truly novel.\n(2) Concerning the possibility of remotely-sensed measurements, studies exist based on a\ncombination of altimetry and firn compaction modelling (e.g. Slater et al. 2021 cited in\nmanuscript). But it is precisely the processes of SIF and ice slab expansion/thickening\nwhich are not yet adequately integrated into these firn compaction models (see also\nhttps://doi.org/10.5194/egusphere-2024-2750) yet which are fundamental to resolving the\ncause of the altimetry signal. We therefore don’t consider these studies to be equivalent to\nour measurement approach here.\nAs we stated in our previous response to reviews, we cannot find any literature which has measured\nrunoff vs retention explicitly in ice slab areas, but we would be happy to receive suggestions.\nIf “these areas” means the ice slabs then say so, otherwise I read it to mean “the lower accumulation\narea” as the previous sentence states. The phrases “at first order” and “most likely” are in there; what\ndo these mean ? Why say “the majority y of superimposed ice forms during the melt season” instead of\nstating the estimated value? This is really quite elementary stuff and should be picked up before\nsubmission by authors.\nWe have made changes as suggested. These changes include adding two statistics to the\nabstract, both of which are found in the section ‘Localised SIF continues into autumn’.\nThe paper structure doesn’t distinguish adequately between the observational and modelling aspects\nof the work, which are quite distinct and should be clearly separated to the reader. Instead there are\nsub-headings related to processes and conclusions which is not helpful at all. I suggest laying the\nresults out in a more tradiational manner; observations then models etc. There is also a tendency to\ndraw firm conclusions from the model (also in the sub-headings) when I am not so sure that this is truly\ndiagnostic as the models are primed to deliver results related to the proximity of water, for example.\nOverall I think too much of the paper is given over to the modelling when this is the least reliable aspect\nof the work, more emphasis could be given to the observations.\nAt this late stage and considering that neither R1 nor R2 raised similar concerns (indeed R2 is now very\npositive: “I also appreciated the attention to the written narrative, which is now much clearer than in the\noriginal submission”), we don’t see the basis for making such dramatic changes. Further, our choice of\nsub-headings reflects the NComms preference for ‘The Results and Methods section should be divided\nby topical subheadings’ (Guide for Authors).\n5. The in situ SIF data are not as strong as the authors would like us to believe.\nAlthough 3 sites close to rivers show high (0.5-0.6 mwe) SIF, one site also very close shows low (0.2\nmwe) SIF. All are within 100 m of a feature and so there is really not enough information to draw a firm\nconclusion other than that the regional rate of SIF is about 0.2 to 0.3 mwe and that rates can be 2-3\ntimes higher close to surface rivers. This is not the same as saying that “proximity to visible surface\nhydrological pathways strongly determines SIF magnitude”.\nWe agree with the spirit of this comment and have therefore removed the aforementioned phrase.\nWe have revised this section of the manuscript to consider not only distance to surface rivers but\nalso the role of local slope angle, supported by the new Fig. 1d (L75-108).\nI actually think the lateral water flow is the more interesting conclusion that could be drawn and\nexplored from this dataset; one could set up a simple model of a river surrounded by firn and assess the\npercolation.\nThis is a major undertaking which would demand a separate study. Indeed, we have recently published\na study dedicated to a first effort at modelling lateral meltwater flow through a slush field and assessing\nwhere the water table first rises to the surface:\nClerx N, Machguth H, Tedstone A, van As D. Modelling lateral meltwater flow and superimposed\nice formation atop Greenland’s near-surface ice slabs. Journal of Glaciology. 2024;70:e54.\ndoi:10.1017/jog.2024.69\nThe challenges here are also acknowledged by R2, who stated “Overall, I think the authors make a\nreasonable argument for using a simple heat conduction model to estimate bottom-up and top-down\nSIF, given the many challenges in modeling slush and lateral water transport in firn.”\nThe authors state that “at least the majority of SIF can be explained by direct refreezing of local surface\nmelt” but the actual fraction can be estimated from the data; why not do this?\nWe have now quantified this accordingly: “These results are evidence that away from surface\nhydrological features, 94 % of SIF (at 0.16 m w.e. yr-1; 68–100 % considering the standard deviation\nof our observations of ±0.05 m w.e. yr-1) can be explained by direct refreezing…” (L133-135).\n6. The figures don’t showcase the main results.\nFig.1 should be redrawn. The location map (panel a) is the largest component while the field\nmeasurements of SIF (panels b, c, and d) are the smallest. This should be the other way around; the\nfield measurements are the main data collected for this study and should be prominent. Supp Fig 2\nshould be in Fig 1 also and the chart annotated with the site names shown in Fig. 1 panels c,d,e; it’s\nprobably the most valuable new information provided.\nWe are grateful for this suggestion. We have shrunk panel a, enlarged panels b and c (with inset)\nwith site labels and 2 metre contours. The previous Supp. Fig. 2 now appears revised in panel d\nincluding the surface slope angle. We now comment on both distance to drainage and the slope\nangle in the text, particularly concerning SIF measurements 3u3 and FS1 (2022), which otherwise\ngive the appearance of outliers if considering distance alone. (L75-108).\nSupp Fig 3 shows a comparison between the model and observed SIF. This figure should also be in the\nmain text (again with field sites labelled) as it is the main evidence to support use of the model for\nregional and ice sheet wide estimates. Its not great though, is it (r=0.79)? You could use the rms\ndeparture as a very rough estimate of the model certainty, though the in situ measurements themselves\nare subject to uncertainty too. I guess its about 0.1 to 0.2 mwe/yr certainty and that’s comparable to the\nbackground SIF of 0.16 mwe which would mean a ~100% error.\nWe agree with moving this Figure into the main text; it can now be found in Fig. 5c. We added\nlabels for all sites where we could do this legibly. See also our response to the subsequent\ncomment about Fig. 4.\nRegarding the comment ‘It’s not great though, is it (r=0.79)?’. The sample size is small (n=10) so the\nstatistical test’s power is indeed never going to be great. Against the line of best fit most sites show\ngood correspondence; there is only one clear outlier, site 231, which we now comment on in more\ndetail in the text. We also note that SIF stakes are point estimates made using 9 cm diameter cores. In\ncontrast, our SIF calculations incorporate the use of SAR data with 25 m ground resolution, so it is\nunrealistic to expect a perfect model fit given the scale mismatch. We conclude that r=0.79 (R2 = 0.66,\np < 0.01; L177) is a pretty good outcome. Concerning ‘errors’, we refer the reviewer to our\ncomprehensive treatment of uncertainties throughout the manuscript.\nFig 4 shows maps of the wetness duration (b) and the % of melt that refreezes (e). These are interesting\nbut a map of the SIF with the field observations from Fig 1bc,d would be far more valuable thing to show\nas it would show clearly how the model compares to reality. The other panels in Fig. 4 could be moved\nto the supplementary materials.\nOn the basis of these suggestions we have split the previous Fig. 4 into two figures (Figs. 4 and 5).\nFig. 5a contains an overview of SIF calculated by our model (using a 0.6 m autumn snowpack) and\npanel b shows a detailed excerpt over our SIF stake locations. Panel c contains the measured-\nmodelled comparison previously found in Suppl. Fig. 3. Panel d remains similar to its previous\nversion, with modifications according to comments from reviewer 2.\nWe have reworked the main text accordingly to focus mainly on absolute SIF (L170-209),\nwhich\nremoves the need to supply extra maps of % melt refreezing in the supplementary materials.\nFig. 6. The components of SIF (panels b,c,d) can be moved to the supplementary materials; this would\ndeclutter the image and refocus it to just to the total SIF and its spatial / temporal variability.\nRespectfully, we disagree. These panels are fundamental to displaying our sensitivity analysis\nassociated with each SIF component, which R3 asked about in their main comment 2.\nResponse to reviews\nConcurrent superimposed ice formation and meltwater runo3 on\nGreenland’s ice slabs, version 3\nComments from Reviewer 2:\nThe authors have resolved my comments from the last round of review. I have a couple minor\npoints of clarification, but I do not think this requires another round with the reviewers.\nFigure 4, panel c – is the thin black line temperature? Please add to the legend or explain in the\ncaption. Maybe I just lost track in the text, but is the idea of the +/- 1C diNerence in 10 m\ntemperature for the bottom-up SIF uncertainty introduced anywhere in the text before this\ncaption?\nYes, the thin black line is temperature – we now note this in the caption.\nThe reviewer is correct, we forgot to introduce the idea of the +/- 1 C diNerence in the text. We\nnow do this in the final paragraph of the section ‘Localised SIF continues into autumn’.\nLine 146 – it would be helpful to add a sentence here on how you \"combined ... model-informed\nestimate of summer and autumn SIF to derive total SIF\". I assume it's as simple as calculating\nthe residual wetness time period for each pixel and taking the corresponding SIF components\nfrom Figure 4c, but that's not made totally clear in the text.\nIndeed, it is this simple. We have modified the cited sentence to now read “We summed our\npixel-by-pixel estimates of summer and autumn SIF introduced above…”"}, "author_response": ""}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors have resolved my comments from the last round of review. I have a couple minor points of clarification, but I do not think this requires another round with the reviewers.\nFigure 4, panel c – is the thin black line temperature? Please add to the legend or explain in the caption. Maybe I just lost track in the text, but is the idea of the +/- 1C difference in 10 m temperature for the bottom-up SIF uncertainty introduced anywhere in the text before this caption?\nLine 146 – it would be helpful to add a sentence here on how you \"combined ... model-informed estimate of summer and autumn SIF to derive total SIF\". I assume it's as simple as calculating the residual wetness time period for each pixel and taking the corresponding SIF components from Figure 4c, but that's not made totally clear in the text."}, "author_response": "Response to reviews\nConcurrent superimposed ice formation and meltwater runo3 on\nGreenland’s ice slabs, version 3\nComments from Reviewer 2:\nThe authors have resolved my comments from the last round of review. I have a couple minor\npoints of clarification, but I do not think this requires another round with the reviewers.\nFigure 4, panel c – is the thin black line temperature? Please add to the legend or explain in the\ncaption. Maybe I just lost track in the text, but is the idea of the +/- 1C diNerence in 10 m\ntemperature for the bottom-up SIF uncertainty introduced anywhere in the text before this\ncaption?\nYes, the thin black line is temperature – we now note this in the caption.\nThe reviewer is correct, we forgot to introduce the idea of the +/- 1 C diNerence in the text. We\nnow do this in the final paragraph of the section ‘Localised SIF continues into autumn’.\nLine 146 – it would be helpful to add a sentence here on how you \"combined ... model-informed\nestimate of summer and autumn SIF to derive total SIF\". I assume it's as simple as calculating\nthe residual wetness time period for each pixel and taking the corresponding SIF components\nfrom Figure 4c, but that's not made totally clear in the text.\nIndeed, it is this simple. We have modified the cited sentence to now read “We summed our\npixel-by-pixel estimates of summer and autumn SIF introduced above…”"}]

Decision letter