Greenland Ice Sheet runoff reduced by meltwater refreezing in bare ice

Matthew Cooper (  guycooper@ucla.edu ) Paci c Northwest National Laboratory https://orcid.org/0000-0002-0165-209X Laurence Smith Brown University Åsa Rennermalm Department of Geography, Rutgers, The State University of New Jersey Kang Yang University of California, Los Angeles (UCLA) https://orcid.org/0000-0002-7246-8425 Glen Liston Colorado State University Jonathan Ryan Brown University https://orcid.org/0000-0002-9432-1512 Dirk van As Geological Survey of Denmark and Greenland Lincoln Pitcher University of Colorado, Boulder https://orcid.org/0000-0001-8624-9760 Clément Miège University of Utah Sarah Cooley Department of Geography, University of Oregon Brandon Overstreet Department of Geology and Geophysics, University of Wyoming


Main
Greenland Ice Sheet (GrIS) mass loss raised global sea level 10.8±0.9 mm between 1992-2018 and climate models forecast an additional 70-130 mm by year 2100 1 . Most of this recent mass loss can be attributed to increased meltwater runoff from the ablation zone [2][3][4][5][6] .
Within this critically important zone, winter snowpack melts entirely each summer exposing dark, bare glacier ice that absorbs up to three times as much sunlight as bright snow 17 . Increasing temperatures and reduced summer snowfall have exposed larger areas of bare ice in recent decades 2,3,18 , driving enhanced surface melting in the ablation zone [2][3][4] . Understanding the fate of meltwater from Greenland's growing bare-ice zone is therefore critical for accurate modeling of sea levels 16 .
Climate models track energy and mass budgets on the ice sheet surface and are the primary tools available to predict future ice sheet runoff 13,19 . However, there is extensive evidence that climate models overestimate runoff from bare ice surfaces [7][8][9][10][11][12] . In Greenland's melt-intensive southwest sector, river discharges measured at the ice sheet periphery reveal up to 67% less annual meltwater release to surrounding oceans than climate model calculations 7,[9][10][11] . Supraglacial lakes that form on the ice sheet surface fill at slower rates than predicted by climate models 12 , and direct measurements of supraglacial runoff are overpredicted by +21-58% during peak summer melt conditions 8 . Climate models overestimate surface melt by +14-40% relative to satellite laser altimetry measurements 20 , by -17% (at one ablation-zone site) to +10-43% (at four other ablation-zone sites) relative to point ablation stake measurements 21,22 , and overestimate ice sheet mass changes by +21-47% relative to GRACE satellite gravity retrievals 23,24 . These discrepancies are not explained by errors in modeled surface energy balance components, which closely correlate with in situ meteorological observations 8,21,22 .
An explanation for why climate models currently overestimate ice sheet runoff may be their conceptual treatment of bare ice surfaces as an impervious, high-density substrate with no capacity to retain water [13][14][15] . Accordingly, runoff produced on bare ice is instantly credited in its entirety to sea level 14 , despite growing field reports of non-trivial meltwater retention on or within bare ice 7,8,[25][26][27][28][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 spaces 30,31 . Because bare ice comprises >80% of the ablation zone by surface area and generates >85% of meltwater runoff 2 , even small amounts of meltwater retention or refreezing would constitute a substantial overlooked component of the ice sheet surface mass balance and explain consistent reports of climate model runoff overestimation from bare ice surfaces 7,8,12 .
Reconciling climate model calculations with real-world runoff measurements is difficult because direct measurements of runoff on the ice sheet surface, prior to interference from englacial and subglacial processes, are extremely rare 7,11 . Here, we pair numerical surface mass balance modeling with in situ measurements of meltwater runoff 8 and bare ice properties 26,32 from a well-studied surface catchment on the southwest GrIS bare-ice ablation zone. We then reconcile discrepancies between modeled and observed meltwater runoff by quantifying retention and refreezing processes in bare ice over the entire southwest GrIS ablation zone, where the majority of GrIS surface mass loss originates 2,3 .

Climate models overestimate meltwater runoff on the bare-ice ablation zone
Rio Behar is a supraglacial river on the southwest GrIS ablation zone that drains an ~60-63 km 2 catchment, depending on year 8,12,33 ( Figure S1). Our previous field measurements during July 2015 found that regional and global climate models overpredicted Rio Behar runoff by +21-58% 8 . During a 6-13 July 2016 field experiment we revisited the site to collect a seven-day record of 168 consecutive hourly discharges 34 , concurrent with three-hourly ice surface lowering measurements from a network of ablation stakes installed near the gauging site (Methods and Figure S2). Shallow ice cores collected at the same site in July 2016 revealed that the bare-ice surface was porous and saturated with meltwater, with an average bare-ice density of just 690 kg m -3 within the upper 1.1 m of ice 26 . Here, these measured discharges (hereafter 'runoff'), surface lowering rates, and physical ice properties are compared with simulated meltwater production and runoff from two regional climate models, a global climate reanalysis model, and a numerical model of spectral radiation and thermodynamic heat transfer developed for bare glacial ice 35 (Methods).
Consistent with previous studies 7-12 , we find that climate model values of runoff overpredict measured meltwater runoff ( Figure 1). By the end of the 6-13 July 2016 field experiment, climate model runoff ranges from -13% lower (MERRA-2) to +53% higher (RACMO2.3) than observations, with all but MERRA-2 exceeding observed runoff values ( Figure 1). Among the climate models examined here, RACMO2.3 most closely reproduces observed albedo and net radiative heat fluxes, yet severely overpredicts runoff despite good representation of these critical surface energy balance components ( Figure S4 and Figure S5).
To explain these discrepancies between observations and climate models, we developed a process-based numerical model of ice sheet meltwater runoff that we call 'SkinModel' (Methods). SkinModel represents the ice sheet surface as an infinitely thin, hydrologically impermeable 'skin' layer of high-density (900 kg m -3 ) ice with zero heat capacity, emulating the assumption used in climate models for bare glacial ice (Methods Equation 1) 14 . SkinModel is forced on an hourly timestep with meteorological observations recorded by the PROMICE/GAP KAN_M surface-based automatic weather station 36 located ~2 km from the Rio Behar catchment. To first establish whether SkinModel effectively emulates climate model meltwater runoff predictions, we force SkinModel with surface albedo values output from each climate model (rather than KAN_M albedo values), and find its predictions are virtually identical to climate model predictions (Figure 1; dashed vs solid lines). Note that all model forcings besides albedo are kept consistent across these emulator simulations, signifying that differences in ice surface albedo alone drive virtually all differences between climate model predictions of runoff shown in Figure 1.
Next, we force SkinModel with values of albedo recorded at KAN_M, and find that modeled runoff is +43% higher than measured runoff ( Figure 1; green solid line). Weather station forcings thus yield runoff predictions much like RACMO2.3, the model with the most accurate albedo for this time and location. Similar results are found when the same models are tested using the same methods against our earlier field measurements from 2015 8 ( Figure S6).
Even larger runoff overestimation is simulated if calibrated MODIS satellite albedo values 37 ( Figure S7) are used as model forcing. This reveals that among the climate models examined here, accurate representation of ice surface albedo leads to overpredicted runoff, whereas accurate predictions of runoff are a fortuitous result of unrealistically high albedo ( Figure 1).

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 'IceModel' (Methods). IceModel updates an earlier model of spectral radiative and thermodynamic heat transfer in glacier ice 35 with a field-calibrated constraint 32 on shortwave radiation absorption enhancement 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 content (Methods Equation 3), informed by measurements of physical ice surface properties collected during a multi-year series of field campaigns on the GrIS bare-ice ablation zone 7,8,11,26,32,33,38 . A critical feature of IceModel is that it allows sunlight to penetrate bare ice, thus allowing melt within the ice subsurface 39 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 isothermal ice column ~1.2 m thick that stores latent heat in the form of liquid meltwater, a phenomenon independently validated by our field observations of saturated porous ice extending at least 1.1 m below the ice surface 26 . Lateral transport of meltwater generated within isothermal bare ice is constrained by its low horizontal hydraulic conductivity  37 , observed ocean-going meltwater runoff is overpredicted by +51% (Figure 4c). In contrast, when IceModel is forced with observed surface albedo, runoff is +5% higher than observations, cumulative over four years, with upper and lower estimates contained entirely within ±15% measurement uncertainty. This +5% overprediction of ocean-going runoff is consistent with observations of englacial and subglacial water storage in this region of the ice sheet 27,28 , an additional ablation-zone retention sink not represented in current climate models.
Our catchment-scale conclusions are further validated by comparison with discharge measured over seven years in the Akuliarusiarsuup Kuua River's northern tributary (AK4) ( Figure S8) and with satellite estimates of supraglacial lake volume (SLV) infilling rate (a proxy for cumulative meltwater runoff 12 ) at two sites near the Rio Behar catchment (Figure 4b). For these SLV sites, climate model runoff is +70% higher than satellite-derived runoff, whereas IceModel runoff closely tracks SLV observations (Figure 4d).
All of this suggests that the conversion of energy (input) to mass (output) is less efficient in bare ice than currently thought, thus decreasing the true export of meltwater runoff from the ice sheet surface to surrounding oceans. Once refrozen, meltwater must thaw again to become runoff, an additional energy expenditure not currently accounted for in climate models. To demonstrate the magnitude of this unrepresented energy sink, IceModel simulations indicate an annual average runoff reduction of 11 Gt a -1 in Greenland's southwest sector alone due solely to refreezing of subsurface liquid meltwater in bare ice during July and August. This bare-ice refreezing estimate is equivalent to 18% of the ~63 Gt a -1 annual refreezing in snow and firn predicted by MAR3.10 and RACMO2.3 for this sector of the ice sheet, or 10% of the ~116 Gt a -1 annual runoff. This estimate is conservative due to our use of minimum plausible bare-ice extent (Methods) and exposure (July-August) in its calculation, but is commensurate with current lower-bound estimates of ice sheet runoff uncertainty 43,44 . Longer periods of below-freezing air temperatures during seasonal transitions from melt to freeze-up similarly consume melt energy, potentially increasing runoff retention on bare ice up to 23 Gt a -1 (~20% of sector runoff) when considering June-September as an upper limit on seasonal bare ice exposure.

Implications for sea level rise predictions
Our comparison of in situ meltwater runoff measurements collected on Greenland's ablating bare-ice surface with runoff predictions from a suite of global and regional climate models finds that the latter overestimate runoff, here by up to +53%. The two models that do match observed runoff with reasonable accuracy (MERRA-2, MAR3.10) do so only incidentally, by overestimating albedo (i.e., by error cancellation). The underlying reason for this consistent tendency of climate models to overestimate runoff generation in the ablation zone 7-12 has not previously been identified, signifying a critical gap in predictive capacity concerning future ice sheet surface mass loss. Here, we show that over-prediction of runoff by climate models is explained by nocturnal refreezing of meltwater generated in the bare-ice ablation zone, a process observed on mountain glaciers 45,46 that has eluded detection in Greenland and is not currently included in climate model predictions of ice sheet runoff contributions to sea level rise 14 , nor their uncertainty 43,47 .
Regional climate models are the primary tools used to estimate the amount of runoff exported from Greenland's ablation zone to the global ocean (and the only tool for predicting future ice mass changes), but these models currently lack representation of the important nocturnal refreezing/retention process described here. Refreezing in snow and firn has potential to delay sea level rise by ~10-17 Gt a -1 at the centennial timescale 30 , yet refreezing in bare ice, currently set to zero in regional climate models 14 , may already exceed these levels, with unknown response to continued surface warming and associated expansion of bare ice areas 2 .
With nearly all (85-93%) of Greenland's meltwater runoff currently sourced from bare ice, these findings suggest the need to consider mechanisms of meltwater retention within Greenland's bare-ice ablation zone, in addition to snow and firn processes, to accurately inform future sea level predictions.    predict +5% higher runoff than observations. (d) Climate model runoff is +70% higher than satellite-derived runoff estimated from SLV infilling rate 12 , but are reconciled by IceModel simulations.

Field datasets of ice sheet surface discharge, ablation rate, and ice density
Hourly Rio Behar catchment ( Figure S1 sub-hourly measurements and associated quality measures were recorded. These raw data were converted to channel flow rate [m 3 s -1 ] in post-processing following data quality-control workflows optimized for the supraglacial environment 8,34 . The sub-hourly measurements were used to estimate ~±15 m 3 s -1 measurement uncertainty (error bars in Figure 1 representing one standard deviation in sub-hourly discharge).
Measurements of ice surface elevation change were provided by a network of bamboo ablation stakes we installed at our field camp (67.049 o N, 49.022 o W), and by pressure transducer ice ablation assembly installed on the PROMICE/GAP KAN_M automated weather station located ~2 km from the Rio Behar catchment 36 . Twelve stakes were installed within an area covering ~0.5 km 2 . 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 ( Figure S2). Stakes were drilled 3 m deep into the ice and allowed to freeze-in for 24 hours prior to initiation of ablation stake measurements recorded at 3hour intervals continuously from 12:00 on 6 July 2016 to 23:00 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 occurred). Prior to each measurement, an 8-inch square wooden ablation board was placed at the base of the stake and oriented to true north. This board operated as a datum from which the stake height was measured. Cumulative changes in stake height were converted directly to ice surface elevation change for comparison with simulated melt rates.
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 reference 26 .

Catchment boundaries and surface classification from satellite and airborne datasets
Catchment topography and surface classification was provided by 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 vehicle 49 ( Figure S1). These aerial images were used to reconstruct the ice sheet surface topography using Agisoft PhotoScan Pro stereo-photogrammetry software and to perform surface classification of snow, water, and bare ice with a k-Nearest Neighbors algorithm yielding 2.7% snow cover during the 6-13 July 2016 field experiment, and 6.5% for the 20-23 July 2015 experiment 8 ( Figure S1). WorldView-1 and WorldView-2 satellite imagery and associated high resolution stereo-photogrammetric digital elevation models were used to delineate the Rio Behar contributing catchment area following methods in reference 8 . Briefly, digital elevation based methods were supplemented with manual delineation of surface stream networks, flow direction, and channel heads following reference 50 .
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 that exist within the catchment boundary were removed from this lower catchment area estimate ( Figure S1 and Figure S3). Outer channel heads (initiation points of channels that drain away from the catchment) yields an upper maximum estimate of catchment area. The optimal "best guess" catchment area was delineated by tracing the inner and outer channel heads in the highresolution WorldView imagery and adding back areas that undoubtedly flow into the catchment and subtracting areas that undoubtedly flow out of the catchment. This approach yields contributing catchment area confirmed by actively flowing water tracks, a further improvement to digital elevation model-based methods 51 . Identical methods were used to delineate the supraglacial lake contributing catchment areas for satellite lake volume (SLV1 and SLV2) retrievals, provided by reference 12

SkinModel and IceModel description
SkinModel and IceModel solve the zero-dimensional surface energy balance 58 :

Model simulations for Rio Behar and Southwest Greenland
IceModel and SkinModel were solved on an hourly timestep at 100 m horizontal grid spacing for the experimental periods 6-13 July 2016 and 20-23 July 2015 for the ~60 km 2 Rio Behar catchment, with model input provided by the KAN_M automatic weather station 36 . following references 2,43 to isolate bare-ice areas i.e., to remove firn/snow-covered grid cells and permanent land-surface. The 75% bare-ice frequency threshold is more restrictive than the 50% threshold used in the GrIS SMB Intercomparison (GrSMBMIP) 43 , and yields a conservative 37,175 km 2 bare-ice surface area for the Southwest sector, 22% lower than GrSMBMIP.

Data availability
The data required to reproduce all figures in this manuscript are available without restriction from https://github.com/mgcooper/icemodel. Supraglacial river discharges for the Rio Behar catchment 64

Code availability
The numerical models and code required to reproduce all figures in this manuscript are available without restriction from https://github.com/mgcooper/icemodel.