Greenland ice sheet runoff reduced by meltwater refreezing in bare ice

Domain: earth_science
Published: 2025-09-12
Editorial outcome: major_revision

Abstract
The contribution of Greenland Ice Sheet meltwater runoff to global sea-level rise is accelerating due to increased melting of its bare-ice ablation zone. There is growing evidence, however, that climate models overestimate runoff from this critical area of the ice sheet. Climate models traditionally assume that all bare-ice runoff enters the ocean, unlike porous firn, in which some meltwater is retained and/or refrozen. We used field measurements and numerical modeling to reveal that extensive retention and refreezing also occurs in bare glacier ice. We found that, from 2009 to 2018, meltwater refreezing in bare, porous glacier ice reduced runoff by an estimated 11–17 Gt a−1 in southwest Greenland alone, equivalent to 9–15% of this sector’s annual meltwater runoff simulated by climate models. This mass retention explains evidence from prior studies of runoff overestimation on bare ice by current generation climate models and may represent an overlooked buffer on projected runoff increases. Inclusion of bare-ice retention and refreezing processes in climate models therefore has immediate potential to improve forecasts of ice sheet runoff and its contribution to sea-level rise.

Introduction
Greenland Ice Sheet (GrIS) mass loss raised global sea level approximately 10.8 ± 0.9 mm from 1992 to 2018, and climate models forecast an additional increase of 70–130 mm by year 21001. Most of this predicted mass loss derives from increased meltwater runoff from the ice sheet’s ablation zone2–6. Within this critically important zone, winter snowpack melts entirely each summer, exposing dark, bare glacier ice which absorbs up to three times more sunlight than bright snow7. Warming air temperatures and reduced summer snowfall have exposed larger areas of bare ice in recent decades2,3,8, driving enhanced surface melting in the ablation zone2–4. Understanding the fate of meltwater from Greenland’s growing bare-ice zone is therefore critical for accurate modeling of sea levels9.

Climate models track energy and mass budgets on the ice sheet surface and are the primary tools available to predict future ice sheet runoff10,11. However, there is substantial observational evidence that these models overestimate runoff from bare ice surfaces12–17. In the ablation zone of Greenland’s melt-intensive southwest sector, proglacial and supraglacial river discharge measurements reveal up to 67% less annual meltwater release to surrounding oceans than climate model calculations12,14–16. Supraglacial lakes which form on the ice sheet surface fill at slower rates than predicted by climate models17, and direct measurements of supraglacial runoff are overestimated by 21–58% during peak summer melt conditions13. Ice sheet mass changes are overestimated by 21–47% relative to GRACE satellite gravity retrievals18,19, and satellite laser altimetry measurements indicate that surface melt rates are overestimated by 14–40%20. Similar discrepancies arise from comparisons with point ablation stake measurements, ranging from an underestimation of 17% at one ablation-zone site to overestimations of 10–43% at four other sites21,22. These discrepancies do not appear to be explained by errors in modeled surface energy balance components, which largely match in situ meteorological observations13,21,22.

A more plausible explanation for overestimation of ice sheet runoff by climate models may be their conceptual treatment of bare ice surfaces. Climate models traditionally treat bare ice as an impervious, high-density substrate with no capacity to retain water11,23,24. Accordingly, runoff produced on bare ice is instantly credited in its entirety to sea level23, despite growing field reports of non-trivial meltwater retention on or within bare ice12,13,25–29. These studies suggest that bare ice surfaces are not necessarily impervious, but can behave rather like snow and firn, with some portion of generated meltwater retained either through refreezing or liquid storage in pore spaces30,31. Because the majority (>78%)24 of Greenland’s meltwater runoff is generated from bare ice, even small amounts of meltwater retention or refreezing would constitute a substantial (and overlooked) component of the ice sheet surface mass balance that may explain consistent reports of climate model runoff overestimation12,13,17.

Here, we formally test this hypothesis of retention and refreezing by pairing field measurements of meltwater runoff with a numerical surface mass balance model explicitly designed to simulate spectral radiation and thermodynamic heat transfer in bare glacier ice (Supplementary Materials). We compare meltwater runoff simulated by our model with outputs from two regional climate models, a global climate reanalysis, and in situ measurements from a well-studied surface catchment on the southwest GrIS ablation zone—including supraglacial river discharge13, surface ablation, and physical bare ice properties26,32. We then demonstrate that discrepancies between modeled and observed ice sheet runoff can be attributed to retention and refreezing processes in bare ice of the southwest GrIS ablation zone, where the majority of GrIS surface mass loss originates2,3.

Results
Field evidence of meltwater refreezing in bare glacier ice
In July 2016, we conducted a field experiment at the “RB” catchment in southwest Greenland (Fig. 1, see also Fig. S1). This catchment covered about 60–63 km2 of ice sheet surface and an elevation range of about 1200 to 1360 meters above sea level (m a.s.l.), depending on the year13,17,33. Meltwater runoff generated within this catchment drained through a supraglacial river network leading to a terminal moulin—a vertical shaft where surface water enters the subglacial system (Fig. 1a). During a previous field study13 in July 2015, we measured supraglacial river discharge at a cross-section upstream from the terminal moulin. Note that discharge refers to the volumetric flow rate, measured as the volume of water passing through a channel cross-section per unit of time, whereas runoff represents the volume of water produced over the upstream contributing area per unit of time. In this paper, ‘observed runoff’ denotes measured discharge (in m3/s) from upstream contributing areas, which is then converted to cumulative volume (m3, and ultimately Gt) for comparison with climate model runoff (Methods).

Our field measurements of RB catchment runoff in July 2015 revealed a substantial overestimation in predicted runoff by regional and global climate models, ranging from +21–58%13. To investigate this discrepancy, we revisited the site from 6 to 13 July 2016, collecting 168 consecutive hourly measurements of supraglacial river discharge34 concurrent with three-hourly measurements of ice surface lowering35 from a network of ablation stakes (Fig. 1b–d) we installed near the gauging site, which comprise the key field datasets of this study (Methods). Shallow ice cores obtained at the same site during July 2016 revealed a porous and meltwater-saturated weathering crust on the bare-ice surface, with an average density of just 690 kg m−3 within the top meter of ice26. Critically, we observed refrozen meltwater within this weathering crust (Figure. S2), leading to our hypothesis that the refreezing of meltwater in bare ice could explain observed discrepancies between climate model predictions and runoff observations13,26.

Climate models overestimate meltwater runoff from the bare-ice ablation zone
Consistent with previous studies12–17, we find evidence that climate models overpredict observed meltwater runoff (Fig. 2). By the end of the 6–13 July 2016 field experiment, climate model runoff ranges from 7% lower (MERRA-2) to 58% higher (RACMO2.3p3) than observations (Fig. 2), similar to +21–58% we previously reported for 201513. Among the climate models examined here, RACMO2.3p3 most closely reproduces observed albedo (Fig. S4) and net radiation (Figs. S5–S6), yet severely overpredicts cumulative runoff despite accurate representation of these critical surface energy balance components.

To explain these discrepancies we found between climate model predictions and runoff observations, we developed ‘SkinModel’, a process-based numerical model of ice sheet meltwater runoff (Methods). SkinModel represents the ice sheet surface as an infinitely thin, hydrologically impermeable ‘skin’ layer of high-density (900 kg m−3) ice, emulating assumptions used in the current generation of climate models (Methods Eq. 1)23. Note that SkinModel is not a statistical emulator, but rather a simplified numerical model designed to reproduce the surface energy balance method of calculating ice sheet meltwater runoff used in climate models.

SkinModel forcings for RB catchment were provided by meteorological observations from the PROMICE/GAP KAN_M automatic weather station36 located ~2 km from the RB catchment’s northeastern boundary (67.067oN, −48.835oW; 1270 m.a.s.l.) (Fig. S1). To establish SkinModel’s ability to effectively emulate climate model meltwater runoff predictions, we deliberately force it with surface albedo values from each climate model (rather than KAN_M albedo values), and find its predictions are virtually identical to climate model predictions (Fig. 2; dashed versus solid lines). Critically, all model forcings besides albedo are kept consistent across these emulator simulations, signifying that ice surface albedo alone can explain virtually all differences between climate model runoff predictions shown in Fig. 2.

Next, we force SkinModel with values of albedo recorded at KAN_M and find that modeled runoff is 42% higher than measured runoff (Fig. 2; green solid line). Weather station forcings thus yield runoff predictions much like RACMO2.3p3, the model with the most accurate albedo for this time and location (Fig. S4). Similar results are found using these same methods to reinterpret our earlier field measurements from 201513 (Figs. S6–S8). Even larger runoff overestimation is simulated if MODIS satellite albedo values37 (Fig. S4) are used as model forcing. This reveals that among the climate models examined here, accurate representation of ice surface albedo lead to overpredicted runoff, whereas accurate predictions of runoff result from overestimated albedo.

Runoff overestimation explained by meltwater refreezing in bare ice
We developed a second process-based numerical model of ice sheet meltwater runoff that we call IceModel38 (Supplementary Materials). IceModel updates an earlier model of spectral radiative and thermodynamic heat transfer in glacier ice39 with a field constraint32 on shortwave radiation absorption by dark impurities present within Greenland’s bare, ablating ice. In contrast to SkinModel, IceModel simulates an ice column with time-varying ice, air, water vapor, and liquid water contents (Methods Eq. 4), informed by observations of physical ice surface properties made during a multi-year series of field campaigns on the GrIS bare-ice ablation zone12,13,16,26,32,33,40. A critical feature of IceModel is that it allows sunlight to penetrate bare ice, thus allowing melt within the ice subsurface41 rather than expending all available energy upon an infinitely thin ‘skin’ surface layer.

IceModel simulations reveal that, during daylight hours, penetration of shortwave radiation produces an approximately isothermal ice column nearly 1 m thick that stores latent heat in the form of liquid meltwater (Fig. S9), a phenomenon independently validated by our field observations of saturated porous ice extending at least one meter below the ice surface26. Lateral transport of meltwater generated within isothermal bare ice is constrained by its low horizontal hydraulic conductivity (~10−2–10−4 m h−1)42,43 and impermeable cold-ice lower boundary. At night, the cold polar boundary layer cools the ice surface to −5°C on average during July at this location, as indicated by IceModel simulated surface temperatures which closely track observations (Fig. 3). These subfreezing temperatures drive refreezing of subsurface liquid meltwater at rates approaching 1 mm h−1 between 02:00–04:00 localtime, when the ice surface temperature drops as low as −6°C. The runoff overestimation by SkinModel (Fig. 2) can therefore be explained by nocturnal refreezing of liquid meltwater stored within the upper decimeters of the porous bare ice matrix. Inclusion of this process enables IceModel to closely reproduce observed cumulative runoff, well within catchment boundary uncertainty (Fig. 2).

Climate model runoff over-estimations in the RB catchment (Fig. 2) are mirrored in ice surface lowering rates measured at our ablation stake network and the KAN_M automatic weather station (Fig. 4a). Unlike SkinModel and the climate model surface energy balance it emulates, IceModel predicts internal mass loss below the ice surface (Fig. 4b) caused by subsurface absorption of transmitted solar radiation. Below ~0.4 m, IceModel over-predicts ice density by up to ~150 kg m−3 relative to core measurements, suggesting deeper light penetration and/or lateral meltwater delivery may increase porosity at depth more than the model allows. This bias yields a column-averaged simulated density of 770 kg m−3—13% higher than the observed mean of 681 kg m−3 yet still well below the canonical 900 kg m−3 assumed in most climate models (vertical green line, Fig. 4b). Importantly, these simulations and field measurements reveal that substantial melt occurs beneath the ice surface, where meltwater is retained within porous bare ice, undetected by climate model simulations of the ‘skin’ energy balance, and available for nocturnal refreezing.

Methods
Field datasets of ice sheet surface discharge, ablation rate, and ice density
Hourly supraglacial river discharges from RB catchment (Fig. S1) were measured during two field campaigns (20–23 July 2015 and 6–13 July 2016)13,34 with a SonTek River Surveyor M9 Acoustic Doppler Current Profiler (ADCP). The M9 was mounted on a SonTek HydroBoard II and escorted across the RB main-stem channel with a bank-operated Tyrolean system (Fig. 1). Channel geometry and flow velocity measurements were transmitted in real time to a bank-operated computer running the River Surveyor software. During each measurement hour, between 3 and 9 sub-hourly measurements were recorded. These raw data were converted to channel flow rate [m3 s−1] in post-processing following data quality-control workflows optimized for the supraglacial environment13,34. The standard deviations of these sub-hourly measurements were used to estimate ±0.6 m3 s-1 measurement uncertainty. Error bars in Fig. 2 depict two standard deviations of cumulative discharge obtained by accumulating hourly discharge uncertainty.

Measurements of ice surface elevation change were provided by a network of twelve bamboo ablation stakes we installed near our field camp (67.0496oN, 49.0201oW, 1215 m a.s.l.) (Fig. 1). Stake locations were selected by generating random distance-direction pairs from a common center until a representative set of sites including bright white ice and dark impurity-laden ice were selected within an area covering ~0.1 km2 (Fig. 1). Stakes were drilled 3 m deep into the ice and allowed to freeze-in for 24 hours prior to initiation of ablation stake measurements, which were then recorded at nominal 3-hour intervals continuously from 12:00 localtime on 6 July 2016 to 23:00 localtime on 12 July 2016. Freeze-in was confirmed prior to each measurement to infer potential vertical displacement of these stakes due to melt-out at their base; none was observed. Prior to each measurement, a 24 × 24 cm square wooden ablation board was placed at the base of the stake and oriented to true north (see Fig. 1). This board operated as a datum from which the stake height above the ice surface was measured. Cumulative changes in stake height were converted directly to ice surface elevation change and summarized in terms of their mean and standard deviation across stakes for comparison with simulated melt rates (see Fig. 4). Field datasets of ice density from shallow ice cores, ice porosity, and ice liquid water saturation within excavated boreholes used to supplement this analysis are described in reference26. These direct ablation stake measurements were supplemented by hourly changes in ice surface elevation from a pressure transducer ice ablation assembly53 installed on the PROMICE/GAP KAN_M automatic weather station36 (67.0667o N, 48.8327o W). During July–September 2016, KAN_M station was located ~2 km from the RB catchment at an average elevation of ~1270 m a.s.l36. (see Fig. S1). These pressure transducer measurements have an estimated uncertainty of ±2.5 cm liquid water equivalent53, depicted as vertical error bars in Fig. 4a. To compare simulated melt (m liquid water equivalent) with these pressure transducer and ablation stake ice surface elevation changes, melt was converted to ice thickness change using the conversion: Δh = M(ρl/ρ), where Δh is ice thickness change, M is melt, ρl = 1000 kg m−3 is the density of liquid water at atmospheric pressure and the triple point temperature of water (273.16 K), and ρ is glacier ice density. To address known variations in near-surface ρ in the study area (see Fig. 4b), uncertainty ranges on these conversions were estimated using lower and upper assumed values of ρ (600–900 kg m−3), depicted as shaded bounds in Fig. 4a. For IceModel simulations, a central estimate was obtained directly from the modeled vertical column ice density (ρ = 770 kg m−3), depicted as a dotted line in Fig. 4a.

Catchment boundaries and surface classification from satellite and airborne datasets
Catchment topography and surface classifications for RB catchment were obtained from WorldView-1 and WorldView-2 satellite images of the ice sheet surface and aerial imagery collected with an RGB sensor mounted on an uncrewed aerial vehicle54 (Fig. S1b). These aerial images were used to reconstruct the ice sheet surface topography and to perform surface classification of snow, water, and bare ice as described in reference13. Surface topography was reconstructed using Agisoft PhotoScan Pro stereo-photogrammetry software. Surface classification was performed with a k-Nearest Neighbors algorithm yielding 3.1% snow cover during the 6–13 July 2016 field experiment, and 6.5% for the 20–23 July 2015 experiment13 (Fig. S1b).

WorldView-1 and WorldView-2 satellite imagery and associated high resolution stereo-photogrammetric digital elevation models were used to delineate the RB contributing catchment area following methods in reference13. These methods combine traditional automated surface delineation from digital elevation with visual delineation of surface stream networks, flow direction, and channel heads visible in WorldView imagery following reference55. Interior channel heads (initiation points of channels that drain into the catchment) yield a minimum estimate of catchment area. Areas of internal drainage to moulins and crevasses within the catchment boundary were removed from this lower catchment area estimate (Figs. S1 and S12). Outer channel heads (initiation points of channels that drain away from the catchment) yield an upper, maximum plausible estimate of catchment area. The optimal “best guess” catchment area was delineated by tracing the digital elevation-based boundary and adjacent inner and outer channel heads in the high-resolution WorldView imagery and adding back areas that flow into the catchment and subtracting areas that flow out of the catchment. This approach yields contributing catchment area confirmed by the direction of actively flowing water tracks, an improvement upon digital elevation-based methods56.

To complement our detailed experimental study of RB catchment, we synthesized several publicly available datasets of supraglacial and proglacial discharge measurements, for comparison with modeled meltwater runoff across the southwest sector of the ice sheet (Table S1 and Table S2). Here, supraglacial refers to discharge measured on the ice sheet surface, and proglacial to discharge measured in terrestrial rivers emerging from the ice sheet periphery which drain supraglacial runoff routed through englacial and subglacial pathways. Proglacial catchments include LG (Leverett Glacier)44 and AK4 (Akuliarusiarsuup Kuua River’s northern tributary)46. Supraglacial catchments include two SLV (Supraglacial Lake Volume) catchments (SLV1 and SLV2)17, and the 660 catchment57, a supraglacial catchment near the western ice sheet margin nested within AK4. Catchment boundaries and discharge gauge installations for these sites are depicted in Fig. S1.

Catchment boundaries for SLV1 and SLV2 were obtained using identical methods to those described above for RB catchment, provided by reference17. The LG and AK4 catchments boundaries were provided by the PROMICE Greenland Liquid Water Discharge dataset58. These delineations were produced with the GRASS GIS software single-flow direction from eight neighbors (SFD-8) drainage basin delineation algorithm, with inputs from ArcticDEM v7 100m gridded ice sheet surface topography and BedMachine v359 150m gridded ice thickness. The upper, middle, and lower contributing area boundaries for the LG and AK4 catchments correspond to upper, middle, and lower prescribed values of the flotation factor k (1.1, 1.0, and 0.9, respectively) used in the GRASS basin delineation algorithm. These values account for the evolving subglacial water pressure and hydraulic routing throughout the summer melt season58.

Climate model data
Climate model outputs were provided by the Modèle Atmosphérique Régional version 3.11 (MAR3.11)60, the polar (p) version of the Regional Atmospheric Climate Model version 2.3p3 (RACMO2.3p3)47 and version 2.3p2 (RACMO2.3p2)22, and the global climate model Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2)61 (Table S2). MERRA-2 data are publicly available on a global 0.5o by 0.625o latitude-longitude grid (~56 × 28 km for the grid cell intersecting RB catchment), with land ice albedo and runoff provided on a 3-hour timestep, and land surface forcings provided on a 1-hour timestep. MAR3.11 data were provided at ~15 km horizontal resolution and 1 h timestep forced by European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis62. RACMO2.3p3 data were provided at ~11 km horizontal resolution and 3 h timestep forced with ECMWF ERA-Interim reanalysis63. Note that RACMO2.3p3 is the latest version of the regional climate model RACMO2 and data for the period 2009–2012 were not acquired for this study. Runoff outputs from RACMO2.3p2 were obtained from the PROMICE Greenland Liquid Water Discharge dataset58, exclusively for comparison with LG catchment discharge during the period 2009–2012.

To address the different grid resolutions and formats of the climate models, we adopted a common grid with 5 km horizontal resolution projected onto the National Snow and Ice Data Center WGS84/NSIDC Sea Ice Polar Stereographic North (SIPSN) coordinate system (EPSG:3413). This common grid was provided with the PROMICE Greenland land-ice albedo product37 (discussed further in Methods Section “Models simulations for individual catchments and southwest Greenland”), which was the highest-resolution gridded product used in this study. Catchment area-weighted values of climate model output were computed by intersecting the climate model horizontal grids with a bounding box surrounding each catchment boundary. The coordinates of bounded grid cell centroids and the associated gridded climate model outputs were resampled onto the common 5 km grid using Delaunay triangulation and natural neighbor interpolation. The 5 km re-gridded values were then converted to catchment area-weighted values using a conservative remapping of the overlapping area of each grid cell and the catchment polygon (Fig. S12 and Section S6). Catchment areas were computed in the National Snow and Ice Data Center Equal-Area Scalable Earth version 2.064 projection (Table S1).

Comparing climate model runoff with measured discharge
The climate models evaluated in this study do not represent the processes of lateral flow, channel routing, or vertical infiltration of meltwater within bare ice. Consequently, these processes were not represented by either IceModel or SkinModel. Although the MAR model includes a routing delay designed to simulate the passage of meltwater through the ice sheet to its periphery13, this mechanism is not applicable at supraglacial catchment scales, nor was it applied to the MAR runoff evaluated here. To compare modeled runoff with measured discharge—which inherently includes real-world catchment routing delays65—we focused our analysis on cumulative values [m3] of runoff and discharge, obtained by summing hourly fluxes [m h−1] multiplied by upstream contributing area [m2]. This approach is necessary due to the intrinsic discrepancy between instantaneous runoff and discharge rates and aligns with our primary objective of assessing the impact of meltwater refreezing in bare ice on total runoff from the ablation zone.

Our methodology for comparing cumulative runoff to cumulative discharge is reflected in our primary evaluation metric—the percent difference (over/underestimation) between cumulative values of modeled runoff and measured discharge (Table S3). To account for diurnal variations in runoff and routing delays, we calculated these percent differences not at a single final timestep, but as averages over the last 24 h of each observation period, which vary in length from approximately three days (RB catchment in 2015) to seven years (AK4 catchment) (Table S2). This averaging mitigates the chance alignment (or misalignment) of modeled runoff and measured discharge when assessing their agreement on the final timestep, while retaining a focus on the total meltwater runoff over the observation period.

Operational climate models including MAR3.11, MERRA-2, and RACMO2.3p2 currently do not include mechanisms for subsurface (internal) radiative heating in bare glacier ice and its consequent effects on meltwater runoff and refreezing. This simplification of bare ice persists despite emerging interest from the climate modeling community, as evidenced by a recent early inclusion of refreezing in RACMO2.3p347, which has yet to undergo empirical validation. Although we emphasize our independent implementation of these processes in IceModel and its empirical validation, for this study, we obtained two configurations of the RACMO2.3p3 model output: a standard version incorporating internal energy absorption in bare ice, provided for the period 2012–2018, and a control version without these processes referred to as RACMO2.3p3-WIE47 (Without Internal Energy), provided exclusively for the period 2012–2015. Runoff outputs from these two configurations are compared with field observations in Section S5 and Table S3 (see also Table S2).

SkinModel and IceModel description
SkinModel and IceModel both solve the zero-dimensional surface energy balance66:

χQsi(1 − α) + Qli − εσT4sfc + Qh + Qe + Qc = Qm (1)

where χ [unitless] allocates the incoming shortwave solar radiation flux, Qsi [W m−2], into a ‘skin’ surface component and a subsurface component, α [unitless] is ice surface albedo, Qli [W m−2] is incoming longwave radiation flux, εσT4sfc [W m−2] is longwave radiation flux emitted by the ice surface, ε = 0.98 [unitless] is ice surface emissivity, σ [W m−2 K−4] is the Stefan-Boltzmann constant, Tsfc [K] is ice surface temperature, Qh [W m−2] is sensible heat flux, Qe [W m−2] is latent heat flux, Qc [W m−2] is conductive heat flux at the surface, and Qm [W m−2] is energy available for meltwater production. For SkinModel simulations, we set χ = 1 following a traditional surface energy balance approach39, meaning all shortwave solar radiation is absorbed at the surface. For IceModel simulations, χ was computed directly from the fraction of solar radiation absorbed in the top layer of the one-dimensional ice column, updated at each timestep with a two-stream radiative transfer model described in the Supplementary Text (Section S3). Monin-Obukhov similarity theory was used to obtain Qh and Qe as described in ref.39.

SkinModel calculates Qc by solving the one-dimensional heat transfer equation:

ρi ci (∂Ti / ∂t) = ∂ / ∂z (ki ∂Ti / ∂z) (2)

where ρi = 900 [kg m−3] is a reference value for glacier ice density, ci = 2093 [J kg−1 K−1] is specific heat capacity of ice at 273.16 K and constant pressure, Ti [K] is ice temperature, t [s] is time, z [m] is the vertical coordinate (positive down), and the glacier ice thermal conductivity, ki [W m−1 K−1], was modeled as a function of ice temperature67:

ki = 9.828 exp(−0.0057 Ti) (3)

IceModel calculates Qc and subsurface meltwater production/refreezing by solving the one-dimensional thermal and spectral radiative heat transfer equation39,68:

∂H / ∂t = ∂ / ∂z (ke ∂T / ∂z) − ∂q / ∂z (4)

where H [J m−3] is enthalpy per unit volume of the ice, liquid water, and water vapor mixture, ke [W m−1 K−1] is the effective thermal conductivity of the ice, liquid water, and water vapor mixture, and q [W m−2] is the subsurface net solar radiative flux at depth z. The model updates an earlier version39 with field-measured values for the spectral absorption coefficient of glacier ice69 and a conservative enthalpy formulation70,71. Detailed descriptions of IceModel and its numerical implementation are available in the Supplementary Text (Section S3).

Equation 2 and Eq. 4 were resolved on a uniform mesh with 4 cm node spacing to a depth of 20 m, using a fully implicit time integration and a 15 min timestep. The upper boundary condition, Tsfc, was estimated by solving Eq. 1 by Newton-Raphson iteration39. At the lower boundary, a zero heat flux condition was assumed. Hourly values for near-surface air temperature, wind speed, relative humidity, air pressure, α, and Qsi and Qli to solve Eq. 1 were provided by the KAN_M automatic weather station for RB simulations and by MAR3.11 for regional simulations with daily MODIS (Moderate Resolution Imaging Spectroradiometer) albedo α from the PROMICE Greenland land-ice albedo product37. These hourly forcings and daily albedo values were linearly interpolated to the 15 min model timestep.

The net solar radiative flux divergence, ∂q/∂z [W m−3], was evaluated with a two-stream radiative transfer model with 118 spectral bands39,72 (Supplementary Text Section S3). A key input to this model is the spectral solar radiation extinction coefficient, kλ, which controls the vertical distribution of absorbed solar radiation within subsurface ice layers. To account for the melt-enhancing effect of dark impurities present within glacier ice, values for kλ were obtained from direct measurements of spectral flux extinction within glacier ice in Greenland’s western ablation zone32. Absorption by impurities was further constrained by KAN_M automatic weather station and MODIS albedo values prescribed at the upper boundary of the two-stream model.

Models simulations for individual catchments and southwest Greenland
IceModel and SkinModel were solved on a 15 min timestep for the experimental periods 6–13 July 2016 and 20–23 July 2015 for the RB catchment, with model inputs from the KAN_M automatic weather station36, supplemented by catchment area-weighted climate model output (e.g., climate model albedo for SkinModel emulator simulations). These simulations informed our initial catchment-scale experimental analysis depicted in Figs. 2–4.

Regional-scale simulations for Greenland’s southwest sector depicted in Fig. 5 were solved on a 15 min timestep for the decade 2009–2018 by running IceModel and SkinModel at 5 km horizontal grid spacing, with each grid cell representing an independent simulation of one-dimensional vertical heat transfer (i.e., without horizontal heat transfer). The model domain was Ice Sheet Mass Balance Intercomparison Experiment (IMBIE)1 GrIS southwest sector, which covers ~12% of the ice sheet surface area but accounts for ~34% of modeled surface runoff24. Model forcings for these simulations were provided by MAR3.11 on an hourly timestep (linearly interpolated to 15 min) for the decade 2009–2018 (RACMO data were provided on a three-hour timestep for 2012–2018; relative humidity, wind speed, and surface pressure were not provided). Quality-controlled daily ice sheet surface albedo from MODIS was provided by the PROMICE Greenland land-ice albedo product37 on a 5 km horizontal grid, which was used as a common grid for both regional-scale simulations and for resampling climate model outputs (see also Methods Section “Climate model data”).

Two sets of IceModel and SkinModel simulations were performed for the southwest sector domain: one in which MAR3.11 data were used exclusively as model forcing, and a second in which MODIS albedo was used in place of MAR3.11 albedo (see Fig. 5 and Table S3). These simulations were designed to isolate the critical role of albedo in modulating runoff predictions, and to mirror our catchment-scale simulations in which climate model albedo values were used in place of weather station albedo values (see Fig. 2). Note that IceModel forced with observed surface albedo from MODIS serves as our benchmark simulation.

Southwest sector IceModel and SkinModel simulations were restricted to grid cells contained within the Greenland Ice Mapping Project (GIMP)73 ice mask with MODIS albedo between 0.29 and 0.55 on at least 75% of days during July and August from 2000 to 2018, following refs. 2, 50. This procedure was used to isolate bare-ice grid cells from firn-covered grid cells and terrestrial land. The 75% bare-ice frequency threshold is more conservative than the 50% threshold used in the GrIS SMB Intercomparison (GrSMBMIP)50, and yielded a 61,975 km2 bare-ice surface area for the southwest sector, used as a common mask for computing regional-scale bare-ice runoff and refreezing estimates (e.g., Figs. 5a and S11).

IceModel and SkinModel runoff estimates for RB, LG, AK4, 660, SLV1, and SLV2 catchments depicted in Fig. 5 and Section S5: Supplementary Runoff Comparison were computed from these gridded regional-scale simulations using the catchment area-weighted conservative remapping described in Methods Section “Climate model data”. Regional-scale runoff and refreezing estimates reported in the Main text (e.g., Fig. 5) and in Fig. S11 were estimated for two time periods to illustrate a range of plausible values of bare ice refreezing and associated runoff reduction: a “lower” (15 June to 15 August) conservative proxy for seasonal bare-ice exposure duration, and an “upper” (1 June to 31 August) traditional one2. IceModel runoff and refreezing estimates for these two time periods were obtained directly from these regional-scale gridded simulations, restricted to grid cells within the common bare-ice mask.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_2": "Overall comments:\nThe paper proposes a very interesting and potentially significant message that the refreezing of meltwater in the near\nsurface weathered crustal bare ice is a major reason for the consistent overestimation of runoff from bare ice regions of the\nGreenland Ice Sheet by regional climate models.\nThe paper goes further by calculating that refreezing in bare ice could reduce runoff estimates from RCMs by between 11\nand 17% and make their outputs more consistent with observations of runoff from both internal ice catchments and from\nproglacial rivers.\nThe idea is conceptually appealing as anyone who has walked across a freezing bare glacier surface late at night or early\nmorning will attest to. The concept has been previously proposed (e.g. in Smith et al 2017). The challenge for this paper is to\nconvincingly model and quantify the impact, then assess its wider significance.\nThe paper is very well written. It does rely on a large amount of Supplementary material which really does need to be read in\ndetail alongside the main paper. This makes it a much longer read than is apparent at first but this is often the case for short\nformat journals.\nFor the most part the work undertaken is very rigorous and the conclusions drawn are well supported by the lines of\nevidence presented. However, it is fiendishly difficult to keep track of which model outcomes support the overall conclusions\nand which outcomes are more ambiguous.\nOverall, I would say that there are some outcomes at the larger scales that undermine the neat logical argument that holds\nup at the smaller scale. These need to be tackled.\nI do think the case for the superiority of IceModel as a more realistic process model has to be backed up with some\nsensitivity testing and independent calibration.\nTo help guide a way through this review I have proposed a few key assumptions that I think need to be upheld for the\nquantitative results of the paper to stand.\n1) That the measurements of runoff are the best indicator available of what runoff from the different catchments actually is.\n2) That the SkinModel is a good emulator of different RCMs when the different RCM albedo schemes are used with\nSkinModel (evidence from Fig 1).\n3) If the only substantive difference between IceModel and SkinModel is that the former includes subsurface melt and\nrefreezing and the latter does not (as explained in Methods), then if IceModel uses the same albedo scheme as SkinModel\nand performs better in comparison against runoff measurements, then its improved performance is due to the inclusion of\nbare ice refreezing.\n4) If SkinModel emulates RCMs then we can use the improvement of IceModel over SkinModel (i.e. (3) above) as a reliable\nindicator of the improvement to RCMs if they were to include bare ice refreezing.\n5) If we accept that IceModel is the most accurate and the most physically realistic model then the quantitative amounts of\nrefreezing that it predicts are a reliable estimate for the actual effect that bare ice refreezing has on surface mass balance.\nAddressing point 1 above, the authors refer to a substantial body of previous work (e.g. refs 13 and 34) in measuring surface\nrunoff from the internal catchments within the bare ice ablation area. Although the possibility that some runoff has been\nmissed from these internal catchments is not ruled out in this earlier work, the team have done a very thorough job and make\na convincing case that their measurements are as good a data set as we are likely to get to test modelled runoff against.\nMeasuring the huge discharges that emerge from the ice sheet into proglacial streams is a very different challenge that is\noutlined in ref 43. Here the most significant uncertainty is whether the runoff measured in the proglacial stream actually\nrelates to the catchments being modelled. This is very hard to quantify but a decent attempt to account for this uncertainty is\nincluded through varying the modelled catchment boundaries according to different assumptions about subglacial water\nrouting given the surface and bed topography of the upstream ice sheet. So point (1) is fine.\nAddressing point (2) above, SkinModel emulates the MAR and RACMO models very well at the small scale over a short time\nperiod (as shown in Fig1). Different albedo schemes, (RACMO, MAR, MODIS) are applied to the “RCM emulator” SkinModel\nto show that differences in albedo alone cannot account for the discrepancy between modelled and measured runoff. So for\nthe RB catchment from July of 2015 and 2016, RACMO albedo is closest to the MODIS derived values (and the AWS\nmeasurements) whilst the MAR and MERRA-2 RCM albedos are much higher (Fig S9). Even though RACMO has more\naccurate albedo and likely does a good job of estimating surface melt (as shown in Figs S5 and S6), its continued persistent\noverestimation of runoff must be due to some other process omission/error. The argument is that MAR and MERRA-2 only\nhave a better match with observations because they use unrealistic albedo schemes. SkinModel forced with AWS albedo\nmeasurements does a slightly better job than the “RACMO emulator” but it is still not a good agreement to measured runoff.\nIceModel, driven by the same AWS albedo but also including bare ice refreezing, improves the results significantly and\ntherefore the inclusion of this process is quite logically proposed as being the key point of difference.\nSo SkinModel emulates the MAR and RACMO models very well at the small scale over a short time period, but there is a\nsignificant discrepancy between SkinModel (MAR albedo) and MAR at the larger LG catchment scale over period 2009-\n2012 (Fig 4c). So it could be that SkinModel may not actually be a particularly good emulator of MAR (or RACMO?) at larger\nspatial and temporal scales. Why does this discrepancy arise? Albedo at different temporal and spatial scales may be\nrelevant here.\nAt the small scales shown in Fig S9, MAR and MERRA-2 have the highest albedos and RACMO has lower albedo.\nHowever, Fig S11 shows that for areas of the SW sector above c.900m (which is most of the ablation area), RACMO albedo\non average is actually higher than MAR albedo. For the LG catchment from 2009-2012, RACMO overestimates runoff by\n30% compared to a MAR overestimate of only 8%. This is despite RACMO quite probably having higher average albedo for\nthis overall catchment and time period (based on the long-term, sector-wide data shown in FigS11). So the differences\nbetween RACMO and MAR are significant and are not all about different albedo. This starts to mess up the assumptions that\nheld at the smaller scales above, and indicates that processes other than albedo and refreezing may explain differences. We\nsee results at the larger scales that don’t fit neatly to the patterns seen at small scale. In summary, point (2) holds at smaller\nscale but is more ambiguous at larger scales.\nPoint (3) above is challenged by the fact that for the larger catchment scale shown in Fig 4c, SkinModel (MODIS albedo)\ndoes just about as good a job as IceModel (MODIS albedo) over these larger scales. The first tracks the upper error band on\nobservations of runoff, the second tracks the lower error band and all three overlap if error bands on model outputs are\nconsidered.\nIf SkinModel does not emulate RCMs consistently then point (4) above is more uncertain and it is possible that the apparent\nimprovements that IceModel shows over SkinModel may be due to other factors that are not explored. The improved\nperformance of IceModel is evidenced by the data shown in Table S3 and is represented by the average percentage\ndifference between modelled and measured runoff, µ. I think there are other ways to assess the different model performance\nthat does not simply take averages of values that are calculated over very wide ranging spatial and temporal scales. As a\nfirst pass the absolute differences should be considered rather averaging positive and negative numbers. If you do that\nIceModel – AWS albedo comes out best. However, once you start applying some kind of weighting/significance criteria to\nthe averages (e.g. where model runs over larger areas and over longer time periods have greater weighting - see attached\nfile for an example), then the assessment of best model performance changes. Excluding runs that didn’t include the full\n2009-2012 LG catchment run, then the best model is SkinModel with MAR3.11albedo, followed by MAR 3.11. All this means\nthat IceModel doesn’t necessarily lead to improved model performance which undermines confidence in points (4) and (5)\nabove.\nAnother issue for point (5), i.e. using IceModel to try to quantify the bare ice refreezing process, is that there is no sensitivity\ntesting of IceModel shown or referred to. For the paper to stand, we have to be convinced that IceModel is the most accurate\nand that it is the most physically realistic. What parameters are IceModel (and for that matter, SkinModel) most sensitive to?\nThe ice core data is helpful in demonstrating the improved process accuracy of IceModel over other models. However, such\ndata is also small scale. The paper would be strengthened if we could see that IceModel was calibrated using larger scale\ndata sets that convinced us that it could be upscaled to the catchment and ice sheet sector scale. More specifically, are\nSkinModel and IceModel calibrated using data that is independent of the model validation data? There may not be space in\nthe paper to include all this, but there are 55 pages of Supplementary material but nothing included that answers these quite\nfundamental points.\nSorry for such a long review but I find myself torn. I like the concept of the paper and it’s great that the authors have obtained\ndata at the smaller scales that really backs up the overall idea. However, it seems to unravel at the larger temporal and\nspatial scales which makes the quantification of this effect much harder to have confidence in. It might involve too much work\nfor a revision for the authors to turn this around unless they can demonstrate how I’ve misunderstood and overlooked some\nkey points that alleviates these concerns. It will make a great paper if they can!\nSmaller points:\nFigure 2,3 and 4 captions start with a statement that reads a bit like a paper conclusion. These should be avoided and just\nexplain to reader in very factual terms what the Figures show.\nIf SkinModel (MAR albedo) referred to in Fig 4 is essentially the same model as MAR Emulator (SkinModel) referred to in\nFig1, then they should be named consistently to avoid confusion.\nWhy isn’t SkinModel (RACMO albedo) shown in Fig4? If it does a better job of emulating RACMO then this would alleviate\nthe issue that emulators don’t transfer well across scales."}, "author_response": "Dear Editors and Reviewers,\nWe thank the reviewers for their thorough and insightful comments, which have greatly helped\nus improve the manuscript. Below, we provide a summary of revisions made in response to\nreviewer comments, followed by a point-by-point response to each comment.\nIn the point-by-point response, reviewer comments are included in plain font, followed by our\nresponses in bold italics.\nMatt Cooper\nSummary of revisions made in response to reviewer\ncomments\n● Reorganization of Results and Discussion:\n○ Moved the regional-scale comparison with climate model runoff from the Results\nsection to an expanded Discussion section.\n○ Ensured that the Results section focuses on our core findings from field\nobservations and numerical modeling at RB catchment.\n○ Expanded discussion of uncertainties in catchment delineation and their potential\ninfluence on model comparisons at proglacial sites.\n○ Added new Supplementary Discussion (Section S4) to clarify uncertainties in\nproglacial discharge comparisons.\n● Clarifications on Model Assumptions and Performance:\n○ Clarified that IceModel was not calibrated but is a physics-based model with\nparameters constrained by literature values.\n○ Clarified the role of SkinModel as an emulator at small scales but not intended for\nlarge-scale quantitative emulation of RCM runoff.\n○ Addressed concerns about SkinModel’s consistency with MAR and RACMO at\ndifferent spatial scales, in particular the role of albedo in controlling runoff biases.\n● Supplementary Material Updates:\n○ Moved a key figure from the Supplementary Materials to the main article.\n○ Added a new discussion of model uncertainties to Supplementary Section S4.\n○ Revised Supplementary Figures to better illustrate the role of albedo in runoff\noverestimation.\n● Formatting Improvements:\n○ Revised figure captions for clarity and consistency.\n○ Standardized terminology (e.g., “SkinModel (MAR albedo)” instead of “MAR\nEmulator (SkinModel)”).\nReviewer #1 (Remarks to the Author):\nReview of: Greenland Ice Sheet runoff reduced by meltwater refreezing in bare ice\nSubmitted to: Nature Communications\nReviewer: Benjamin Hills\nSummary.\nCooper et al. observe that meltwater runoff measured in the Greenland ablation zone is\nlower than would be predicted by climate models. They hypothesize that the reduction in runoff\nis associated with meltwater retention and refreezing in a near-surface layer of\nweathered ice, much like the relatively well-studied water retention in firn (Harper et al.,\n2012). They test their hypothesis with a pair of 1-dimensional thermodynamic numerical\nmodels they developed. They find that the model which considers light penetration and\nsubsurface melting is a better representation for their measurements of meltwater runoff. The\narticle is well written and a great fit for Nature Communications after the minor revisions I\nsuggest below.\nGeneral Comments.\nIn my opinion, the difference in model output between SkinModel and IceModel is one of the\nmore important results in this body of work. If you agree, I believe that some reorganization of\nthe Results section could help emphasize that result to all readers. Currently, “Results” is broken\nout into subsections on: fieldwork, model overestimation of meltwater (SkinModel), explanation\nof overestimation (IceModel), upscaling to SW Greenland. The separation between explanations\nof SkinModel and IceModel could make it hard for the reader to appreciate what you are doing\nhere (even though I do like your clarifying statement on L137-139). I would consider\nreorganizing as follows:\n● Intro (pretty much as you have it)\n● Results (field results only)\n○ You could consider moving Figs 2 and 3 up here, although I understand why you\nwant to lead with what you now have as Fig 1.\n● Analysis\n○ Bring your explanation of SkinModel and IceModel\n○ Compare to climate model output at the local scale (currently Fig 1)\n○ Compare to climate model output at the regional scale (currently Fig 4)\n● Discussion/Conclusions (pretty much as you have it)\nPerhaps this organization does not work well for Nature Communications but that is how I would\nwrite it to emphasize your two models and the important difference between them (take it or\nleave it).\n● Thank you for the suggestion and for clearly summarizing our study. We agree\nthat emphasizing the differences between SkinModel and IceModel is important.\nWe implemented a modified version of your suggested re-formatting, restructuring\nthe Results and Discussion sections to better emphasize our core results at RB\ncatchment versus our regional analysis.\n● Specifically, we moved the retrospective runoff analysis and sector-scale runoff\nreduction to the Discussion section, which adheres to Nature Communication’s\nformatting guidelines which request a dedicated Discussion section that offers\n“an extended analysis of the results and their comparison to the literature rather\nthan a short summary or conclusion”.\nOur implemented version of your suggested re-formatting is summarized below. Please\nnote that we have included an additional component of the Discussion section\n(highlighted in bold) that was not included in your list:\n● Intro (remains largely unchanged)\n● Results (field results only)\n○ Compare to climate model output at the local scale (currently Fig 1)\n○ Explanation using SkinModel and IceModel (Fig 2 and 3)\n● Discussion/Conclusions\n○ Compare to climate model output at the regional scale (currently Fig 4)\n○ Contextualize in terms of sector-scale total runoff predictions (9-15%)\nWhile this does not exactly adhere to your suggested re-formatting, we believe it aligns\nwith the overall goal of better highlighting our core results derived from our detailed field\nstudies and associated numerical modeling (SkinModel and IceModel), now formally\nseparated from the extended sector-scale analysis in the Discussion section.\nMy second large take-away from your article is the difference between the local scale (and\nsupraglacial discharge) as opposed to regional scale (and proglacial discharge). In Figure 1 we\nsee that the measurements are at the lower end of the model spread, whereas in Figure 4 the\nmeasurements are closer to the middle of the spread (especially in c). There are several\nimportant points which come to mind for me which I do not think were mentioned in the article.\nFirst is the fundamental difference in the catchments. The physical processes important for your\nsurface discharge measurements at RB will be entirely captured by your measurements and\nmodeling efforts (as far as I can tell), but the proglacial discharge measurements have water\nthat has moved through englacial and subglacial environments which you do not consider.\nCould basal melting provide an additional component to discharge (you may be able to rule that\nout with a sentence or two)? What about groundwater? Could there be some groundwater flux\nnot accounted for in the discharge measurements?\n● We acknowledge the distinction between supraglacial and proglacial discharge\nand agree that additional clarification is needed.\n● To address this, we have added a new opening paragraph to the Discussion\nsection explicitly addressing potential contributions from basal melting,\nsubglacial routing, contributing area uncertainty, and groundwater flux.\nAdditionally, we have made several edits throughout the Discussion section to\nhighlight uncertainties in the proglacial runoff comparisons more clearly. Please\nsee the tracked-changes document.\n● To further address this, we have added a new discussion of uncertainty to the\nSupplementary Materials. Please see Section S4: Supplementary Discussion.\nI also wonder whether the uncertainty in defined catchment areas are sufficiently captured (I\nknow they came from a different study). Compared to your surface catchment at RB, the\nregional catchments are more difficult to constrain with confidence and if subglacial water gets\nrouted in unexpected ways it could elevate (or suppress) the measured proglacial discharge. I\ndon’t think this discredits your study, but worth mentioning some of these uncertainties.\n● Thank you for highlighting this important source of uncertainty. To address this,\nwe made three changes:\n○ We primarily incorporated this into the revisions mentioned in the previous\ncomment, where we emphasize the additional uncertainty in proglacial\ndischarge comparisons.\n○ We restructured the Results section by moving the comparison with\nproglacial discharge (LG and AK4 catchments) into the Discussion section.\nThis clarifies that our core results derive from our detailed field studies at\nRB catchment.\n○ We edited the first paragraph of the revised Discussion (where the\nproglacial analyses are presented) to better clarify that our intention is to\ndemonstrate the possible broader implications of our findings over larger\nspatial and temporal scales.\n● Please see these revisions in the tracked-changes document.\nI am generally confused by the use of the term “emulator”. That term is often used to describe\nthe statistical representation of a complex physical model which is highly computationally\nexpensive (e.g. Conti et al., 2009; Gu et al., 2018). Here, I believe you are running the forward\nmodels (IceModel and SkinModel) directly. By that I mean that you are actually calculating the\nnumerics on your prescribed 1-D mesh (which would not be done for an emulator). If I am\nmistaken, please correct me and add some context in the methods on the statistical emulation.\n● Thank you for identifying this potential source of confusion. To address this, we\nmade two changes:\n○ We added a statement explicitly stating that SkinModel is not a statistical\nemulator (see L126 of the tracked changes document).\n○ We changed the legend text of Fig. 1 (revised Fig. 2) to read “SkinModel\n(MAR albedo)” rather than “MAR Emulator (SkinModel)”. While we retained\nour use of emulator in the text, this change makes the legends of Fig. 1 and\nFig. 4 (revised Fig. 2 and Fig. 5) consistent.\n● To clarify, it is correct that SkinModel (and IceModel) are run as forward models,\nand there is no statistical emulation involved. We use the term emulator in its\ndictionary-definition sense—to reproduce the function or action of another\nsystem. Specifically, SkinModel emulates the surface energy balance method\nused by climate models to calculate meltwater runoff. We believe that emulate and\nemulator remain appropriate terms in this case, and we did not find suitable\nreplacements.\n● Additionally, in the second paragraph of the section “Climate models overestimate\nmeltwater runoff from the bare ice ablation zone”, we describe SkinModel as a\n“process based numerical model of ice sheet meltwater runoff ... emulating\nassumptions used in the current generation of climate models”. While emulator\nhas gained a specific meaning in statistical modeling, our usage is conceptually\nsimilar, except that we’ve designed a numerical rather than statistical emulator.\nSkinModel is much simpler to run than a climate model, but effectively emulates\none for the narrow purpose of calculating meltwater from glacier ice.\nThe supplementary material is referenced many times in the main article. I appreciate the\nthoroughness of the authors here, but so many references can disrupt the readability of the\narticle. My personal opinion is that only the most important references to the supplementary\nmaterial should be used, and that any figures which are critical to explain the main narrative\nshould be included in the article itself (see my comment under Figures below).\n● Thank you for highlighting the impact of supplementary material references on\nreadability.\n● To address this, we carefully reviewed all supplementary references and removed\nor consolidated those that were not essential for understanding the main text.\n● Additionally, we have moved key supplementary figures into the main article, as\nsuggested previously.\n● Please see the tracked-changes document for specific modifications.\nSpecific Comments.\nL21 – Currently not clear whether the statement starting with “From 2009-2018…” is based on\nyour work. I would start with something like “We found that, from 2009-2018, meltwater…”\n● Edited as suggested.\nL24 – Same as above, I believe this statement could be more clear if you are explicit that it is\nbased on prior work, something like “The mass retention explains [the evidence from prior\nstudies] that runoff could be overestimated…”\n● Edited as suggested.\nL43 – Some readers may not immediately understand that SW Greenland is the same “ablation\nzone” that you had been talking about above.\n● Edited as suggested.\nL56 – is “substrate” the correct word here?\n● We retained “substrate”, which we believe is correct because it refers to both a\nsurface and its underlying material, which accurately describes how the ice sheet\nis represented in models (all of which include some representation of subsurface\ntemperature to compute the conductive flux).\nL65 – I would remind the reader again with “…this hypothesis [of retention and refreezing]…”,\nmostly because folks who skim the article may jump to this paragraph.\n● Edited it as suggested.\nL76 – Is “RB” defined somewhere? This is the first time I see it.\n● In previous work by Smith et al., 2019, we referred to this catchment as “Rio\nBehar”, in honor of our late colleague Alberto Behar who participated in previous\nfield studies in this region of Greenland. However, in comments on an earlier draft\nof this paper, we were asked to refrain from assigning place-names which could\nbe offensive to Greenlanders. We replaced “Rio Behar” with “RB” to reduce\npossible offense. We’ve cited prior work each time RB is introduced so that\nreaders can surmise this catchment is the same one we previously reported on.\nL86 – Somewhat confused by the definition of runoff. Discharge/area would give units of m/s but\nyou discuss runoff in Gt/a in the abstract and plot cumulative runoff in m3 (e.g. in fig 1).\n● Thank you for highlighting the potential confusion between discharge and runoff.\nTo clarify, we have revised the manuscript text as follows:\n○ “Note that discharge refers to the volumetric flow rate, measured as the\nvolume of water passing through a channel cross section per unit of time,\nwhereas runoff represents the volume of water produced over upstream\ncontributing areas per unit of time. In this paper, ‘observed runoff’ denotes\nmeasured discharge (in m3/s) from upstream contributing areas, which is\nthen converted to cumulative volume (m3, and ultimately Gt) for\ncomparison with climate model runoff (Methods).”\n○ Please see L89 of the tracked changes document.\n● Please also see Methods section “Comparing climate model runoff with measured\ndischarge” for additional clarification on units and model-observation\ncomparison.\nL106 – “To investigate these discrepancies…”? explain makes it sound like you knew the\nanswer coming into the study.\n● The discrepancies mentioned here refer to those reported in the preceding\nparagraph, which are the first results of the study. To ensure clarity, we edited the\nsentence to read: “To explain the discrepancies we found …” rather than “To\nexplain these discrepancies …”.\nL131 – You could just cite your Zenodo doi for IceModel instead of the supplement.\n● As suggested, we added a citation for the Zenodo doi. We retained the\nparenthetical reference to the supplementary material because it contains a formal\ndescription of the numerical model.\nL137-139 – This is one of the most important sentences in my opinion, and I fear it may get lost\nin the current organization. See my general comment on organization of the results section\nabove.\n● Thank you for highlighting this opportunity to enhance clarity. While we were\nunable to fully implement the suggested reorganization, we believe the distinction\nbetween SkinModel and IceModel is clearer now that the Results section is\nfocused entirely on RB catchment.\n● To explain why we organized our results in this manner, we believe it is necessary\nto first explain the variation between RCM predictions (using SkinModel), then\nexplain the variation between predictions and observations (using IceModel).\n● We believe this is essential in light of the background given in the Introduction,\nwhich describes a consistent picture of runoff overestimation, whereas our July\n2016 field study revealed a more nuanced (and potentially more impactful) picture,\nwhere runoff overestimation could be masked by albedo biases.\n● By first explaining this, we establish consistency with the background material,\nand then go further by showing how refreezing could resolve the overestimation.\nL147 – Not immediately clear why “July”? I would say “summer months” or “during the time of\nthe field investigation” or similar.\n● We write “July” because our two field studies were conducted in July 2015 and\nJuly 2016. If we write “summer” then we need to specify which months and also\nthe Northern hemisphere. We decided to keep our current phrasing since we\nspecifically refer here to the average nighttime temperature during July from the\nKAN-M air temperature data, and then point to the model versus observation\ncomparison in Fig. 2 (revised Fig. 3), which is specific to the field study period in\nJuly.\nL167 – Are you truly applying your models to “all bare ice areas of SW Greenland”? or just to\nthe LG and AK4 catchments?\n● Yes, we applied our models to all bare ice areas of SW Greenland and then\nextracted runoff for LG and AK4 (and 660, SLV1, and SLV2) from these gridded\nsimulation outputs, using the conservative remapping method described in the\nMethods and in the Supplementary Materials (Fig. S33-49). We did not run\ndedicated simulations for LG or AK4 catchments. The only dedicated\ncatchment-scale simulations we ran were for RB catchment, to conduct the\ndetailed albedo-sensitivity emulator simulations. Please see L518 in the Methods\nsection where these simulations are described in detail.\nL244 – Why two standard deviations? One would be more typical and I think appropriate here.\nEither way, it is obviously going to be more certain than the models which is true even for two.\n● We used two standard deviations for consistency with the catchment-area runoff\nuncertainty depicted by the shaded bounds in Figs. 1 and 4 (revised Figs. 2 and 5),\nas we consider this conceptually comparable to the confidence level implied by\ntwo standard deviations. While this choice is subjective, our method for\nestimating catchment-area uncertainty is highly conservative, as described below\n(also see Methods L370).\n● In brief, we delineated the “middle” catchment using traditional DEM-based\nmethods (with custom adjustments based on the channel heads method\ndescribed below), while the “lower” and “upper” bounds were delineated by\nmanually tracing channel initiation points (heads) of actively flowing streams\nvisible in WorldView imagery. The upper estimate was created by connecting the\nheads of channels that were visually confirmed to drain away from the RB\ncatchment divide. This likely overestimates the true contributing area, as it\nincludes some surface area that drains to those channels, away from the RB\ncatchment.\n● Conversely, the lower estimate was created by connecting channel heads that\nwere visually confirmed to drain into the RB catchment. This likely underestimates\nthe true contributing area, as it excludes surface area above these channel heads\nthat in reality drains to them, into the RB catchment.\n● Thus the upper estimate prioritizes capturing all possible contributing areas at the\ncost of including some non-contributing areas (Type I error; false positives), while\nthe lower estimate ensures that all included areas drain into the catchment at the\ncost of excluding some contributing areas above the identified channel heads\n(Type II error; false negatives).\n● Overall, we subjectively consider this to be conceptually similar to the notion of\n(at least) two standard deviations, so we used two standard deviations on the\ndischarge for consistency.\nL306 – In the caption of supplemental figure S2 you give the total area of the stake survey as\n~0.1 km2. I believe that is a relevant number to include in the main document, could be in this\nparagraph which is why I mention it here.\n● Edited as suggested.\nL330 – Giving water density to 5 significant digits is implying a certainty that you do not have,\nbut I do not feel strongly about changing it (certainly changing to 1000 will not change your\nanalysis at all).\n● Edited as suggested.\nL380 – I am assuming that the “flotation factor” is an input to the GRASS algorithm? Would be\nuseful to state that.\n● Correct. Revised as suggested.\nL434 – Can you strengthen the statement by dropping “to the best of our knowledge”? Seems to\nme like this list of authors would know if it were otherwise.\n● Edited as suggested.\nL478 – I am assuming that there are many 1-d models run in parallel? All with different surface\nboundary conditions? Might be worth clearing that up unless it is already clear and I missed it.\n● Correct. We edited the sentence by adding “... by running IceModel and SkinModel\nat 5 km horizontal grid spacing, with each grid cell representing an independent\nsimulation of one-dimensional vertical heat transfer (i.e., without horizontal heat\ntransfer).”\nL502 – “5 km horizontal grid spacing”. Should be made clear (as in previous comment) that this\nis not a model mesh spacing (i.e. including horizontal energy transfer) but parallel 1d models. I\ndon’t think it matters as the horizontal temperature gradient would be negligible, but worth\nmaking clear what you did.\n● As requested, we have made it clear that horizontal energy transfer is not\nrepresented. Please see our reply to the previous comment.\nL546 – I had trouble with the Leverett Glacier doi:\nhttps://doi.org/10.5285/17c400f1-ed6d-4d5a-a51f-aad9ee61ce3d. But its possible that the BAS\nserver was having issues when I checked. Anyway, double check me.\n● The link works for us. Perhaps there is a line-break or extra space. If problems\npersist, please let us know and we will rectify it.\nL551 – it seems like the PROMICE data portal may have changed to:\nhttps://promice.org/download-data.\n● Edited accordingly. Thank you for catching this.\nFigures\nI like Figure S2. I think it could elevate figures in the main article. That is, consider including (a)\nin Fig 1 of the main article and (b-d) in Fig 3 of the main article. If not this, then you should do\nsomething to make clear to the reader that Fig 1 is local scale compared to Fig 4c and 4d.\n● As suggested, we moved Fig. S2 into the main article. It is now Fig. 1. Since\nNature Communications allows five display items, we made it a dedicated figure.\nThe study area map in Fig. 4 (revised Fig. 5) displays each catchment boundary in\na regional context which we believe conveys the scale differences appropriately.\nReferences\n● Conti, S., Gosling, J. P., Oakley, J. E., & O’Hagan, A. (2009). Gaussian process\nemulation of dynamic computer codes. Biometrika, 96(3), 663–676.\nhttps://doi.org/10.1093/biomet/asp028\n● Gu, M., Wang, X., & Berger, J. O. (2018). Robust Gaussian stochastic process\nemulation. Annals of Statistics, 46(6A), 3038–3066. https://doi.org/10.1214/17-AOS1648\n● Harper, J. T., Humphrey, N. F., Pfeffer, W. T., Brown, J., & Fettweis, X. (2012). Greenland\nice-sheet contribution to sea-level rise buffered by meltwater storage in firn. Nature, 491,\n240–243. https://doi.org/10.1038/nature11566\nReviewer #1 (Remarks on code availability):\nThe model code for IceModel and SkinModel is extremely well documented, including examples\nto run.\nWe thank the reviewer for their insightful feedback. These revisions enhance the clarity\nand impact of our manuscript.\nSincerely,\nMatt Cooper\nReviewer #2 (Remarks to the Author):\nOverall comments:\nThe paper proposes a very interesting and potentially significant message that the refreezing of\nmeltwater in the near surface weathered crustal bare ice is a major reason for the consistent\noverestimation of runoff from bare ice regions of the Greenland Ice Sheet by regional climate\nmodels.\nThe paper goes further by calculating that refreezing in bare ice could reduce runoff estimates\nfrom RCMs by between 11 and 17% and make their outputs more consistent with observations\nof runoff from both internal ice catchments and from proglacial rivers.\n● We appreciate this summary of our study and the recognition of its significance.\nTo clarify, our study finds that refreezing, as estimated with our own models,\naccounts for 11–17% of the total runoff from snow and firn predicted by the MAR\nregional climate model (with nearly identical values for RACMO). We demonstrate\nthat incorporating refreezing into our models reconciles discrepancies between\nmodel predictions and observations. However, the extent to which RCMs would\npredict lower runoff and/or better align with observations if they accounted for\nrefreezing in bare ice depends on multiple factors beyond the scope of this study.\nThe idea is conceptually appealing as anyone who has walked across a freezing bare glacier\nsurface late at night or early morning will attest to. The concept has been previously proposed\n(e.g. in Smith et al 2017). The challenge for this paper is to convincingly model and quantify the\nimpact, then assess its wider significance.\nThe paper is very well written. It does rely on a large amount of Supplementary material which\nreally does need to be read in detail alongside the main paper. This makes it a much longer\nread than is apparent at first but this is often the case for short format journals.\nFor the most part the work undertaken is very rigorous and the conclusions drawn are well\nsupported by the lines of evidence presented. However, it is fiendishly difficult to keep track of\nwhich model outcomes support the overall conclusions and which outcomes are more\nambiguous.\n● We appreciate the positive assessment of our work’s rigor. To improve clarity, we\nhave restructured the Results section to more clearly distinguish between model\noutcomes that directly support our conclusions and those that introduce\nadditional complexity. Specifically, we have moved the section “Independent\nvalidation across the southwest Greenland Ice Sheet” into an expanded\nDiscussion section, ensuring that the Results section remains focused on our RB\ncatchment field study. This restructuring clarifies that our core results derive from\nour own field experiments, while the retrospective analysis of previously\npublished data serves as an extended assessment of the broader significance of\nour findings.\n● This restructuring aligns with Nature Communications formatting guidelines,\nwhich request a dedicated Discussion section which offers “an extended analysis\nof the results and their comparison to the literature rather than a short summary\nor conclusion”.\n● Additionally, we detail the ambiguity of the proglacial comparison in a new\nSupplementary Discussion section (revised Section S4).\nOverall, I would say that there are some outcomes at the larger scales that undermine the neat\nlogical argument that holds up at the smaller scale. These need to be tackled.\n● We have addressed this concern through specific responses to each comment\nbelow. In general, we have expanded the Discussion section to explicitly\nacknowledge and contextualize the differences between small- and large-scale\nresults, with additional discussion of key uncertainties at both scales.\nI do think the case for the superiority of IceModel as a more realistic process model has to be\nbacked up with some sensitivity testing and independent calibration.\n● We would like to clarify that IceModel was not calibrated, and we are uncertain\nwhat is meant by “independent calibration.” IceModel is a physically based model\nwhose predictions depend on initial and boundary conditions, forcing data, and\nparameters constrained by established values from the published literature, as\ndetailed in the Methods section and Supplementary Information. The only\nparameter derived from our field measurements is the ice absorption coefficient\nspectrum, which differs fundamentally from model calibration, where parameters\nare randomly varied to maximize agreement with observations. No such\ncalibration was performed, nor would it be appropriate given the research\nquestions we are addressing.\n● To address this comment, we added an acknowledgment in the Discussion\nsection (see L354 of the tracked changes document) that further sensitivity testing\nis needed, both for our models and the climate models we examined. While\nprocess-level sensitivity studies similar to ours exist, to our knowledge, the\nclimate models we examined have not undergone parameter sensitivity testing or\n“independent calibration” (again, the meaning of this term is unclear). However, if\nsuch studies exist, we would be grateful to learn of them and properly\nacknowledge them.\nTo help guide a way through this review I have proposed a few key assumptions that I think\nneed to be upheld for the quantitative results of the paper to stand.\n1) That the measurements of runoff are the best indicator available of what runoff from the\ndifferent catchments actually is.\n● We agree that our supraglacial runoff measurements are among the best available.\nHowever, proglacial runoff measurements come with inherent uncertainties. To\naddress this, we have moved our comparison with proglacial measurements at\nAK4 and LG into the Discussion section. We also note that this comparison\noccupies only a portion of a single figure (Fig. 4c-d, now 5c-d in the revised text).\nOur conclusions are primarily based on our detailed field studies and numerical\nmodeling of a physical process that we previously identified but which, until now,\nhas received little direct numerical simulation or empirical verification in terms of\nits broader significance.\n2) That the SkinModel is a good emulator of different RCMs when the different RCM albedo\nschemes are used with SkinModel (evidence from Fig 1).\n● To clarify, our core conclusions do not rely on the emulator assumption at large\nscales. To address this, we have:\n○ Removed the statement from the caption of Fig. 4 (revised Fig. 5) that\ndescribed SkinModel as a proxy for climate model predictions.\n○ Added a new opening sentence to the Discussion section which clearly\nstates that our sector-scale simulations and retrospective runoff analysis\nserve to explore the broader significance of our core finding that refreezing\nreduces runoff from bare ice.\n○ Added a discussion of limitations to the emulator approach to the\nSupplementary Materials, please see Section S4.3.\n● While SkinModel emulates RCMs well at small scales when using their respective\nalbedo schemes, it is not expected to exactly replicate RCM runoff predictions at\nall scales, times, or locations. As described in the paragraph beginning on L110,\nSkinModel is a process‑based numerical model of ice sheet meltwater runoff that\nemulates assumptions used in current-generation climate models. Specifically, we\ndescribe it as a “climate model surface energy balance emulator”. At the RB\ncatchment scale, this emulation allows us to explain the variation among RCMs in\nterms of albedo for that time and location.\n● The ability of SkinModel to emulate different RCM predictions at the RB catchment\nscale is critical for interpreting our experimental results—otherwise, one might\nmistakenly conclude that MAR and MERRA-2 are accurate despite their clear\nalbedo biases. However, this emulator function is not necessary to establish the\nbroader significance of refreezing over larger spatial and temporal scales. This\ndistinction is why we use the term “emulator” in the legend of Fig. 1 (revised Fig.\n2) and in the text when describing the RB catchment-scale simulations but not in\nFig. 4 (revised Fig. 5), where SkinModel and IceModel are used to quantify the\nrunoff reduction due to refreezing rather than to interpret RCM spread.\n3) If the only substantive difference between IceModel and SkinModel is that the former includes\nsubsurface melt and refreezing and the latter does not (as explained in Methods), then if\nIceModel uses the same albedo scheme as SkinModel and performs better in comparison\nagainst runoff measurements, then its improved performance is due to the inclusion of bare ice\nrefreezing.\n● We agree with this reasoning and have strengthened the text to emphasize that\nrefreezing is the primary process differentiating IceModel from SkinModel. Please\nsee Section S4.3.\n4) If SkinModel emulates RCMs then we can use the improvement of IceModel over SkinModel\n(i.e. (3) above) as a reliable indicator of the improvement to RCMs if they were to include bare\nice refreezing.\n● Our results support the conclusion that including refreezing improves runoff\npredictions by IceModel relative to both RCMs and SkinModel. However, other\nprocesses also influence model performance. We have expanded the Discussion\nsection to clarify that while our study isolates refreezing as a critical factor, it is\nnot the only one.\n● Our primary conclusion is that refreezing in bare ice reduces runoff from bare ice.\nWhile we also show that incorporating this process improves agreement between\nmodel predictions and observations, the degree to which individual RCMs would\nimprove by incorporating the processes in IceModel is beyond our scope and\ncannot be predicted. We do not believe our conclusions or their broad\nsignificance depend on assumption (4).\n● To address potential confusion, we have clarified in Section S4 that the\nimprovement of IceModel over SkinModel does not necessarily serve as a reliable\nindicator of the improvement to RCMs if they were to include bare ice refreezing.\nThis logical argument was proposed by the reviewer, but it was not a claim made\nin our article. Instead, our study documents the physical process, models its\nimpact, and deduces that RCMs would be improved by including it. However, this\nproposed improvement is one of physical realism, not a quantitative estimate of\nthe improvement to RCMs.\n● The SkinModel emulator simulations allow us to explain the spread among RCMs\nand show that it is driven by albedo differences, which is critical for interpreting\nour experimental results at RB catchment. However, SkinModel is not intended to\nserve as a basis for quantifying the improvement RCMs would experience if they\nincorporated the additional physics in IceModel. Such an assessment is beyond\nour study’s scope.\n5) If we accept that IceModel is the most accurate and the most physically realistic model then\nthe quantitative amounts of refreezing that it predicts are a reliable estimate for the actual effect\nthat bare ice refreezing has on surface mass balance.\n● We agree that IceModel provides a physically reasonable estimate of refreezing\nbut acknowledge that additional processes affect surface mass balance. We have\nexpanded the Discussion section to include a discussion of model uncertainties\nand potential improvements in future work.\nAddressing point 1 above, the authors refer to a substantial body of previous work (e.g. refs 13\nand 34) in measuring surface runoff from the internal catchments within the bare ice ablation\narea. Although the possibility that some runoff has been missed from these internal catchments\nis not ruled out in this earlier work, the team have done a very thorough job and make\na convincing case that their measurements are as good a data set as we are likely to get to test\nmodelled runoff against. Measuring the huge discharges that emerge from the ice sheet into\nproglacial streams is a very different challenge that is outlined in ref 43. Here the most\nsignificant uncertainty is whether the runoff measured in the proglacial stream actually\nrelates to the catchments being modelled. This is very hard to quantify but a decent attempt to account for this uncertainty is\nincluded through varying the modelled catchment boundaries according to different assumptions about subglacial water\nrouting given the surface and bed topography of the upstream ice sheet. So point (1) is fine.\n● We appreciate the reviewer’s recognition of the strength of our supraglacial\ndischarge measurements. We also agree that proglacial discharge measurements\ncome with inherent uncertainties.\n● To address this, we have moved our comparison with proglacial measurements at\nAK4 and LG to the Discussion section, where we explicitly discuss the associated\nuncertainties. Our conclusions remain primarily based on detailed field studies\nand numerical modeling of a physical process that we have documented\nempirically in previous publications but that remains largely unexplored in climate\nmodels.\n● To further address this, we have added a new discussion of uncertainty to the\nSupplementary Materials. Please see Section S4: Supplementary Discussion.\nAddressing point (2) above, SkinModel emulates the MAR and RACMO models very well at the\nsmall scale over a short time period (as shown in Fig1). Different albedo schemes, (RACMO,\nMAR, MODIS) are applied to the “RCM emulator” SkinModel to show that differences in albedo\nalone cannot account for the discrepancy between modelled and measured runoff. So for the\nRB catchment from July of 2015 and 2016, RACMO albedo is closest to the MODIS derived\nvalues (and the AWS measurements) whilst the MAR and MERRA-2 RCM albedos are much\nhigher (Fig S9). Even though RACMO has more accurate albedo and likely does a good job of\nestimating surface melt (as shown in Figs S5 and S6), its continued persistent\noverestimation of runoff must be due to some other process omission/error. The argument is that MAR and\nMERRA-2 only have a better match with observations because they use unrealistic albedo\nschemes. SkinModel forced with AWS albedo measurements does a slightly better job than the\n“RACMO emulator” but it is still not a good agreement to measured runoff. IceModel, driven by\nthe same AWS albedo but also including bare ice refreezing, improves the results significantly and\ntherefore the inclusion of this process is quite logically proposed as being the key point of\ndifference.\n● We appreciate this summary of our results and confirm that it accurately reflects\nour findings.\nSo SkinModel emulates the MAR and RACMO models very well at the small scale over a short\ntime period, but there is a significant discrepancy between SkinModel (MAR albedo) and MAR\nat the larger LG catchment scale over period 2009-2012 (Fig 4c). So it could be that SkinModel\nmay not actually be a particularly good emulator of MAR (or RACMO?) at larger spatial and\ntemporal scales. Why does this discrepancy arise? Albedo at different temporal and spatial\nscales may be relevant here.\n● The reviewer raises a valid point about the discrepancy between SkinModel (MAR\nalbedo) and MAR at the LG catchment scale. However, we do not expect\nSkinModel to perfectly replicate MAR at all spatial and temporal scales, nor is this\nnecessary for our conclusions. Besides, the LG catchment presents significant\nuncertainties in drainage basin delineation, making it an imperfect test for model\nemulation.\n● Most importantly, our study does not rely on SkinModel being a perfect emulator\nof MAR or RACMO at all scales. Our core conclusion—that refreezing reduces\nrunoff—is supported by direct field observations and process-based numerical\nmodeling, particularly at the RB catchment, which is the primary focus of our\nstudy. The focus on SkinModel emulation at LG overlooks this broader validation.\nTo clarify this, we have:\n● Acknowledged in the Discussion that while other factors such as spatial\nheterogeneity in albedo may contribute to discrepancies at LG, these do not affect\nour primary conclusions.\n● Emphasized that uncertainties in catchment delineation at LG make it difficult to\nuse as a definitive test for model emulation.\n● Reinforced that our main contribution is the identification and quantification of\nrefreezing in bare ice, demonstrated through field experiments and robust\nnumerical modeling, rather than the large-scale emulation of RCMs.\nWhile we appreciate the reviewer’s focus on LG, we respectfully note that this particular\ndiscrepancy does not alter our findings.\nAt the small scales shown in Fig S9, MAR and MERRA-2 have the highest albedos and RACMO\nhas lower albedo. However, Fig S11 shows that for areas of the SW sector above c.900m\n(which is most of the ablation area), RACMO albedo on average is actually higher than MAR\nalbedo. For the LG catchment from 2009-2012, RACMO overestimates runoff by 30% compared\nto a MAR overestimate of only 8%. This is despite RACMO quite probably having higher\naverage albedo for this overall catchment and time period (based on the long-term, sector-wide\ndata shown in FigS11). So the differences between RACMO and MAR are significant and are\nnot all about different albedo. This starts to mess up the assumptions that held at the smaller\nscales above, and indicates that processes other than albedo and refreezing may explain\ndifferences. We see results at the larger scales that don’t fit neatly to the patterns seen at small\nscale. In summary, point (2) holds at smaller scale but is more ambiguous at larger scales.\n● We respectfully disagree. Fig. S11 (revised Fig. S10) does not support the\nreviewer’s claim, nor did we make assumptions about albedo at any scale.\n● The July–August average albedo in LG catchment from MAR, RACMO, and\nMERRA-2 is shown in Figure 1 below. It is clear that the albedo patterns in LG\ncatchment are in fact quite consistent with those at RB catchment, with MERRA\nand MAR exhibiting high albedo, biased high relative to MODIS, and RACMO\nhaving lower albedo, in closest agreement with MODIS.\n● Therefore, the assertion that Fig. S11 “messes up the assumptions that held at the\nsmaller scales” is incorrect, as no such assumptions were made, and the results\nare consistent at both scales.\n● To address this, we have added the information in Figure 1, shown below, as a\nnew subpanel in Fig. S11 (revised Fig. S10).\nFigure 1: July-August average albedo in LG catchment from MAR3.11, RACMO2.3, and\nMERRA-2 over the period 2008-2018. RACMO values are omitted for 2008-2011 due to\ndata availability.\nPoint (3) above is challenged by the fact that for the larger catchment scale shown in Fig 4c,\nSkinModel (MODIS albedo) does just about as good a job as IceModel (MODIS albedo) over\nthese larger scales. The first tracks the upper error band on observations of runoff, the second\ntracks the lower error band and all three overlap if error bands on model outputs are\nconsidered.\nThe reviewer suggests that at larger spatial and temporal scales, IceModel and\nSkinModel (both using MODIS albedo) yield similar results, challenging the claim that\nIceModel’s improved performance is solely due to refreezing.\n● We respectfully disagree. While there is some overlap in error bands, a closer\nanalysis of biases reveals that SkinModel (MODIS albedo) systematically\noverpredicts runoff at every site, much like RACMO and (depending on albedo)\nMAR and MERRA. In contrast, IceModel (MODIS albedo) consistently predicts\nlower runoff and remains within observational error bounds. Specifically, in Fig. 4c\n(revised Fig. 5c), SkinModel runoff is above the upper error bound (+21% higher\nthan observations), while IceModel runoff falls on the lower error bound (-15%\nlower than observations). IceModel is therefore more accurate than SkinModel,\nand the 32% reduction in runoff from SkinModel to IceModel is due to refreezing,\ndemonstrating that this process is the primary difference between the two models.\n● While Fig. 4c (revised Fig. 5c) suggests that additional research questions remain\nregarding the observational runoff record at LG, it does not undermine our\nconclusions. One key factor that likely plays a role at LG is the uncertainty in the\ndefined catchment area. If the true LG catchment is larger than inferred from\navailable surface and bed topography, this would shift all model predictions\nupward, reconciling some of the discrepancy. However, this uncertainty does not\naffect our primary conclusion: IceModel predicts systematically lower runoff than\nSkinModel at all sites, including LG, and this difference is explicitly due to\nrefreezing.\n● To address this comment, we have clarified in the Discussion and Section S4 that\nwhile refreezing is a dominant factor in reducing runoff, additional uncertainties in\ncatchment delineation may also impact the comparison at LG.\nIf SkinModel does not emulate RCMs consistently then point (4) above is more uncertain and it\nis possible that the apparent improvements that IceModel shows over SkinModel may be due to\nother factors that are not explored.\n● We respectfully disagree. We believe the reviewer is placing undue emphasis on\nthe assumption that the quantitative performance improvement of IceModel over\nSkinModel should serve as a proxy for the improvement that RCMs would\nexperience if they incorporated refreezing. This overlooks the extreme complexity\nof the climate models examined in our study.\n● It is not possible to predict how a climate model will perform if it hypothetically\nadds a new physical process. However, what can be said with confidence is that\nclimate models should include physical processes that: 1) occur in reality, and 2)\nare known to be critical factors in the predictions being made by them.\n● Refreezing in bare ice is a physical process that affects meltwater retention and\nrunoff. Therefore, climate models should include light penetration and the\nsubsequent melting and refreezing that occurs in bare ice, as “anyone who has\nwalked across a freezing bare glacier surface late at night or early morning will\nattest to”, to quote the reviewer. Whether the exact quantitative improvement in\nIceModel over SkinModel translates directly to other models is beyond our study’s\nscope, but that does not diminish the fundamental need for models to incorporate\nthis process.\nThe improved performance of IceModel is evidenced by the data shown in Table S3 and is\nrepresented by the average percentage difference between modelled and measured runoff, µ. I\nthink there are other ways to assess the different model performance that does not simply take\naverages of values that are calculated over very wide ranging spatial and temporal scales. As a\nfirst pass the absolute differences should be considered rather than averaging positive and\nnegative numbers. If you do that IceModel – AWS albedo comes out best. However, once you\nstart applying some kind of weighting/significance criteria to the averages (e.g. where model\nruns over larger areas and over longer time periods have greater weighting - see attached\nfile for an example), then the assessment of best model performance changes. Excluding runs that\ndidn’t include the full 2009-2012 LG catchment run, then the best model is SkinModel with\nMAR3.11 albedo, followed by MAR 3.11. All this means that IceModel doesn’t necessarily lead\nto improved model performance which undermines confidence in points (4) and (5) above.\nThe reviewer suggests that different methods of evaluating model performance (e.g.,\nusing absolute differences instead of averaging percentage errors) may yield different\ninsights. They also propose applying weighting schemes to emphasize longer time\nperiods and larger catchments.\nWe appreciate this suggestion and have reassessed our evaluation approach.\n● When considering absolute differences, the conclusions remain\nunchanged—IceModel reduces runoff overprediction relative to SkinModel, and\nbetter matches observations, particularly when forced with MODIS albedo.\n● The weighting scheme proposed by the reviewer assigns 97% weight to LG\ncatchment, effectively eliminating the evidence from AK4, SLV1, SLV2, 660, and\nRB. Given the uncertainties in LG’s subglacial drainage pathways and\ncontributing area, it is inappropriate to allow LG to dominate the model evaluation.\nInstead, we have adopted a controlled comparison approach:\n● We compare SkinModel and IceModel under identical forcing conditions to isolate\nthe impact of refreezing.\n● We emphasize the well-constrained supraglacial runoff measurements at RB,\nwhich provide the most reliable validation data.\n● We moved the proglacial and SLV comparisons to the Discussion section to\nclarify their role as supplementary assessments rather than primary validation\nsites.\n● We expanded the discussion of uncertainty in both the main text and the\nSupplementary Materials (Section S4).\nFinally, the reviewer points out that considering models with known albedo biases (MAR\n3.11, see revised Fig. S10) or missing physical processes forced with biased albedo\n(SkinModel forced with MAR albedo) should be taken to mean “that IceModel doesn’t\nnecessarily lead to improved model performance”. We respectfully disagree, for reasons that\nare thoroughly explained in this response letter.\nAnother issue for point (5), i.e. using IceModel to try to quantify the bare ice refreezing process,\nis that there is no sensitivity testing of IceModel shown or referred to. For the paper to stand, we\nhave to be convinced that IceModel is the most accurate and that it is the most physically\nrealistic. What parameters are IceModel (and for that matter, SkinModel) most sensitive to? The\nice core data is helpful in demonstrating the improved process accuracy of IceModel over other\nmodels. However, such data is also small scale. The paper would be strengthened if we could\nsee that IceModel was calibrated using larger scale data sets that convinced us that it could be\nupscaled to the catchment and ice sheet sector scale. More specifically, are SkinModel and\nIceModel calibrated using data that is independent of the model validation data? There may not\nbe space in the paper to include all this, but there are 55 pages of Supplementary material but\nnothing included that answers these quite fundamental points.\n● No calibration was performed in IceModel—it is a physics-based model with no\ntuning parameters. Its performance is dictated entirely by physical laws and\nprescribed boundary, initial, and forcing conditions. This is why no calibration\ninformation appears in the 55 pages of Supplementary material or the Methods\nsection—if calibration had been performed, it would have been clearly\ndocumented.\n● Rather, IceModel’s uncalibrated performance at six independent sites across a\nrange of spatial and temporal scales is strong evidence that it is an accurate\nmodel with physically realistic parameter values, none of which were calibrated.\n● However, we appreciate the importance of understanding model sensitivity and\ncalibration. We agree that this is an important and often overlooked issue in ice\nsheet modeling. While sensitivity tests were not a primary focus of our study, we\nrecognize the value of such analyses and hope our findings will motivate further\nresearch in this direction. We plan to pursue this in future work as part of our\nbroader efforts to improve process representation in ice sheet models.\nTo address this comment, we have expanded the Discussion to acknowledge the\nimportance of future sensitivity testing of process-based models including RCMs.\nSorry for such a long review but I find myself torn. I like the concept of the paper and it’s great\nthat the authors have obtained data at the smaller scales that really backs up the overall idea.\nHowever, it seems to unravel at the larger temporal and spatial scales which makes the\nquantification of this effect much harder to have confidence in. It might involve too much work\nfor a revision for the authors to turn this around unless they can demonstrate how I’ve\nmisunderstood and overlooked some key points that alleviates these concerns. It will make a\ngreat paper if they can!\nWe appreciate the reviewer’s thoughtful engagement with our study and their recognition\nof its significance.\nWe respectfully disagree that our findings “unravel” at larger scales. Rather, we\nacknowledge that quantification at large spatial and temporal scales comes with inherent\nuncertainties—uncertainties that we explicitly discuss in our manuscript. The\nfundamental process we identify (bare ice refreezing) is well-supported by our field data\nand physics-based modeling, and our broader analyses demonstrate that its effect\nremains substantial across a range of spatial and temporal scales.\nWe hope these clarifications reinforce confidence in our conclusions and demonstrate\nthat our study makes a meaningful contribution to understanding the role of refreezing in\nice sheet hydrology.\nSmaller points:\nFigure 2,3 and 4 captions start with a statement that reads a bit like a paper conclusion. These\nshould be avoided and just explain to readers in very factual terms what the Figures show.\n● The statements at the beginning of Fig. 1–4 are figure caption titles, formatted\naccording to the journal’s style guidelines. These titles, along with the captions,\nprovide factual summaries of what each figure communicates. However, we have\nmade the following refinements:\n○ Edited the title for Fig. 2 (revised Fig. 3) to better reflect its content.\n○ Revised the caption for Fig. 4 (revised Fig. 5) to improve clarity.\nIf SkinModel (MAR albedo) referred to in Fig 4 is essentially the same model as MAR Emulator\n(SkinModel) referred to in Fig1, then they should be named consistently to avoid confusion.\n● We appreciate this suggestion and have revised the legends of Fig. 1 (revised Fig.\n2), Fig. 3 (revised Fig. 4), and Figure S6 to use “SkinModel (MAR albedo)” instead\nof “MAR Emulator (SkinModel)” for consistency with Fig. 4 (revised Fig. 5).\nHowever, we continue to use “emulator” in the main text to describe the role of\nSkinModel in explaining RCM differences. This ensures clarity in both figure\ncaptions and the broader discussion.\nWhy isn’t SkinModel (RACMO albedo) shown in Fig4? If it does a better job of emulating\nRACMO then this would alleviate the issue that emulators don’t transfer well across scales.\n● As detailed in the Methods section, simulations using SkinModel with RACMO\nalbedo were not performed at the SW sector scale because RACMO data were\nunavailable for the full simulation period. Additionally, the available RACMO data\nwere not adequate for model forcing as the wind speed, relative humidity, and\nsurface pressure were not available.\nWe greatly appreciate the reviewer’s detailed feedback and thoughtful suggestions.\nThese revisions have further strengthened our manuscript.\nSincerely,\nMatt Cooper\nREVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nOverall, this was a strong review. The authors adequately addressed any minor\nconcerns that I had with the original submission. In particular, I think that the slight\nreorganization that they did makes a big difference in readability and gives added\nemphasis to the model differences between SkinModel and IceModel. I also appreciate\nthat the authors added some nuance to their discussion for uncertainty associated with\nthe larger basin analyses (e.g., groundwater and subglacial hydrology uncertainty). I\nhave some additional minor points below, but no major concerns. The article is well\nwritten and a great fit for Nature Communications.\nI should have brought this up in the first review, but it became more striking to me that\nIceModel over-predicts the ice density at depth (i.e., Figure 3b at 0.5-1 m depth). Is this\nbecause light is penetrating deeper below the surface (more subsurface melt) than your\nmodel predicts? Are there implications to that for your model performance? A few\nsentences pointing this out in the discussion would be sufficient in my opinion.\nThank you for highlighting the depth-dependent bias in simulated density. To address\nthis, we have:\n1) Revised the paragraph to acknowledge the bias first (L183 in the tracked changes\ndocument).\n2) Added a short paragraph to the Supplementary Information (S4.2 Uncertainty in\nice core density measurements) noting that melt/drainage during core handling\ncould have lowered the measured densities at depth.\nThe revised sentence in the main text now reads:\n“Below ~0.4 m, IceModel over-predicts ice density by up to ~150 kg m‑3 relative to\ncore measurements, suggesting deeper light penetration and/or lateral\nmeltwater delivery may increase porosity at depth more than the model allows.”\nThe paragraph then states that, despite this depth-specific bias, the column-average\ndensity remains close to (but above) observations, and well below the canonical 900 kg\nm⁻³ assumed in many climate models.\nThis is a minor concern, but many journals are now including the references from the\nsupplementary material in the main article.\nTo the best of our knowledge, and based on this document, SI references are not\nincluded in the main text reference list at Nature Communications. However, we also\nunderstand that these types of formatting issues will be addressed by the copy editor if\nthe paper is accepted for publication.\nL69-72 – I got a little lost in this sentence, it is long with many commas. Consider\nmoving the final clause closer to “in situ records”, something like:\nWe compare meltwater runoff simulated by our model against outputs from two regional\nclimate models and a global climate reanalysis [as well as] in situ records from a\nwell-studied surface catchment on the southwest GrIS ablation zone, including:\nsupraglacial river discharge, surface ablation, and physical bare ice properties.\nRevised as suggested.\nL208-210 – I like this framing around inefficient melt. Consider adding some sort of\nphysical intuition like this to the abstract.\nThe abstract is 20 words longer than requested by the journal guidelines and we were\nunable to reduce its word count enough to accommodate this additional framing.\nL254-255 – I am not sure you need to credit the photos to yourself, it is implied.\nThe attribution is included to distinguish the aerial imagery shown in Fig. 1a from that\ncollected by co-author J.C. Ryan, which was used for surface classification of the entire\nupstream contributing area, as described in the Methods.\nReviewer #1 (Remarks on code availability):\nAs before, the code is freely available and well documented.\n—\nWe thank the reviewer for their feedback, which has greatly improved the clarity of our\nmanuscript.\nSincerely,\nMatt Cooper\nReviewer #2 (Remarks to the Author):\nThis manuscript uses field observations and surface energy balance modeling to\ndemonstrate potential liquid water storage within west Greenland ice sheet bare ice.\nThe results suggest that overestimation of runoff by regional climate models and the\nMERRA-2 global reanalysis may likely be attributable to storage and refreezing of water\nwithin the bare ice. These are important findings and hold potential for improving\nestimates of ice sheet mass change.\nI’ve evaluated the manuscript as well as the previous second reviewer’s evaluation and\nthe authors’ response. In general, I find the manuscript to be well written and supported,\nthough I have a few concerns documented below that should be addressed before the\nmanuscript should be considered further for publication.\nSpecific comments:\nRegarding R2’s point: “Fig S11 shows that for areas of the SW sector above c.900m …,\nRACMO albedo on average is actually higher than MAR albedo. … So differences\nbetween RACMO and MAR are significant and not all about different albedo. This starts\nto mess up the assumptions held at the smaller scales above, and indicates processes\nother than albedo and refreezing may explain differences.”\nAnd the authors response:\n“Fig. S11 (revised Fig. S10) does not support the reviewer’s claim”\nMy take: When I look at S10 (former S11) panels a and b that show the whole SW\nsector, I tend to agree with R2. For example, at 1500m, the MAR albedo trend line is\n~0.66, whereas RACMO is ~0.71. As such, RACMO albedo is higher than MAR, on\naverage, across the SW basin. This differs from the LG catchment (S10c) and RB\ncatchment (S8), where RACMO has lower albedo than MAR. This then supports R2’s\nclaims that assumptions at the smaller scale regarding attributing differences in the\nmodels’ runoff to albedo and refreezing may not apply to the entire SW basin. I don’t\nthink this detracts from the main points of the paper, but it should be acknowledged\n(in particular in the text ranging from L123-L134).\nWe agree with the reviewer that the RACMO albedo trend line is higher than MAR above\n~900 m in Fig. S10. To address this, we made the following revisions:\n1. Acknowledged explicitly in the main text Discussion section (L278 in the tracked\nchanges document) that factors beyond albedo and refreezing, particularly biases\nin snowline elevation, could influence interpretations of model differences at the\nSW sector scale (see our extended clarification below for why the RACMO trend\nline is related to snowline elevation).\n2. Added a new paragraph to the supplementary Discussion (Section S4.1, final\nparagraph) to clearly explain RACMO’s albedo bias as indicated by the trend lines\nin S10 and its linkage to our MODIS-derived bare ice mask. We also added a note\nin the Fig. S10 caption directing readers to this paragraph.\n3. Revised the caption of Fig. 5, directing readers to Fig. S10 for context on albedo\npatterns specifically at the LG catchment scale.\nWe also clarify several points to help the reviewer interpret Fig. S10 and the revisions we\nmade to address this:\n1. Most importantly, we do not attribute runoff differences between MAR, RACMO,\nand MERRA-2 at the SW sector scale specifically to albedo in our analysis. In fact,\nwe do not report RACMO or MERRA-2 runoff at the SW sector scale. MAR runoff at\nthis scale is used strictly to contextualize refreezing estimates from IceModel,\nspecifically highlighting that IceModel refreezing equates to ~9–15% of total (bare\nice + firn) runoff simulated by MAR. In contrast, Fig. S10 panel c shows that the\nalbedo patterns observed at the smaller scale catchments do in fact hold up at the\nlarger LG catchment scale. This is what we referred to when we stated that “Fig.\nS11 (revised Fig. S10) does not support the reviewer’s claim”, since R2’s\nargument centered on LG catchment. Their assumption that the SW sector trend\nlines could be extrapolated to LG catchment was insightful but ultimately\nincorrect, and allowed us to clarify that albedo patterns at the LG catchment scale\nare in fact consistent with patterns at the smaller-scale catchments, like RB.\n2. To understand why the RACMO trend line is steeper, recall that the values shown\nin Fig. S10 are confined to grid cells within our SW sector bare ice mask, which\nwas constructed from MODIS observations. Thus, RACMO's higher trend line\nslope reflects its documented bias toward greater snow accumulation and larger\nsnow-covered area relative to MODIS (Ryan et al., 2019, 2020). This bias creates a\ncluster of high-albedo (snow-covered) grid cells above ~1100 m that steepens\nRACMO’s trend line. Although MERRA-2 snowlines were not analyzed in Ryan et\nal., 2019, 2020, MAR snowlines were shown to be more consistent with MODIS,\nthus explaining the relative absence of high-albedo (snow-covered) MAR grid cells\nat the higher elevations in Fig. S10.\n3. RACMO’s snow-cover and associated albedo bias does not propagate into our\nprimary analyses, as neither IceModel nor SkinModel is forced by RACMO at the\nSW sector scale.\nAdditional concerns:\nL107: The authors cite Fig S4-S5 to claim an evaluation of albedo in various models.\nHowever, these figures do not show albedo directly, but rather SW . Since albedo is the\nratio of SW /SW¯, SW alone is insufficient to assess whether albedo is accurately\nrepresented. For example, we don’t know whether discrepancy in SW­ is a function of\nthe model’s albedo or rather, the amount of SW¯ that is impacted primarily by the\natmosphere.\nTo support claims about fidelity in assessing albedo, the authors should present both\nSW and SW¯ to compute albedo or revise their text to accurately state they are not\ndirectly evaluating albedo.\nThank you for identifying this oversight. Please note that Fig. S8 (Fig. S4 in the revised\ntext) does in fact directly evaluate albedo, as requested. On L107, we should have\nreferenced this figure.\nTo address this, we have now included a reference to revised Fig. S4 (formerly Fig. S8),\nand we apologize for any confusion caused by omitting this reference originally. Please\nnote that the sentence in question (L107) also refers to net radiation, which remains\nsupported by Fig. S5 and S6 (formerly Fig. S4 and S5), thus the revised text now\nreferences all three figures (Fig. S4 and S5–S6), which together directly evaluate net\nradiation and albedo across models and observations.\nThis issue is particularly evident in the authors’ claim that RACMO closely reproduces\nobserved albedo, citing Figure S4. However, a closer examination of this figure reveals\nthat although RACMO represents SW¯ and Rnet reasonably well, the agreement is\nlikely not because RACMO is getting albedo correct. Specifically, RACMO significantly\nunderestimates LW¯ during a multi-day period ~July 9-11, suggesting it\nunderrepresents cloudiness. Despite this, RACMO matches Rnet, likely through a\ncompensating error – because it has similar SW and seemingly misses this cloudy\nperiod (evidenced through lack of LW¯), the only way it can match Rnet is by having a\nlower albedo than the surface observations. In other words, more SW would be\nincoming due to lack of clouds in RACMO, and to match the SW and Rnet, it must be\nabsorbing significantly more SW radiation and thus have a lower albedo. It gets the right\nanswer for the wrong reasons.\nThis lack of critical evaluation raises significant concern that the authors have not\nclosely considered the data upon which they are basing their statements. At minimum, I\nwould suggest the authors should reexamine their evaluation framework and revise\noverly strong claims that are not directly supported by the figures or data presented.\nWe thank the reviewer for bringing this discrepancy to light. However, our original\ninterpretation remains accurate, with respect to both the 2-day period highlighted by the\nreviewer and the full 7-day experimental period, as described below. Importantly, we\nbelieve the reviewer’s perception of this discrepancy may have been inadvertently\ninfluenced by the y-axis limits, which were not held constant across panels, making the\nLW¯ discrepancy appear larger than it is. While the actual LW¯ discrepancy (relative to\nAWS values) does reach 66 W/m2 during the period in question, the 7-day mean error\nis -9.6 W/m2, whereas the mean error in SW¯ and net radiation are -8.6 W/m2 and -18.2\nW/m2, respectively.\nWe’ve included a figure below with approximately identical y-axis limits across panels\n(except panel (a) which shows SW¯), to clarify how the y-axis limits may have impacted\nthe reviewer’s perception (see panels (c) and (d) in particular). This figure differs slightly\nfrom the one in question (revised Fig. S5) to specifically highlight differences between\nRACMO values and weather station observations, showing both downward and net\ncomponents to fully clarify discrepancies (or lack thereof) in the modeled radiation\nbudget. We hope this figure (together with revised Fig. S4) clarifies our interpretation that\nRACMO has the most accurate albedo and net radiation for this time and location.\nTo address this, we adjusted the y-axis limits of Fig. S5 and Fig. S6 to reduce potential\nfor similar confusion among readers. Please note that we also adjusted the color scheme\nto be consistent with the color scheme used in other figures.\nFigure 1: Surface energy balance comparison. Hourly values of (a) downward shortwave\nradiation, (b) net shortwave radiation, (c) downward longwave radiation, (d) net longwave\nradiation, and (e) net radiation during the July 2016 field experiment from RACMO2.3p3 and\nKAN_M weather station observations.\nL108: Terminology of “net radiative heat fluxes” is strange – should probably say “net\nradiation” or “turbulent heat fluxes and net radiation” . Seems like maybe this was a\nconsequence of word cutting so the assessment of the turbulent heat fluxes and net\nradiation got merged.\nRevised as suggested.\nL129: Again, the authors have not demonstrated that RACMO has the most accurate\nalbedo. Revise.\nPlease see our response to related comments above—L129 (L143 in the revised tracked\nchanges document) states that RACMO is “the model with the most accurate albedo for\nthis time and location” referring specifically to the RB catchment during our July 2016 field\nexperiment. This statement remains accurate, as supported by Fig. S4 (formerly S8), which\nshows that RACMO is the model with the most accurate albedo relative to both AWS and\nMODIS for this time and location. Again, we apologize for omitting reference to this figure\nin the previous version, which we have now addressed.\nL148 / Figure S9: Monthly labels in the inset figure are very hard to interpret – for\nexample which line relates to which month is almost impossible to assess. I would\nsuggest the authors consider revising this figure to potentially include using a legend\nthat links the line color/style or month number.\nRevised as suggested.\nFigure 4: Referring to panel b, the ice core density profile suggests that the subsurface\nis even more porous than IceModel, and thus there’s even more potential for storage of\nmeltwater. I am curious how the authors have accounted for this discrepancy in their\nmodel or results.\nThank you for highlighting this discrepancy. To address this, we have:\n1) Revised the paragraph to acknowledge the discrepancy (L183 in the tracked\nchanges document).\n2) Added a new paragraph to the Supplementary Information (S4.2 “Uncertainty in\nice core density measurements”) noting that melt/drainage during core handling\ncould have lowered the measured densities at depth.\nThe revised sentence in the main text now reads:\n“Below ~0.4 m, IceModel over-predicts ice density by up to ~150 kg m‑3 relative to\ncore measurements, suggesting deeper light penetration and/or lateral\nmeltwater delivery may increase porosity at depth more than the model allows.”\nAs mentioned elsewhere, IceModel is an uncalibrated physically-based model, and we\nwould not expect it to perfectly explain the ice core density observations without\ncalibration. We believe the model produces a realistic subsurface density profile without\nthe risk of overfitting to field measurements, and we now acknowledge this explicitly in\nthe revised text and SI.\nL182 / Figure 5: Panels c and d both show 8 lines each that are very difficult to\ndifferentiate. Part of this is the small size of each subplot and part of this is that the\nweight of the lines and dash patterns in the legend do not equal the weight and patterns\nof the lines in the actual plots. For example, which green line is RACMO and which is\nIceModel? Which dashed/dotted line is SkinModel and which is MERRA-2?\nWe acknowledge the challenge of displaying multiple time series in a compact figure. To\ndistinguish the time series, we varied not just color, but also line style (solid, dashed,\ndotted) and line weight, with colors generated from a maximally perceptive colormap.\nWhile we cannot control the legend rendering without modifying the plotting software\nbackend, the differing line weights are clearly visible in the plot itself and map to the\ncolors and symbols in the legend.\nSpecifically, RACMO is the thicker bright green line; IceModel forced with MAR albedo is\nthe thinner dark green line. They can also be distinguished by their positions relative to\nobserved runoff—RACMO overpredicts, IceModel underpredicts—as reported in the text\nand Table S3. SkinModel is the dotted grey line, and MERRA-2 is the dashed line, as\nindicated in the legend. We’ve taken considerable care to make each line distinct and\nbelieve the figure is interpretable as presented.\nTo address this, we added a note to the Fig. 5 caption that Table S3 can be used to aid\ninterpretation of panels (c) and (d) by reference to the reported mean errors.\nL186 (and elsewhere): This is a stylistic thing, but “underpredicts by -15%” is a double-\nnegative. Suggest using the directional verb or signed value, not both.\nRevised as suggested.\nL406 onward (Method section related to climate models):\n• I suggest the authors should carefully assess which modeled products they used\nand whether the citations here are correct. For example, the authors state they used\nRACMO2.3p3, citing reference 22 (Noël 2018). The Noël paper discusses\nRACMO2.3p2 (not p3).\nThank you for identifying this. We have corrected the citations to accurately reflect the\nRACMO versions used.\n• Beyond this, I am concerned about inconsistencies in the RCMs used in this study.\nMAR data are ~15 km resolution and forced with ERA5. The RACMO data are\napparently ~11 km and forced with ERA-Interim. Ideally, the authors should be using\nmore recent RACMO data forced by ERA5 for a more fair assessment of the relative\nability of MAR and RACMO.\nWe agree that differences in forcing data and resolution are important to consider. Our\nmanuscript clearly describes the specific versions, resolutions, and forcings of the\nmodels used, allowing readers to assess these differences. We believe it is fair to\nevaluate the models as they were provided to us by the modeling teams.\n• The authors state they use both RACMO2.3p3 and RACMO2.3p2. The citation for the\np2 data is for the 1km downscaled version of RACMO2.3p2, yet the authors state\nthey used ~11km RACMO products.\nAs requested, we have corrected the citation for RACMO2.3p2 to Noël et al. (2018), which\ndescribes the ~11 km product used in our analysis.\n• Both MAR and RACMO have been statistically downscaled to ~1km resolution. Why\nare the authors not using these higher resolution products that would in theory\nbetter represent catchment-scale meltwater processes?\nUsing statistically downscaled products would conflate the impact of model physics with\nthe impact of the downscaling method, which would interfere with our study objective:\nisolating the influence of model physics—specifically radiative transfer in bare ice and\nsubsurface melt/refreezing—revealed by our field observations. While we believe\nstatistical downscaling has a useful purpose, it cannot substitute for physically based\nprocesses and may not remain valid under future climate conditions. Moreover, the\ndownscaling methods for MAR and RACMO likely differ in unforeseen ways, and no\nequivalent product exists for MERRA-2. Using downscaled inputs would also imply a\nneed to run SkinModel and IceModel at 1 km resolution, increasing runtime by ~25x from\n>1 month to >2 years. Finally, it is unclear if all required forcing variables needed for\nSkinModel and IceModel simulations are available at 1 km resolution, or if only SMB and\nits components are available. For these reasons, we evaluate all models at their native\nresolutions, as they were provided to us by the respective modeling teams.\n—\nWe thank the reviewer for their feedback, including their consideration of the previous\nsecond reviewer’s feedback, which has greatly improved the clarity of our manuscript.\nSincerely,\nMatt Cooper\nReferences\nNoël, B. et al. Modelling the climate and surface mass balance of polar ice sheets using RACMO2 –\nPart 1: Greenland (1958–2016). The Cryosphere 12, 811–831 (2018).\nRyan, J. C. et al. Greenland Ice Sheet surface melt amplified by snowline migration and bare ice\nexposure. Sci. Adv. 5, eaav3738 (2019).\nRyan, J. C. et al. Evaluation of CloudSat’s Cloud-Profiling Radar for Mapping Snowfall Rates Across the\nGreenland Ice Sheet. Journal of Geophysical Research: Atmospheres 125, e2019JD031411 (2020)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Overall, this was a strong review. The authors adequately addressed any minor\nconcerns that I had with the original submission. In particular, I think that the slight\nreorganization that they did makes a big difference in readability and gives\nadded emphasis to the model differences between SkinModel and IceModel. I also appreciate\nthat the authors added some nuance to their discussion for uncertainty associated with\nthe larger basin analyses (e.g., groundwater and subglacial hydrology uncertainty). I\nhave some additional minor points below, but no major concerns. The article is well\nwritten and a great fit for Nature Communications.\nI should have brought this up in the first review, but it became more striking to me that\nIceModel over predicts the ice density at depth (i.e., Figure 3b at 0.5-1 m depth). Is this\nbecause light is penetrating deeper below the surface (more subsurface\nmelt) than your model predicts? Are there implications to that for your model performance?\nA few sentences pointing this out in the discussion would be sufficient in my opinion.\nThis is a minor concern, but many journals are now including the references from the\nsupplementary material in the main article.\nL69-72 – I got a little lost in this sentence, it is long with many commas. Consider\nmoving the final clause closer to “in situ records”, something like:\nWe compare meltwater runoff simulated by our model against outputs from two regional\nclimate models and a global climate reanalysis [as well as] in situ records from a\nwell-studied surface catchment on the southwest GrIS ablation zone, including:\nsupraglacial river discharge, surface ablation, and physical bare ice properties.\nL208-210 – I like this framing around inefficient melt. Consider adding some sort of\nphysical intuition like this to the abstract.\nL254-255 – I am not sure you need to credit the photos to yourself, it is implied.", "reviewer_3": "Review of: Greenland Ice Sheet runoff reduced by meltwater refreezing in bare ice\nSubmitted to: Nature Communications\nReviewer: Benjamin Hills\nSummary.\nCooper et al. observe that meltwater runoff measured in the Greenland ablation zone is\nlower than would be predicted by climate models. They hypothesize that the reduction in\nrunoff is associated with meltwater retention and refreezing in a near-surface layer of\nweathered ice, much like the relatively well-studied water retention in firn (Harper et al.,\n2012). They test their hypothesis with a pair of 1-dimentional thermodynamic numerical\nmodels they developed. They find that the model which considers light penetration and\nsubsurface melting is a better representation for their measurements of meltwater runoff.\nThe article is well written and a great fit for Nature Communications after the minor\nrevisions I suggest below.\nGeneral Comments.\nIn my opinion, the difference in model output between SkinModel and IceModel is one of the\nmore important results in this body of work. If you agree, I believe that some reorganization of\nthe Results section could help emphasize that result to all readers. Currently, Results is broken\nout into subsections on: fieldwork, model overestimation of meltwater (SkinModel),\nexplanation of overestimation (IceModel), upscaling to SW Greenland. The separation between\nexplanations of SkinModel and IceModel could make it hard for the reader to appreciate what\nyou are doing here (even though I do like your clarifying statement on L137-139). I would\nconsider reorganizing as follows:\n- Intro (pretty much as you have it)\n- Results (field results only)\no You could consider moving Figs 2 and 3 up here, although I understand why you\nwant to lead with what you now have as Fig 1.\n- Analysis\no Bring your explanation of SkinModel and IceModel\no Compare to climate model output at the local scale (currently Fig 1)\no Compare to climate model output at the regional scale (currently Fig 4)\n- Discussion/Conclusions (pretty much as you have it)\nPerhaps this organization does not work well for NatureCommunications but that is how I\nwould write it t"}, "author_response": "We thank the reviewer for their feedback, which has greatly improved the clarity of our\nmanuscript.\nSincerely,\nMatt Cooper"}, {"round": 3, "reviewer_comments": {"reviewer_3": "I thank the authors for addressing my previous comments and have no further concerns."}, "author_response": ""}]

Decision letter