Predicting mean depth and area fraction of Antarctic supraglacial melt lakes with physics-based parameterizations

Domain: earth_science
Published: 2025-07-15
Editorial outcome: major_revision

Abstract
Despite the importance of supraglacial melt lakes to the future evolution of polar ice sheets, they are not represented in current large-scale climate and ice sheet models. In this study, we use ICESat-2 satellite surface elevation measurements to show that roughness on the Antarctic Ice Sheet surface is largely self-affine. Estimation of ice sheet surface roughness statistics then enables the development of a set of simple mathematical expressions parameterizing the average supraglacial melt lake area fraction and lake depth from statistical fitting of large simulation ensembles of water flow on random, self-affine surfaces. These parameterizations provide predictions that are generally consistent with observations, with some exceptions. Finally, we predict that on large portions of Antarctic ice shelves supraglacial lakes are likely to, on average, stay less than one meter deep and occupy less than 40% of the ice area, absent changes in ice shelf surface roughness.

Introduction
Over the past three decades, observations have shown widespread supraglacial melt lakes alter the latent and sensible heat fluxes through surface melting and supraglacial lake formation on Antarctic ice air-ice sheet interactions12,13. shelves1, driven by atmospheric warming2. Atmospheric temperatures Meltwater from supraglacial lakes may also fill surface fractures in over Antarctica are expected to continue increasing in the future, the ice sheet, causing them to propagate deeper9,14–16. Incorporating along with an increasing prevalence of melt lakes on the ice sheet supraglacial melt lakes into large-scale models is thus vital in properly surface3–6. Supraglacial melt lakes have been implicated in the thin- simulating the potential interaction between surface melt and ice sheet ning and collapse of floating ice shelves in Antarctica, with the pro- fracturing and calving. Supraglacial melt lakes are currently omitted minent example of the Larsen B ice shelf collapse during the from large-scale models because they are typically small (101–103 2001–2002 melt season7–9. In the days before the collapse, many meters across)17 compared to the grid spacing in climate (104–105 supraglacial melt lakes formed on the ice shelf surface, which is meters) and ice sheet models (103–104 meters). Existing numerical hypothesized to have caused hydrofracturing across the ice shelf and models of supraglacial hydrology are computationally expensive to its eventual disintegration7,9,10.

Ice sheet surface melt affects the surface energy balance of glacial run over entire ice sheets at decadal or longer timescales11,18,19. Fur- surfaces. Darkening of ice by water (i.e., decreasing surface albedo) thermore, no physically-based theories or empirical parameterizations causes increased absorption of incoming shortwave radiation by the exist that capture the area-averaged effects of supraglacial melt lakes ice and snow. Ice sheet surface melting initiates a positive feedback, on the surface energy balance and ice sheet fracturing. This study aims whereby a lower surface albedo allows for more shortwave absorption, to remedy these gaps through the development of supraglacial melt further enhancing meltwater production and deepening lakes. After lake parameterizations that can be implemented in large-scale models.

In this study, we develop simple equations parameterizing the refreezing, albedo is permanently altered from its original value due to spatially averaged size of supraglacial melt lakes, which can be added snowpack and firn density changes from meltwater saturation, which to large-scale models at negligible computational expense. The core promotes future melt lake formation at these locations11. Additionally, idea behind these parameterizations derives from prior theoretical advances from percolation physics20,21 and terrestrial hydrology studies which find robust statistical relationships between lake size and parameters describing surface roughness (on land or otherwise)22. In this work, we define this surface roughness through the self-affinity of the glacier surface. Self-affinity is a fractal property exhibited by most of the Earth’s surface22,23. This property, quantified by the Hurst exponent, measures how the vertical scale of topography varies across horizontal scales23.

We start by using satellite altimetry measurements to show that surface roughness is self-affine over the vast majority of the Antarctic Ice Sheet. We quantify relevant roughness parameters for the grounded and floating portions of the ice sheet. We then develop a set of simple mathematical parameterizations through a large Monte Carlo ensemble of physics-based simulations of supraglacial meltwater flow over self-affine surfaces. We validate these parameterizations by comparing melt lake area fraction and depth predictions to observations from satellite imagery on two Antarctic ice shelves. Finally, we estimate the maximum attainable mean melt lake depth and area fraction for the current roughness characteristics of the Antarctic Ice Sheet. The methods used throughout this study are described in more detail in the “Materials and Methods” section. In short, this study confirms that ice sheet surfaces, like other geomorphic surfaces, are self-affine and provides a set of supraglacial melt size parameterizations that can be implemented in large-scale models.

Results
Self-affinity of Antarctic surface roughness from ICESat-2
Self-affinity describes a repeating pattern within an object, sequence, or surface that occurs over a wide range of scales (also referred to as “fractal” or “self-similar”) (Fig. 1). Analysis of land surface elevation measurements shows Earth’s surface roughness is largely self-affine23,24, extending even to the land surface beneath ice sheets25. The self-affinity of Earth’s land surface has been used to develop elegant predictions for the size distribution of terrestrial lakes filling depressions on this surface22,24,26. In this section, we investigate whether the Antarctic ice sheet surface also exhibits such self-affinity, and use satellite altimetry measurements to quantify two statistical roughness properties of the ice sheet surface: Hurst exponent (H) and standard deviation (\sigma).

We estimate surface roughness properties across the Antarctic ice sheet, using the entire Antarctic ICESat-2 ATL06 land-ice elevation data catalog from June 2021 to June 2022 (described further in “Materials and Methods”). This one-year period includes four repeats for each track27, and we found little temporal variation in surface roughness parameters over that period. All sub-tracks with calculated Hurst exponents less than zero or greater than one are discarded, since they cannot be reliably described as self-affine (i.e., outside the acceptable confines of the range for the Hurst exponent). Typically, such cases are related to spurious elevation data causing large deviations in the power spectral density. Additionally, any linear regressions with R2 values less than 0.7 are discarded to ensure that the only sub-tracks considered are those in which the data strongly indicate roughness is self-affine. Of the 1.8 million sub-tracks that were analyzed, 7.3% of the sub-tracks were discarded during the analysis process. 7.1% were discarded due to poor R2 values when fitting, and the remaining 0.2% were discarded due to the calculated Hurst exponents being outside the acceptable range of [0,1]. As a result, over 92% of all the subtracks analyzed exhibit the qualities of a self-affine surface. Consequentially, we conclude that the Antarctic Ice Sheet surface is largely self-affine.

Figure 2 also shows little spatial variation of the Hurst exponent across the ice sheet. The average Hurst exponent across the continent is \approx 0.41 with 95% of Hurst exponents estimated to be between 0.35 and 0.48. Floating ice has an average Hurst exponent of 0.47, and grounded ice has an average Hurst exponent of 0.41 (Supplementary Fig. 4). This difference between the mean of the two distributions is statistically significant as assessed by Welch’s t-test, with a reported P value of 0. With higher mean Hurst exponent, ice shelf surfaces have larger scale roughness features, since basal melt rates tend to be higher on steep small-scale basal roughness features28,29 and the buoyant flexure of floating ice filters small scale roughness at the base from surface expression. The Hurst exponent is less than previously reported Hurst exponents for subglacial topography in Greenland as observed via ice penetrating radar in25, which has a mean of 0.65. This indicates either a consistent difference between roughness in Greenland and Antarctica, or that processes related to snow redistribution and densification may influence surface topography, which creates smaller-scale topographical variability at the ice sheet surface than at the base30.

The standard deviation, \sigma, quantifies the average amplitude of roughness along the ice sheet surface within each ICESat-2 sub-track and is plotted in Fig. 2b for the whole of the Antarctic ice sheet surface. \sigma mainly varies over grounded ice due to the varying subglacial bed topography, over which ice flows. Greater roughness is notable over subglacial topographic features including the Trans-Antarctic mountains, Palmer Island, and Roosevelt Island. 74% of sub-tracks have relatively smooth ice sheet surfaces with \sigma < 10 m. Grounded ice has an average \sigma of 13.0 meters, which is considerably rougher than floating ice with an average \sigma of 4.7 meters (Supplementary Fig. 5).

We also investigated whether the estimated roughness parameters are biased by the track orientation or “capture” angle of the ICESat-2 satellite. We found that the Hurst values do not vary much based on the orientation angle of the satellite, and most of the measurements center around the value of \approx 0.4, which the average difference in Hurst exponent values at intersection points being 0.0218. For \sigma, the values between the two different sets of capture angles were also found to center around a single value. Details of this crossover analysis can be found in the supplementary material (Supplementary Fig. 1).

Parameterizations for mean lake depth and area fraction on self-affine surfaces
In the prior section, we confirmed that roughness on nearly the entirety of the Antarctic ice sheet surface is self-affine. Since water gathers in surface depressions, the statistics of roughness on the surface can be used to determine the statistics of lake size on that surface, assuming that the surface topography does not change temporally22,31. While ice sheet surface melt may deepen depressions in which melt lakes form over many years, at first order, the pre-existing topography plays a strong role in determining melt lake size as validation with observations confirms in the validation section. Prior work in the statistical physics community has tackled the related problem of percolation on uncorrelated potential surfaces20, producing analytical predictions for the size of clusters (melt lakes in our case) using mean-field mathematical approaches. However, such analytical approaches have been shown21 to be intractable for strongly correlated surfaces (H > 0; i.e., self-affine) and surfaces that are not completely flooded (i.e., at their maximum capacity for water retention), such as those commonly found on Earth. For self-affine and partially flooded surfaces, numerical methods are necessary to calculate the statistics of lake size.

In this section, we use a cellular automaton (details described in “Materials and Methods”), a discrete algorithm that divides a larger surface into state-dependent cells, to simulate the progressive development of lakes on randomly generated self-affine surfaces. We run this cellular automaton for thousands of randomly generated self-affine surface, as a function of \sigma, H, and meltwater available for ponding (referred to as “supply” hereafter). From these numerical results, we calculate two melt lake size metrics, average area fraction (F_l) and average supraglacial lake depth (w_l). Once we have a set of model results encompassing the full parameter space, we seek to find a simple mathematical relationship between each of the melt lake size metrics, and the three input parameters: surface roughness Hurst exponent (H), and amplitude (\sigma), and the averaged depth of meltwater supplied to the surface (w_s). However, first, we fit a relationship between the surface roughness parameters and the averaged depth of water over the entire surface when it is at its retention capacity, which we call w_d^* . This quantity is frequently referred to as the “maximum depression storage” in terrestrial hydrology32. We calculate the maximum depression storage in a separate set of simulations using the full range of \sigma and H by provided an excess of water supply for the generated surface and extracting the mean depth across the entire surface. We then fit a simple mathematical expression that can determine the maximum average water depth of any self-affine surface at its maximum water retention capacity utilizing only surface roughness parameters:

w_d^*(H, \sigma) = \sigma(0.2 - 0.12H^{0.6}) (1)

Consistent with prior studies in hydrology and percolation physics21,33, this equation predicts a linear increase (with a similar coefficient) in retained water depth with roughness amplitude and a sub-linear decrease in retained water depth with Hurst exponent.

We use this maximum depth to define S, a dimensionless “supply ratio” metrics, the relative difference throughout the parameter space

S(H, \sigma) = w_s / w_d^*(H, \sigma) (2)

Defining this ratio normalizes differences between surfaces with strongly different retention capacities, and thus greatly improves the parameterization fit to numerical simulations over a wide range of roughness statistics.

We fit spatially averaged supraglacial melt lake size metrics to the topographical Hurst exponent (H), roughness amplitude (\sigma), and supply ratio (S). We know that average lake depth (w_l) should be linearly dependent on \sigma since the depth of lakes is set by the depth of depressions on the surface which are described by the roughness amplitude. The numerical results from the cellular automaton confirmed this linear dependency. By definition, the average horizontal extent of vertical topographic relief on a self-affine surface does not depend on the average roughness height22. Therefore, we expect that the area fraction of supraglacial melt lakes (F_l) should have no dependence on \sigma. The numerical results from the cellular automaton also confirmed that F_l is not dependent on \sigma despite the slight variation in Fig. 4(c) across the \sigma space, which may be due to the sample size of surfaces leading to spurious residuals from the individual surfaces. With these results in mind, we were able to use standard least squares regression tools to fit a set of two parameterizations that describe a relationship between each melt lake characteristic (w_l and F_l) and the three parametric inputs (\sigma, H, S). The form of these equations is chosen to enforce certain basic physical constraints (described in “Materials and Methods”). We found the following parameterizations for supraglacial melt lake size metrics:

w_l(H, \sigma, S) = 0.6\sigma erf(67S) * (1 - 0.41H^{0.6}) (3)

F_l(H, S) = 0.13 erf(55S) * (1 - 0.13H^{1.3} + 2 erf(0.08S)) (4)

These parameterizations all have an R2 value greater than 0.82 compared to the generated numerical results from the cellular automaton. Figures 3 and 4 show the numerical predictions for melt lake metrics, the predictions from parameterizations, and the relative difference between the numerical results and parameterizations. The relative difference between the numerical and parameterized predictions is generally less than 25% for the lake depth parameterization and less than 10% for the area fraction parameterization. For both melt lake size metrics, the relative difference throughout the parameter space falls under 10% when the Hurst value is equivalent to 0.4, i.e., the mean Hurst exponent of the Antarctic ice sheet surface.

Two random self-affine surfaces with the same Hurst exponent and standard deviation can be quite different from each other and produce fairly different supraglacial melt lake distributions. Thus, when applying these parameterizations to a particular self-affine ice sheet surface area, they may depart from the melt lake statistics due to random surface variation alone. In the supplementary material, we plot the root-mean-square deviation between the parameterizations and each realization from the numerical results, to quantify the expected natural variation among randomly generated self-affine surfaces (Supplementary Fig. 3). Since we have generated enough random self-affine surfaces in our Monte Carlo ensemble to reach statistical convergence in mean melt lake size metrics, the resulting parameterizations represent the typical size statistics expected for a given ice sheet area. Therefore, these parameterizations are useful for efficiently predicting spatially averaged melt lake metrics over large ice sheet areas. To predict the exact location and depth of melt lakes over a surface with known topography, a high-fidelity model of ice sheet hydrology should be used e.g., 19, but at considerably greater computational expense which may be infeasible over larger spatial and temporal scales.

Validation of lake parameterizations with observations
After developing a set of physics-based parameterizations, we aim to validate the parameterizations against observations of melt lakes in Antarctica. To do so, we focus on two ice shelves with regular summer surface melt in the past 10 years: the Larsen C ice shelf in the Antarctic Peninsula and the Amery ice shelf in East Antarctica. To provide inputs to the parameterizations, we use measurements of H and \sigma for these ice shelves from the ICESat-2 analysis in the first Results subsection, and meltwater runoff supply as predicted by the RACMO regional climate model (described in “Materials and Methods”). We aim to validate our developed parameterizations by computing a range of mean melt lake depth and area fraction predictions for the Amery and Larsen C ice shelves, with inputs taken from the ICESat-2 derived topographical parameters (H and \sigma) and the simulated meltwater runoff supply from the RACMO v2.3p2 model34. Due to the difference in time steps and spatial grid of the observational products and predictions, we include all ICESat-2 derived topographical parameters and RACMO seasonal accumulated runoff over the regions of each of the two ice shelves where melt lakes were observed in ref. 35. From RACMO meltwater runoff, we calculate the accumulated mean water runoff depth over the region where melt lakes are observed from the beginning of each melt season and then reset the accumulated supply to zero at the end of each melt season under the assumption that refreezing occurs. Additionally, we assume that the surface roughness has no changes and remains static over time. This results in a distribution of monthly predictions calculated from equations (3) and (4), with inputs taken over the area where supraglacial lakes are observed in Landsat, and including: surface roughness parameters estimated in the first Results subsection, and surface melt supply from RACMO. We compute the interquartile, 5th, and 95th percentiles of this distribution in each month over the years 2014–2019. This prediction assumes that the surface is impermeable and that no refreezing occurs during the melt season. Additionally, we assume that topographical parameters measured by the ICESat-2 mission in 2021–2022 apply to predictions over the 2014–2018 period when observations are available due to the relatively small temporal variations in these parameters during the ICESat-2 mission. Future validation exercises could use new melt lake observations derived directly from ICESat-236 to compare against predictions made from ICESat-2 topographical parameters. Here we use the 35 dataset to match the period covered by existing RACMO simulations. All further details of these validation exercises are detailed in the “Materials and Methods” section.

In Fig. 5, we plot observed (black dots) and predicted (blue shading) melt lake depth and area fraction at Amery and Larsen C ice shelves. Observed mean melt lake depth on both ice shelves (Fig. 5B) generally falls right within the 25%–75% range at the center of the distribution of predictions, though the melt season in RACMO seems to occur consistently later than the observed appearance of lakes (by 1–2 months). In contrast, observations of mean area fraction for both Larsen C and Amery ice shelves are consistently below the lowest quartile of observations. One possible reason for this over-prediction of the area fraction is that, in reality, meltwater may percolate into the snow and firn pack early in the melt season instead of flowing over the surface into lakes. Another possible explanation is that refreezing may still occur early in the season, thus preventing meltwater from accumulating at the ice shelf surface. As discussed above, any particular ice sheet surface may depart from the parameterization prediction due to random natural variability of self-affine surfaces, however, we consider melt lake metrics over a sufficiently large part of these ice shelves and determine that this is an unlikely explanation for these deviations. Additionally, the dependence of the area fraction parameterization on Hurst exponent is weak enough that it is unlikely that relatively small errors in roughness estimation from ICESat-2 can explain the validation mismatch. We further discuss these and other potential explanations for the mismatch between predicted and observed area fraction of supraglacial melt lakes in the Discussion section.

We also derive another set of predictions by inversely solving for the optimal water supply (w_s) to fit the observed mean area fraction using the mean Hurst exponent and mean standard deviation of topography for both Amery and Larsen C ice shelves. We scale the original runoff estimates from RACMO with a ratio of the aforementioned optimal water supply and mean runoff for each melt season under the assumption that either: some runoff percolates into the snow, some runoff refreezes, or RACMO over predicts runoff. With this altered water supply, we perform the same calculations to develop an additional set of predictions (Fig. 6). Notably, the observations of area fraction primarily reside within the lower 25% of predicted area fractions. The opposite is true for the mean lake depth where predictions are near the lower end of observed lake depth at Larsen C, and are generally less than observed lake depths at Amery ice shelf.

Maximum potential spatially-averaged Antarctic lake size under current roughness parameters
The maximum capacity of a rough surface to retain water is purely a function of roughness statistics20,21. As the water supply increases (S -> 1), the size of supraglacial melt lakes, as captured in equations (10) and (11), eventually reaches a maximum:

w_l^* = 0.6\sigma(1 - 0.41H^{0.6}) (5)

F_l^* = 0.13(3 - 0.13H^{1.3}) (6)

From equations (5) and (6), we see that the standard deviation of roughness (\sigma) directly influences the depth of water filling surface depressions. As \sigma increases, the maximum spatially-averaged lake depth of the surface increases linearly. This is evident in the maximum-potential lake depth estimations mapped in Fig. 7 which largely resemble the map of roughness amplitude (Fig. 2b.), with a spatial average across the continent of 5.67 meters. As expected, ice shelves are predicted to have a lower maximum-potential lake depth, having a spatially averaged depth of 2.6 meters. Grounded ice has a maximum-potential lake depth of 5.9 meters averaged over space.

Lower Hurst exponents produce higher maximum area fraction of supraglacial melt lakes as shown at the Larsen C ice shelf (Fig. 8). This is due to a greater abundance of small-scale depressions in the topography. As the Hurst exponent increases, the maximum-potential area fraction decreases as there are fewer smaller depressions where melt gathers on the surface. As Antarctica has a narrow range of Hurst exponents, there is little variation in the maximum-potential area fraction with a spatial average of 0.38 over the entire ice sheet area, which is the same when averaged over grounded ice or floating ice.

Methods
ICESat-2 land-ice elevation dataset (ATL06)
To quantify the roughness properties of the ice sheet surface, we use measurements from the ICESat-2 satellite altimetry mission. The ATL06 product from ICESat-2 provides surface elevation at high spatial resolution covering nearly the entire ice sheet (less a small hole of missing coverage around the South Pole in the ice sheet interior). ICESat-2 repeats a near-polar orbit, repeating each track over Antarctica every 91 days, with an along-track spatial resolution of 20 meters between elevation measurements. The ATL06 data product measures the land-ice elevation and is a processed version of the raw photon data which measures the travel time between the satellite and the Earth’s surface. Each ATL06 file has elevation estimates for the six beams, three weak and three strong beams. The weak beam estimates are best suited for bright surfaces such as snow and ice and so are used in this study. We select the first available weak beam from each data file in order not to double count data from both tracks. For every ATL06 track analyzed, we filter points in the track that have been registered to have an elevation greater than 1030 meters, which account for the unphysical elevation values.

Each ICESat-2 ATL06 elevation track is divided into sub-tracks (h) 10 km in length to remove the influence of long-range ice sheet shape on the elevation calculation. Since almost all supraglacial melt lakes are less than 10 km in horizontal extent, ice sheet roughness below 10 km is the most relevant for setting the average area and depth of supraglacial melt lakes. Two surface roughness parameters are computed for each sub-track: the standard deviation of elevation (\sigma), and the Hurst exponent (H). The standard deviation of elevation (\sigma) measures the square-root of the variance across the elevation of the subtrack.

Computing the self-affinity of an ICESat-2 altimetry track
The self-affinity of a surface can be quantified through the Hurst exponent, H, which for topographic data measures the relative importance of topographic variations at a range of horizontal length scales1. Put another way, a surface, z(x,y), is self-affine if it can be rescaled by some factors according to the following equation:

s^{-H}z(sx, sy) = z(x,y) (7)

The Hurst exponent has a value between zero and one. Hurst exponents near zero indicate that topography at short horizontal length scales has similar vertical relief to topography at long horizontal length scales. Hurst exponents near one indicate that topography at short horizontal length scales has much less vertical relief than topography at long horizontal length scales (i.e., the topography is dominated by smooth, large-scale features).

We calculate the Hurst exponent using the “power spectral decay” (PSD) method of ref. 24. Other commonly used methods to calculate the Hurst exponent include Detrended Fluctuation Analysis44 and Higuchi’s method45. Here we find the PSD method to be computationally efficient and provide similar results to other methodologies used to compute the Hurst exponent46,47. In our approach, a Fast-Fourier transform (h^(f)) is applied to each sub-track of ICESat-2 ATL06 elevation data (h(x)), where the frequency is the spatial wavenumber along track (f)

h^(f) = \int_0^{10km} h(x)e^{-2\pi ixf}dx (8)

The power spectral density is defined as

P(f) = |h^(f)|

We then fit a power law relationship between the power spectral density (P(f)) and the spatial wavenumber (f) by linear regression of -log(P) against log(f). Fig. 9 The slope of this linear regression is known as the power spectral decay coefficient,

\beta = -logP(f) / logf (10)

This \beta coefficient can be written as a simple algebraic function of the Hurst exponent (H)

H = 2 - (5 - 2\beta) / 2 (11)

Monte-Carlo Supraglacial Melt Drainage Simulations
Our primary goal in this study is to find simple mathematical expressions relating mean melt lake area fraction and depth to surface roughness characteristics and melt water supply. Since we have shown in the prior section that the Antarctic ice sheet surface is self-affine, we focus on such surfaces in developing parameterizations. To achieve this goal, we follow a three-step method: (1) synthetically generate many random self-affine surfaces, (2) simulate water flow and ponding on these surfaces, (3) statistically fit simple mathematical relationships between mean melt lake area characteristics and surface roughness parameters and meltwater runoff supply.

We use a random self-affine surface generator which takes as inputs the Hurst exponent and standard deviation of the height profile. These inputs are used to calculate the power density spectrum and height probability distribution of surface48. As described in ref. 46, the surface is then generated from this height probability distribution through an inverse Fourier Transform method with randomized phases. Each generated surface has boundaries 10 km apart along a flat plain with topographic height with grid spacing of 100 meters.

After generating a surface, we simulate the flow of water into depressions on the surface using a “conditioned-walker" (CW) model developed for this study based on cellular automata49. The CW model takes as inputs: the surface elevation field and the total volume of meltwater supply. In the CW model described originally by ref. 49, a small discrete fraction (a “precipiton”) of the total melt supply is added randomly at a location on the surface. The model moves the precipiton to one of the eight surrounding grid spaces with the lowest elevation and repeats this process until the precipiton reaches a grid space surrounded by grid spaces of higher elevation (i.e., a local minimum) or runs off the surface. This is under the assumption that the surface is non-porous with none of the precipitons penetrating the subsurface. The precipiton water thickness is now added to the surface elevation at this grid space until it equals the elevation of a surrounding grid space nearest in elevation. If water is still left over in the precipiton, then the movement process continues until no water is left. After a single precipiton has been deposited, a new iteration begins with an additional precipiton of water added to a random location on the surface, over a prescribed number of fill stages. This occurs until the prescribed total volume of meltwater has been distributed across the surface. At prescribed intermediate water supply levels, the CW model outputs a field of water depth.

The CW model has some features which make it ideal for considering water flow on relatively flat, rough, ice sheets. First, it allows for the accumulation of water at the capacity of the surface, in contrast to “depression-filling” water routing algorithms. Additionally, the CW model can simulate internally drained catchments that cause ponding on ice sheets, unlike many drainage algorithms used in terrestrial hydrology which assume that all water added to the surface is drained off the surface.

Utilizing the self-affine surface generation algorithm and the CW model, we simulate the distribution of water depth on many surfaces with a given Hurst exponent (H) and roughness amplitude (\sigma) and varied water supply (w_s). For each surface, we calculate the average melt lake depth (w_l) over just the water-covered portions of the surface and the average area fraction of melt lakes (F_l). Since the average melt lake depth (w_l) is just the water depth only over the area fraction of the ice sheet surface covered by water (F_l), we can, a priori, write a relationship between both characteristics and the mean water depth across the entire surface (w_d)

w_l = w_d / F_l (12)

Additional work is done to calculate the mean water depth across the entire surface in the Supplementary Material (Supplementary Fig. 2). Our workflow uses a Monte Carlo approach where all the melt lake metrics described above are averaged over 500 randomly generated surfaces to reduce the influence of single outlier surfaces. The result is a numerically generated dataset of average w_l and F_l values, each corresponding to inputs of w_s, H and \sigma. This multivariate dataset is fit using standard nonlinear least squares fitting techniques to derive simple mathematical relationships between the melt lake characteristics. We set the functional form of equations for each melt lake characteristic to meet certain physical constraints (defined below), while using the simplest relationship between each independent parameter and the resulting melt lake parameters. In general, we prioritize forms of the parameterizations with relatively few parameters, and where no single term can be removed without substantially reducing the R2 fit to the numerically generated dataset (i.e., simplicity). We also enforce two constraints:

w_l(w_s=0, H, \sigma) = F_l(w_s=0, H, \sigma) = 0 (13)

and

d/dw_s (w_l) |w_s->1 = d/dw_s (F_l) |w_s->1 = 0 (14)

The first constraint ensures that when there is no meltwater supply (w_s=0), the melt lake has zero area fraction and depth. The second constraint ensures that all melt lake size statistics asymptotically approach a constant value since finite-scale bumpy surfaces have a finite capacity to store melt water. To enforce the second constraint while obtaining smooth parameterizations, all terms involving the water supply are parameterized in terms of the error function, erf(w_s). We emphasize that there are likely other mathematical forms of these functions that would fit the numerical results equally well, but we search for mathematical functions that are relatively simple and easy to implement in large-scale models.

We generate surfaces with roughness varying over the full range of possible H values ([0,1] with \Delta H=0.1) and a wide range of \sigma values encompassing ~90% of estimates for the Antarctic Ice Sheet (1–20 meters with \Delta \sigma=1 m; see Fig. 2b). For each (H, \sigma) pair, we add meltwater until the surface has reached its maximum water retention capacity (i.e., F_l, w_l and w_d no longer change). For every possible (H, \sigma, w_s) value combination, 500 surfaces are generated and processed under the workflow described in the previous section, making a total of 110,000 simulations to make up the entire parameter space. A mean lake depth and area fraction are calculated for each (H, \sigma, w_s) parameter combination, averaging over each surface domain and all 500 generated surfaces. A schematic of the workflow can be seen in Fig. 10.

Observations from Landsat
We validate our parameterizations with observations of melt lake area and depth from the Landsat-derived product of 35. This product uses Landsat 8 Level 1 imagery that has been corrected to account for Top of the Atmosphere (TOA) reflectance and Sentinel 2 Level 1 C products to analyze the surface conditions of ice sheets and shelves35. Utilizes the Normalized Difference Snow Index (NDSI) to differentiate rocks from snow and the Normalized Difference Water Index to determine the presence of water. Finding the lake depth is achieved by calculating the difference between the reflection of optically deep water and the albedo of the lake bed and the top of atmosphere reflectance of the lake respectively and dividing the difference by the attenuation rate (Equation (15)).

z = ln(A_d - R_1) - ln(R_w - R_1) (15)

The area fraction is calculated as the ratio of the number of pixels classified as surface melt to the total number of land pixels, which are pixels classified either as surface melt or ice. This product has been developed for the entirety of Antarctica though here we focus on Larsen C and Amery ice shelves because: (1) this observational dataset indicates recurrent surface melt on these ice shelves in recent melt seasons, and (2) regional climate models also predict some surface melt on these ice shelves. Many regional climate models (e.g., RACMO v2.3p2 which we use here) predict no surface melt runoff on many other Antarctic ice shelves where melt lakes are observed.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Summary\nThe paper focuses on the development of a parameterization for supraglacial lakes in Antarctica than can potentially be\nused by numerical atmosphere and glacier models to account for the effects of surface meltwater impoundment on the ice\nsheet energy balance and stability, respectively. The authors fit the equation for self-affine surfaces to ICESat-2 ATL06\nelevations for Antarctica in order to empirically parameterize Antarctic surface roughness. They use the\nparameterized roughness coefficients to create synthetic ice sheet surfaces, which are then filled with water to varying levels to generate a\nsufficiently large dataset to develop their supraglacial lake parameterizations. The supraglacial lake parameterizations are\nassessed with respect to optical image-based estimates of supraglacial lake depth and fractional area coverage for the\nLarsen C and Amery ice shelves. When RACMO is used for the meltwater input to the parameterization, lake depths are\nreasonably reproduced but the fractional lake area is over-estimated. When the runoff inputs are tuned to better reproduce\nfractional lake area, lake depths are typically under-estimated.\nOverall this is an interesting topic and the parameterization for supraglacial lakes is an important step forward in\nincorporating complex surface melt processes in numerical models without increasing computational effort/time. I think the\nexplanation of the self-affinity and the Hurst equation are terrific. The writing is mostly clear but there are a few places where\na more detailed description or an additional figure is necessary. Most of my comments are fairly general, but I have a few\nspecific comments for typos and other small details.\nMajor Comments\n1. The use of ICESat-2 elevations to characterize surface roughness is a really cool idea. You don’t need to describe ATL06\nin detail because it is used a lot and there are technical documents that the reader can review if they’d like. However, there\nare several data and methods descriptions in this section that require more detail due to the importance of the ICESat-2 data\ninterpretation on your parameterization.\na. line89: “Depending on availability, all three of the weak beam elevation estimates are used.” Does this statement mean\nthat the ATL06 weak beam data were not always available because of clouds or something else?\nb. line 90: “the entirety of the elevation profile of the track is analyzed to remove any elevation points that are a result of\natmospheric back-scattering, which leads to unrealistic elevation values.” You just mean that you used the ATL06 data\nproduct, which filters-out background signal photons, correct? You didn’t do any additional filtering from what I can tell but\nthis sentence doesn’t clearly ascribe the filtering to the ATL06 data product.\nc. lines 119-123: You state that you filter the tracks that don’t support self-affine surface topography but describe the\ntopography as self affine throughout the paper. Your filtering and interpretation of the ICESat-2 data would really benefit from\na figure and a quantitative description of the filtering. How many (number and fraction) sub-tracks were filtered for Hurst\nexponents <0 or >1? How many were filtered based on the R^2 value of the linear regression to the log-log plot? If you\nadded a figure in the main text that showed all the ICESat-2 tracks as lines with colors corresponding to each one of those\nfiltering criteria, that would be really helpful because the reader would immediately get a sense of whether you have\nremoved important data. Your other figures make it look like you have total coverage of Antarctica, so you could potentially\nget by with simply reporting those metrics on filtered tracks. It’s up to you but more detail should be provided.\n2. For the Hurst equation parameters, is sigma always defined as such: “The standard deviation of elevation (sigma) \nmeasures the average root-mean-square deviation of elevation from the mean value within the sub-track.”? Is this measure\nof variability typically prescribed using the standard deviation in elevations but you use something else here and simply\nkeep the sigma notation to preserve convention? It’s a bit confusing to call something the standard deviation but then not use\nthe standard deviation.\n3. For all your different plots (like figures 4-5), there is a strange oscillation in errors between sigma values. This does not\nseem to be discussed anywhere, unless I missed it, and it is worthwhile to include in the discussion at a minimum. Is this\ndue to a relatively small sample size (500) for such a large parameter space?\n4. The validation using remote sensing estimates of lake extent and depth makes it seem like the modeled lakes are too\nshallow and wide. The remotely-sensed area estimates are more accurate than the depth estimates, correct? That should be\nmentioned in the discussion. A disagreement with remote sensing data doesn’t necessarily mean your results are wrong and\nthat needs to be articulated more clearly. You mention that Fricker at al (2021) compared imagery- and ICESat-2-based lake\ndepths. Were the ICESat-2 lakes deeper than suggested by imagery? You could really strengthen your results with a bit\nmore discussion here.\nLine-by-line Comments\n• Lines 31 & 53: I would not make surface energy balance an acronym. You don’t use it enough to warrant it.\n• Lines 31-43: This paragraph starts to wander a bit. When you get into the latent and sensible energy fluxes, it is not really\nclear what you want the reader to get out of the paragraph. Is the point of the paragraph to explain that supraglacial melt\ninfluences air-ice sheet interactions? I think you could streamline this a bit.\n• Lines 47 & 49: References are needed for these statements.\n• Lines 184-185: The use of “surface area” is a bit tricky here. I think you mean the grid is 10km by 10km with 100m grid\nspacing, but this paper focuses on surface roughness and that roughness inherently increases the surface area relative to a\nperfectly flat surface. I’m not sure how to rephrase it, but think about it.\n• Lines 234-235: These two lines are confusing. You describe each H-sigma pair filled to capacity but then say you create\n500 surfaces for each H-sigma-w combination. Is the filled H-sigma combination just an end-member? Do you create 500\nother surfaces or 500 total? Or do you create 500 surfaces for each H-sigma-w combination that is possible (which would be\nthousands or more simulations)?\n• Lines 314-321: Be careful with the phrasing here. You initially state that you assume refreezing occurs in the winter then\nthat no refreezing occurs in the summer. This isn’t contradictory but it makes it hard to decipher if you evolve the surface\ntopography or leave it static over time.\n• Line 318: Change to “each monthly RACMO output time step” or something similar.\n• Figure 6: Panel B is incorrectly labeled. The caption should also make it clear what circles vs shading mean.\n• Section 5 and elsewhere: For each sigma, you look at the “maximum mean lake depth prediction”. Is this the spatial\naverage for a given H-sigma paired that is filled to capacity? Or are multiple surfaces averaged? If a single surface is used,\nI’d describe it initially as “spatially-averaged maximum lake depth” even though that is longer. The combination of multiple\nambiguous descriptors makes it difficult to read. You also switch the order of “mean” and “maximum” repeatedly throughout.\nPlease check.\n• Lines 385-387: How is this possible if there is no water penetration or drainage? As supply increases, do you just get more\novertopping of disconnected water bodies so that the area expands without much of a change in water depth?\n• Supplement line 14: I think you mean “inert” or “intrinsic” or something along those lines, not “inane”.", "reviewer_2": "Please see attached PDF"}, "author_response": "We sincerely appreciate the constructive feedback provided by the reviewers. We have addressed all\nsuggested edits. Below you can find reviewer comments in blue and our responses in black.\nReviewer 1\nSummary\nThe paper focuses on the development of a parameterization for supraglacial lakes in Antarctica\nthan can potentially be used by numerical atmosphere and glacier models to account for the\neffects of surface meltwater impoundment on the ice sheet energy balance and stability,\nrespectively. The authors fit the equation for self-affine surfaces to ICESat-2 ATL06 elevations for\nAntarctica in order to empirically parameterize Antarctic surface roughness. They use the\nparameterized roughness coefficients to create synthetic ice sheet surfaces, which are then filled with water to varying levels to generate a sufficiently large dataset to develop their supraglacial lake\nparameterizations. The supraglacial lake parameterizations are assessed with respect to optical\nimage-based estimates of supraglacial lake depth and fractional area coverage for the Larsen C\nand Amery ice shelves. When RACMO is used for the meltwater input to the parameterization, lake\ndepths are reasonably reproduced but the fractional lake area is over-estimated. When the runoff\ninputs are tuned to better reproduce fractional lake area, lake depths are typically under-estimated.\nOverall this is an interesting topic and the parameterization for supraglacial lakes is an important\nstep forward in incorporating complex surface melt processes in numerical models without\nincreasing computational effort/time. I think the explanation of the self-affinity and the Hurst\nequation are terrific. The writing is mostly clear but there are a few places where a more detailed\ndescription or an additional figure is necessary. Most of my comments are fairly general, but I have\na few specific comments for typos and other small details.\nMajor Comments\n1. The use of ICESat-2 elevations to characterize surface roughness is a really cool idea. You don’t\nneed to describe ATL06 in detail because it is used a lot and there are technical documents that the\nreader can review if they’d like. However, there are several data and methods descriptions in this\nsection that require more detail due to the importance of the ICESat-2 data interpretation on your\nparameterization.\na. line89: “Depending on availability, all three of the weak beam elevation estimates are used.”\nDoes this statement mean that the ATL06 weak beam data were not always available because of\nclouds or something else?\nAs listed in the ATL06 Version 006 known issues document, there are granules that were available in\nprevious releases are no longer available in later releases such as 006, which we use here. Line 89 was\nappended to the following: “Depending on availability on the 006 version of ATL06, at least one of the\nthree of the weak beam elevation estimates are used.”\nb. line 90: “the entirety of the elevation profile of the track is analyzed to remove any elevation\npoints that are a result of atmospheric back-scattering, which leads to unrealistic elevation values.”\nYou just mean that you used the ATL06 data product, which filters-out background signal photons,\ncorrect? You didn’t do any additional filtering from what I can tell but this sentence doesn’t clearly\nascribe the filtering to the ATL06 data product.\nWe did some filtering of the ATL06 product to ensure that there are no unphysical values in the land\nheight elevation track. For instance, there have been reported land elevation values up to 1030 within an\nICESat-2 data track file. This is non-sensical as this would be five times thicker than Earth’s atmosphere.\nLine 90 was appended to address that we filtered ATL06 elevation tracks.\nc. lines 119-123: You state that you filter the tracks that don’t support self-affine surface topography\nbut describe the topography as self-affine throughout the paper. Your filtering and interpretation of\nthe ICESat-2 data would really benefit from a figure and a quantitative description of the filtering.\nHow many (number and fraction) sub-tracks were filtered for Hurst exponents <0 or >1? How many\nwere filtered based on the R^2 value of the linear regression to the log-log plot? If you added a figure\nin the main text that showed all the ICESat-2 tracks as lines with colors corresponding to each one\nof those filtering criteria, that would be really helpful because the reader would immediately get a\nsense of whether you have removed important data. Your other figures make it look like you have\ntotal coverage of Antarctica, so you could potentially get by with simply reporting those metrics on\nfiltered tracks. It’s up to you but more detail should be provided.\nTwo additional lines were added to the manuscript demonstrating that only 7.35% of the sub-tracks\nanalyzed were tossed out of the total 1.8 million sub-track dataset. Of the 7.3% of sub-tracks that were\ntossed from the data set, 7.1% was due to poor R2 fitting of the linear regression, and the remaining\n0.25% was due to the linear regression producing a Hurst Exponent value that was outside of the\nacceptable range [0,1].\n2. For the Hurst equation parameters, is sigma always defined as such: “The standard deviation of\nelevation (sigma) measures the average root-mean-square deviation of elevation from the mean\nvalue within the sub-track.”? Is this measure of variability typically prescribed using the standard\ndeviation in elevations but you use something else here and simply keep the sigma notation to\npreserve convention? It’s a bit confusing to call something the standard deviation but then not use\nthe standard deviation.\nApologies for the confusion. We use the built-in standard deviation function in MATLAB to compute this\nparameter. Line 97 has been rewritten to clarify the exact mathematical calculation performed to calculate\nsigma.\n3. For all your different plots (like figures 4-5), there is a strange oscillation in errors between sigma\nvalues. This does not seem to be discussed anywhere, unless I missed it, and it is worthwhile to\ninclude in the discussion at a minimum. Is this due to a relatively small sample size (500) for such a\nlarge parameter space?\nThis variation is related to weak nonlinear dependencies on water supply that are not possible to resolve\nin the fit without adding substantial complexity to the chosen mathematical expressions. We have opted to\nkeep our expression simple, with the understanding that such high-order dependencies can be resolved by\nusing the more complex numerical model. The fourth paragraph of the results section has been appended\nto provide a small discussion regarding this matter.\n4. The validation using remote sensing estimates of lake extent and depth makes it seem like the\nmodeled lakes are too shallow and wide. The remotely-sensed area estimates are more accurate\nthan the depth estimates, correct? That should be mentioned in the discussion. A disagreement\nwith remote sensing data doesn’t necessarily mean your results are wrong and that needs to be\narticulated more clearly. You mention that Fricker at al (2021) compared imagery- and ICESat-2-\nbased lake depths. Were the ICESat-2 lakes deeper than suggested by imagery? You could really\nstrengthen your results with a bit more discussion here.\nYes, the area estimates are more accurate than the depth estimates. In Fricker et al (2021) measures the\ndepth of the melt lakes using ICESat-2 altimetry and compares them to the depths that were computed\nusing satellite imagery. It was found that the satellite imagery derived depths underestimate the depth of\nthe melt lakes. Both these points have been added to the discussion:\n“Fricker et al. [2021] compares lake depth measurements of the surface between ICESat-2 and satellite\nimagery (though a comprehensive ice-sheet-wide open-access product has not yet been released using\nICESat data) and concludes that using ICESat-2 improves the accuracy of the melt lake depth...\"\nLine-by-line comments\n• Lines 31 & 53: I would not make surface energy balance an acronym. You don’t use it enough to\nwarrant it.\nEliminated all acronyms in the manuscript.\n• Lines 31-43: This paragraph starts to wander a bit. When you get into the latent and sensible\nenergy fluxes, it is not really clear what you want the reader to get out of the paragraph. Is the point\nof the paragraph to explain that supraglacial melt influences air-ice sheet interactions? I think you\ncould streamline this a bit.\nReworded paragraph to focus on how supraglacial melt lakes effect the surface energy balance of glacier\nsurfaces. Paragraph is now shorter and more streamlined.\n• Lines 47 & 49: References are needed for these statements.\nLutz et~al. 2023 – measures the distribution changes of melt lake changes in Northeast Greenland\n• Lines 184-185: The use of “surface area” is a bit tricky here. I think you mean the grid is 10km by\n10km with 100m grid spacing, but this paper focuses on surface roughness and that roughness\ninherently increases the surface area relative to a perfectly flat surface. I’m not sure how to\nrephrase it, but think about it.\nRewrote the sentence to clarify that the horizontal bounds of the plain surface of 10 km apart with a grid\nspacing of 100 meters.\n• Lines 234-235: These two lines are confusing. You describe each H-sigma pair filled to capacity\nbut then say you create 500 surfaces for each H-sigma-w combination. Is the filled H-sigma\ncombination just an end-member? Do you create 500 other surfaces or 500 total? Or do you create\n500 surfaces for each H-sigma-w combination that is possible (which would be thousands or more\nsimulations)?\nWe have revised these lines to specify the total number of simulations that end up making the parameter\nspace. Each H-sigma pair can be used to randomly generate surfaces. Thus, we sweep over X values of H,\nY values of sigma and, and for each pair of parameter values we randomly generate 500 surfaces. Thus,\nwe run X*Y*500=ZZZ simulations using the cellular automaton model. For each simulation, we consider\noutput at N different levels of water supply.\n• Lines 314-321: Be careful with the phrasing here. You initially state that you assume refreezing\noccurs in the winter then that no refreezing occurs in the summer. This isn’t contradictory but it\nmakes it hard to decipher if you evolve the surface topography or leave it static over time.\nWe have clarified this sentence by changing to “...reset the accumulated supply to zero at the end of each\nmelt season under the assumption that refreezing occurs.”\n• Line 318: Change to “each monthly RACMO output time step” or something similar.\nFixed.\n• Figure 6: Panel B is incorrectly labeled. The caption should also make it clear what circles vs\nshading mean.\nPanel B has been corrected.\n• Section 5 and elsewhere: For each sigma, you look at the “maximum mean lake depth prediction”.\nIs this the spatial average for a given H-sigma paired that is filled to capacity? Or are multiple\nsurfaces averaged? If a single surface is used, I’d describe it initially as “spatially-averaged\nmaximum lake depth” even though that is longer. The combination of multiple ambiguous\ndescriptors makes it difficult to read. You also switch the order of “mean” and “maximum”\nrepeatedly throughout. Please check.\nWe have made changes in this section to refer to “Maximum potential spatially-averaged Antarctic lake\nsize”. This should clarify that we mean the first of your suggestions, the “spatial average for a given H-\nsigma paired that is filled to capacity”.\n• Lines 385-387: How is this possible if there is no water penetration or drainage? As supply increases,\ndo you just get more overtopping of disconnected water bodies so that the area expands without\nmuch of a change in water depth?\nYes, this is correct. Every surface has a maximum water retention (Schrenk et~al., 2014) regardless of\npercolation and maximum saturation of the firn. As water is added to the already saturated surface, excess\nwater will flow off the surface into neighboring cells and surfaces until they are fully saturated as well.\n• Supplement line 14: I think you mean “inert” or “intrinsic” or something along those lines, not\n“inane”.\nFixed.\nReviewer #1 (Remarks on code availability):\nThere is very little documentation with the code. Although there is technically a README, it doesn't\ncontain any information. I did not try to reproduce the code.\nGitHub was updated for all scripts including the README, which now has a description of every script\nin the depository. Updated GitHub can be found in the manuscript as well as here.\nhttps://github.com/dgrau13/meltlake-parameterizations\nReviewer #2\nReview of “Predicting Mean Depth and Area Fraction of Antarctic Supraglacial Melt Lakes\nwith Physics-Based Parameterizations”\nGrau et al. 2024\nSummary\nIn this manuscript, the authors provide parameterizations for supraglacial melt lake\ncharacteristics that could be used in future ice sheet modeling work. Overall, this is an\ninteresting manuscript. It is well-written and the figures are generally clear. I have some\nmajor concerns relating to the methodology. Once these concerns are addressed, I believe\nthis work will be a good contribution as parameterizations for supraglacial meltwater will\nbe increasingly important in ice sheet modeling.\nOne general comment: The structure of this manuscript is a bit odd. Sections 2-4 each\nhave individual methods and results subsections. I don’t necessarily mind this structure\nand I found it relatively intuitive and easy to follow, but I’m not sure that it aligns with the\njournal requirements. If the current structure is maintained, I would recommend using\nmore verbose titles for sections 2.2, 2.3, 3.1, 3.2 and 4.3 so there are not multiple\nsubsections titled simply “methods” or “results”. For example, Section 2.2 could be\nrenamed “AIS self-affinity methodology” or something of the sort.\nMajor Comments\nGenerally, I find some methodological details lacking. I also believe that an overarching\ndiagram or schematic of all the methods would be helpful.\nAdditional figures depicting the self-affinity of an icy surface and the methodological workflow described\nin Section 3 have been added.\nSelf-affinity (Section 2)\n• what quality proves that a surface is self-affine? It is unclear to me from this work, what\nactually demonstrates that the AIS is self-affine. Is it the Hurst exponent that proves\nself-affinity?\nA surface is self-affine if the scale of vertical relief with respect to horizontal scales can be described by a\npower law distribution. Thus, the fact that we can statistically fit more than 90% of the surface area of\nAntarctica to such a power law distribution indicates that the AIS surface is largely self-affine. This fact is\nnow stated plainly with reference to these results in section 2.2.\n• The paragraph from L119-123 is confusing. In this section, I thought you were trying to\ndemonstrate that the AIS is self-affine. But in these lines, the manuscript mentions that\ntracks are discarded if they are not described as self-affine or that “the only sub-tracks\nconsidered are those in which the data strongly indicate roughness is self-affine”. How\ncan you demonstrate self-affinity if you self-select only tracks that are self-affine?\nAdditionally, I find some details to be missing. For example: What percentage of the\nsub-tracks are selected? How many have R2 < 0.7?\nWe have modified this section to be much clearer:\n“ Of the 1.8 million sub-tracks that were analyzed, 7.3% of the sub-tracks were discarded during the\nanalysis process. As a result, over 92% of all the subtracks analyzed exhibit the qualities of a self-affine\nsurface. Consequentially, we affirm that the Antarctic Ice Sheet surface is largely self-affine.”\nMean lake depth and area fraction (Section 3)\n• L206 – “with varying water supply depth, Hurst exponent, and roughness amplitude”.\nWhat do these values vary between? I guess this is in L232-233 but this information\nshould come in the methods section. Also, more details are necessary. Are they\nrandomly selected from the distributions found in Section 2 (i.e. Fig 2,3)? Can you elaborate on the\n“water supply depth”? Is this the same as the “meltwater availability”?\nWhere does this information come from?\nWe have modified this for clarity:\n\"Utilizing the self-affine surface generation algorithm and the CW model, we simulate the distribution of\nwater depth on many surfaces with a set range of Hurst exponent (H) and roughness amplitude (𝜎) values\nand varied water supply (w).” Line 206 has been rewritten to specify that the water supply varies from a\ns\nsurface-to-surface basis. It also addresses that the range of the Hurst exponent and roughness amplitude is\ndefined from a set range. This paragraph was also moved from the results section to the methods section\nexpanding on the set-up of the model including the increment steps of the Hurst exponent and roughness\namplitude of the parameter space. This water supply depth is merely the volume of the meltwater supply\ndivided by the planar area of the generated surface. This water supply depth is equivalent to the meltwater\nsupply or “runoff” as calculated in SMB models. Further clarification has been added.\n• L232 – is it realistic to do this over the range of possible H values [0,1] since it seems\nlike >99% of the calculated H values fall between 0.3 and 0.6 (Fig. 2).\nThere are still some outlier values on the Antarctic Ice Sheet that fall outside of this range. We wanted to\nbuild a parametrization that could be applied to a variation of surface roughness conditions not one\nexclusive to the Antarctic Ice Sheet.\nValidation of lake parameterizations (Section 4)\n• Can you provide the equation used to calculate depth? Which bands are used? What\nare the limitations and uncertainties associated with these depth calculations?\nThroughout this section it seems that you treat the ‘observed depths’ as the truth.\nHowever, there are large uncertainties when using optical imagery to determine lake\ndepth. A greater discussion of this is necessary.\nRed and Blue bands are used to compute the Normalized Difference Water Index, which is able to detect\nwhether a pixel contains water. We have added some general details about the calculation made in\nMoussavi et al. (2020), but the complete details for these calculations are in that original paper and\ntherefore we opted not to repeat them in their entirety here. A further discussion of the uncertainties\ninherent in these calculations has been added.\n• It appears from Figure 6 that for the lake depth and area fraction predictions, it is\nassumed that meltwater refreezes before July? In reality, melt refreezes much before\nthen. Can you better temporally constrain the predictions? Maybe using RACMO air\ntemperature?\nThis validation purely depends on the runoff estimates from RACMO. It would be challenging to attempt\nto constrain the predictions using the air temperature estimates, especially when attempting to scale the\namount of depth that would be refrozen due to air temperature. An additional model or parameterization\nwould be needed to implement this, unless we assume once the air temperature reaches a stable value\nbelow freezing that all the melt lakes freeze. Therefore, we have opted here to take RACMO outputs at\nface value as would be available in large-scale models with SMB simulated.\n• Is there a way to include firn air content in this parameterization, thereby accounting for\npercolation?\nThere is a way to account for percolation in this parameterization by reducing the input melt water supply\nby assuming that “x” meters of melt water supply percolate into the firn leaving “x” meters of melt water\nsupply to pool on the surface. This percolation could be considered to be depedent on firn air content.\nHowever, this is beyond the scope of our study, and indeed there are prior studies considering such\npercolation.\n• L344 – what is the scaling factor used?\nModified to:\n“We scale the original runoff estimates from RACMO with a ratio of the aforementioned optimal water\nsupply and mean runoff for each melt season under the assumption that either: some runoff percolates\ninto the snow, some runoff refreezes, or RACMO overpredicts runoff.” The optimal water supply is\nderived by solving inversely for w from the observational area fraction.\ns\nDiscussion (Section 6)\n• L369 –“The parameterizations developed in this study are mostly consistent with\nobservations. ” I’m not really convinced by this. I am especially concerned with Figure\n6b/d. Considering there are huge uncertainties in both the observations (lake depth\nfrom optical imagery) and the parameterizations (not accounting for firn percolation),\nI’m not convinced by the comparison to observations.\nWe have modified this to account for some of the shortcomings identified by the reviewer:\n“Of the parameterizations developed in this study, the mean melt lake parameterization is mostly\nconsistent with observations.” Ultimately, the two points brought up by the reviewer above are\nacknolwedged and dicussed in the text. We feel that the manuscript is stronger with a validation section\nthan without, and absent more concrete suggestions about how we can improve the validation, we feel\nthat it adds to the arguments made in this study.\n• L390 – which areas specifically are most at risk of future lake expansion? Could there\nbe some more advanced future insight included in this manuscript, besides just that\nother areas may be at risk of melt pond development?\nAnswering the question posed here is ultimately best left up to a more detailed regional climate model\nwhich can predict where runoff may occur. We do not purport to answer this question in our study, which\nhas already been tackled by many other studies. Instead, we answer a question that has received much less\nattention in the literature, namely: “If melt were to occur, what would be the size of melt ponds?” This\nquestion gets at different issues related to albedo and potential hydrofracture that previous studies focused\non melt production are not able to answer.\n• L392-393-“This may explain why, … observed lake depths span a relatively narrow\nrange” . Or perhaps, the algorithm used to calculate lake depth has saturated. Previous\nwork (e.g. Lutz et al 2024) has shown that lake depths computed using the red band\nsaturate at low depth, therefore a full range of lake depths cannot be calculated using\nthis method. This limitation should really be discussed in comparison to observations.\nWe have added a discussion and reference to Lutz et al. (2024). Thank you for bringing this new study to\nour attention.\nAbstract – I think you need a concluding sentence in the abstract that summarizes the\nmotivation and importance of this work. This should tie back to the first 2 sentences of the\nabstract. As written now, the significance of this work is not immediately clear from the\nAbstract. \nAn additional sentence was added to the end of the abstract that restates the overall applications and limits\nof the parameterizations.\nMinor Comments\nL21 – I don’t think that the Kingslake et al. (2017) paper showed an increase in surface\nmelting and supraglacial lake formation.\nRemoved increase so sentence now reads: “... observations have shown widespread surface melting and\nsupraglacial lake formation on Antarctic ice shelves...”\nL31-43 – This paragraph goes into a lot of detail about the impact of supraglacial meltwater\non surface energy balance and energy fluxes. This seems like too much detail considering\nthe manuscript is not about ice sheet surface energy balance. I would appreciate more\nbackground introduction on topography theory instead of energy fluxes.\nThe introduction is amended to have a small section discussing the attributes and properties of self-\naffinity and less detail on surface energy balance processes.\nL48 –“Supraglacial melt lakes are omitted from most current large-scale models…” Are there any\nmodels that incorporate them?\nThis is correct, we have deleted the word “most”.\nL53 – The last sentence of this paragraph should be reworded to be more specific and\nstronger regarding the contribution of this work.\nWe expanded the last sentence in this paragraph to be more direct with the expectations and applications\nto come from this work.\nL55-59 – It would be helpful to have more background information here\nThe introduction is amended to have a small section discussing the attributes and properties of self-\naffinity and surface evolution.\nL59-69 – Relate the work done back to the original knowledge gaps.\nAn overarching sentence was added describing the applications of this work and the intent.\nL111 –“… provides similar results to other methods. ” Did you compare to other methods, or\nwas this found in a previous study?\nThis was found in previous studies comparing different methods. Two citations were added at the end of\nthe sentence for readers interested in learning more about the different methods and theory behind\ncalculating the Hurst exponent.\nL130 –“There is a clear difference in Hurst exponent between grounded and floating ice” . Is\nthis difference statistically significant?\nIt is statistically significant. The p-value and t-test between both distributions are respectively zero and 1.\nWe have added a statement to this effect.\nL149 –“It was found that the Hurst values do not vary much…” How much? Can this be\nquantified in the main text?\nWe have added a statement indicating that the mean difference of the value of the Hurst exponents at the\nintersection point of the ICESat-2 track is 0.0218.\nFigure 2,3 – Why do you show a PDF for the Hurst exponent and a CDF for the standard\ndeviation?\nTo show that most of the values for the standard deviation of topography falls within 25 meters. There is a\nfigure in the supplemental that demonstrates the PDF for the Hurst of the grounded ice and floating ice\nrespectively, which exhibits a distribution like that of a Poisson distribution.\nL167 – can you provide a brief description/explanation of “cellular automata”?\nWe have appended sentence to include a brief definition of cellular automata.\nL174 –“meltwater runoff depth” This is the same as “meltwater availability” , right?\nWe have changed the phrase “meltwater runoff depth” to “meltwater supply”.\nSection 3.1 – Somewhere in this section it is important to mention that these methods\nassume a non-porous surface and that there is no infiltration and all precipitons pond on\nthe surface.\nWe added the following sentence:\n“We also assumed that the surface is non-porous with none of the precipiton penetrating the subsurface.”\nL216 – Can you elaborate what is meant by “custom equation forms”?\nWe have modified the sentence to read:\n“We set the functional form of equations for each melt lake characteristic to meet certain physical\nconstraints (defined below), while using the simplest relationship between each independent parameter\nand the resulting melt lake parameters.”\nL234 – when writing the “H – 𝜎 pair” (and later for H – 𝜎 - 𝑤!), is there a better way to write\nthis that does not look like subtraction? Maybe (H, 𝜎)\nNotation changed to (H, 𝜎) and (H, 𝜎,w).\ns\nL245 –“… we create a simple mathematical expression…” How is this created? Can more detail be\nprovided?\nThis part modified to read:\n“ We calculate the maximum depression storage in a separate set of simulations using the full Hurst and\nstandard deviation of topography range by over-saturating the generated surface and extracting the mean\ndepth across the entire surface. We then fit a simple mathematical expression that can determine the\nmaximum average water depth of any self-affine surface utilizing only surface roughness parameters:”\nFigure 4/5 – Please ensure that all color bars are labeled.\nWe have added labels to all color bars.\nAlso, in 4a, for large 𝜎, why is the average lake depth higher for lower total water supply?\nThe deepest depressions also tend to be the most laterally extensive (the result of self-affinity of the\nsurface) and so these fill in first due to capturing more randomly generated surface melt, producing the\ndeepest lakes. As further water supply is added small depressions are filled in, thus lowering the average\nlake depth. This is a relatively minor second-order effect and so we do not prioritize this particular feature\nin our simple parameterizations which are mainly meant to capture the initial rapid increase in lake depth\nat low water supply and the dependence of lake depth on sigma.\nIn Figure 5, the relative difference is almost always negative. This may indicate that F is not\nparameterized properly?\nThe relative difference is almost always positive in Figure 5. We have prioritized meeting the physical\nconstraints of the parameterizations which produces some minor systematic error in this case.\nL309 – Specify the “predictions” you aim to provide (e.g. mean lake depth, lake area\nfraction).\nModified to:\n“We aim to validate our newly developed parameterizations by computing a range of mean melt lake\ndepth and area fraction predictions for the Amery and Larsen C ice shelves...”\nL314 – What is meant by “the regions within the two ice shelves where melt lakes are\nObserved”?\nModified to:\n“Due to the difference in time steps and spatial grid of the observational products and predictions, we\ninstead include all ICESat-2 derived topographical parameters and RACMO seasonal accumulated runoff\nover the regions within the two ice shelves where melt lakes were observed in Moussavi et al., 2020.”\nFigure 6 – what do the different shades of blue represent? Why is it often cut-off at the top\nof the plot? The y-axis should be extended to show the full range of predicted values. Also,\nwould it be possible to include RACMO meltwater production on these plots?\nFigure 6 has been redone to include a legend describing each component of the plot, and an additional\naxis added to incorporate the runoff (melt supply) for the validation.\nL361 – When thinking about maximum lake depth, it does not make sense to me that\n“smaller more frequent bumps” would “correlate with a higher mean maximum water\ndepth” . If the bumps are smaller, then at maximum capacity, wouldn’t the mean lake depth\nalso be smaller? Maybe I misunderstand the interpretation of ‘mean maximum lake depth’.\nIs this the mean depth of all water pixels? Or the average of the maximum lake depth for all\nlakes? Or maybe a change to the word “smaller” would be helpful.\nWe mean the mean depth of all water pixels. The previous changes to indicate that this is a spatial average\nclarify this point. Additionally, we have added the following modification:\n“Lower Hurst exponents indicate smaller more frequent bumps along a rough surface and correlate with\nhigher mean maximum area fraction values as shown at the Larsen C ice shelf (Figure 9).”\nL380 – what, specifically, is meant by “melt lake characteristics”?\nMelt Lake Characteristics would be referring to the melt lake depth. Lined has been changed to clarify.\nSpecific Comments\nL12 – Change “ice sheet surfaces” to “the Antarctic Ice Sheet surface” since you only focus\non Antarctica in this work.\nIn this study, we first use extensive surface elevation measurements from the ICESat-2 satellite altimetry\nmission to show that roughness on the Antarctic Ice Sheet surface is largely self-affine, consistent with\nprior observations of bed roughness beneath ice sheets and geomorphic surfaces more broadly.\nL33 –“As melt causes ice sheet surface warming…” This sentence is a bit strangely worded.\nConsider changing to something like: “Ice sheet surface melting initiates a positive\nfeedback loop, whereby a lower surface albedo allows for more shortwave absorption,\nfurther enhancing meltwater production and deepening lakes. ” Or something similar.\nModified to:\n“Ice sheet surface melting initiates a positive feedback, whereby a lower surface albedo allows for more\nshortwave absorption, further enhancing meltwater production and deepening lakes.”\nL34 – add “surface” before “albedo”\nModified to:\n“Ice sheet surface melting initiates a positive feedback, whereby a lower surface albedo allows for more\nshortwave absorption, further enhancing meltwater production and deepening lakes.”\nL35 – change “causes” to “promote”\nModified to:\n“After refreezing, albedo is permanently altered from its original value due to snowpack and firn density\nchanges from meltwater saturation, which promotes future melt lake formation at these locations...”\nL47 – change “fracture” to “fracturing”\nModified to:\n“Incorporating supraglacial melt lakes into large-scale models is thus vital in properly simulating the\npotential interaction between surface melt and ice sheet fracturing and calving.”\nL84-87 –“The ATLO6 data product…” For readability, reword this sentence to: “The ATLO6\ndata product measures the land-ice elevation and is a processed version of the raw photon\ndata which measures the travel time between the satellite and the Earth’s surface\n(available as the ATLO3 product, thou not used in this study)”.\nModified to:\n“The ATLO6 data product measures the land-ice elevation and is a processed version of the raw photon\ndata which measures the travel time between the satellite and the Earth’s surface.”\nL128 – For “little variation” (twice in this line), please specify temporal or spatial.\nModified to:\n“... we found little temporal variation in surface roughness parameters over that period. Figure 2 also\nshows little spatial variation of the Hurst exponent across the ice sheet.”\nL139 – could “bumps” be changed to something like “topographic variability”?\nModified to:\n“This indicates the possibility that processes related to snow redistribution and densification may\ninfluence surface topography, which creates smaller-scale topographical variability [Lowe et al., 2007].”\nL157 – Add “temporally” to “… assuming the ice sheet surface does not change”.\nModified to:\n“Since water gathers in surface depressions, the statistics of roughness on the surface can be used to\ndetermine the statistics of lake size on that surface, assuming that the surface topography does not change\ntemporally...”\nL208 – Move (𝑤\" ) after “average melt lake depth” (L207) for clarity.\nModified to:\n“For each surface, we calculate the average melt lake depth (w)...”\nl\nL211 – Change “for” to “to calculate”\nModified to:\n“Additional work is done to calculate the mean water depth...”\nL323 – specify temporal or spatial \"variations”\nModified to:\n“Additionally, we assume that topographical parameters measured by the ICESat-2 mission in 2021-2022\napply to predictions over the 2014-2018 period when observations are available due to the relatively small\ntemporal variations in these parameters throughout the ICESat-2 mission.”\nL403 – Change “Due” to “Given”\nModified to:\n“Given to the relatively minor variation in Hurst exponents...”\nL414 – Specify “Greenland subglacial topography”. Also, why might the Hurst exponent for\nsubglacial topography differ from that for terrestrial geomorphic surfaces?\nFixed. This is a complicated question and one that we are hoping to address in ongoing work. It is likely\ndue to the fact that subglacial topography is subject to qualitatively different erosive processes than\nsurface ice sheet topography. However this is beyond the scope of this manuscript.\nReferences\nLutz et al., (2024). Assessing supraglacial lake depth using ICESat-2, Sentinel-2, TanDEM-\nX, and in situ sonar measurements over Northeast and Southwest Greenland, The\nCryosphere. https://doi.org/10.5194/tc-18-5431-2024."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Review of “Predicting Mean Depth and Area Fraction of Antarctic Supraglacial Melt Lakes\nby Grau et al.\nsubmitted for consideration in Nature Communications\nSummary\nThe authors implemented relatively modest revisions in order to address comments made by two reviewers during the\nprevious round of review. While many of the revisions improved the manuscript, I found that a few of them were not\nbeneficial. Below, I call attention to the detrimental revisions and offered some recommendations for how they can be\nimproved. A handful of other revisions are also included below.\nMajor Comments\n1. You do not explicitly define the terms in Eqn. 2. For example, presumably h is the ICESat-2 elevation but that needs to be\nstated. Additionally, the bounds of the integration are listed as infinite but your sub-tracks are 10 km in length and it seems\nas though your frequencies are restricted to 10^-3.5-10^-9.2 (based on Fig. 2) but that I not defined or explained. Please\nprovide more detail in the formulation of this equation.\n2. In my opinion, Figure 1 does not help convey self-affinity. I always imagine snowflakes when thinking about fractals (see\nhttps://personal.math.ubc.ca/~cass/courses/m308-03b/projects-03b/skinner/ex-dimension-koch_snowflake.htm). I like the\nmulti-scale approach, with the zoom insets, but each zoom should show the same level of complexity but at different scales. I\nimagine something in which large chevron-shaped lakes are visible at the largest scale and you zoom in once or twice to\nshow that there are also smaller chevron-shaped lakes that are not as visible at the larger scales. That would be much more\neffective.\n3. For Figures 8 and 9, the color scheme does not appear correct based on the legend. Based on my interpretation of these\nfigures, the darker blue color brackets the 25-50th percentiles and the lighter blue brackets the 5-95th percentiles. The\nlegends suggest the opposite.\n4. I appreciate the inclusion of two different predictions for the validation regions but I think that the discussion of why either\nthe area or depth parameterization is in error is still lacking. One of the major discussion points that is missing is a\ndiscussion of the elevation data sub-track size and the frequency range used for solving the Hurst values. In Figure 2, you\ncan see that the slope of the line is not the same across the full range of frequencies under consideration. If you focused on\nfrequencies from ~10^-5 to ~10^-3 only, you would have a much shallower slope for that best-fit line. You mention the\ninfluence of the Hurst exponent on results in lines 382-388 but that discussion tries to argue that Antarctica has little variation\nin the exponent and does not mention how your interpretation of the data could potentially influence the agreement with\nobservations.\nLine-by-line Comments\n• Lines 20-21: The new sentence added at the end of the abstract is somewhat awkward and does not accomplish the goal\nof articulating the impact of the work. Please revise further to circle back to the first two sentences that focus on albedo and\nice shelf stability. Your penultimate sentence describes potential limits on pond size and depth. What are the potential\nimplications related to albedo and hydrofracture?\n• Lines 91-92: This sentence still reads as though you only select certain weak beams because you cannot expect the\nreader to know that sometimes data is not available from all three beams. Please rephrase to indicate you use all available\nweak beams.\n• Lines 122-129: Why not include the same level of detail about filtering as you described in the response to the reviewers?\n• Line 203: You do not allow percolation into the subsurface so I would not use the term “percolation” to describe movement\nof water across the surface.\n• Line 231-232: I agree that you need to explicitly state that you do not permit percolation into the subsurface but this\nstatement should come in the previous paragraph.\n• Line 306 and onwards: Why is the “s” in Landsat capitalized here when it is typically not in literature?\n• Lines 307-318: I like the additional detail that was added but this paragraph could be condensed. Sentences 2-4\nessentially repeat the same points but with slightly different details.\n• Lines 336-337: Please rephrase. The revisions to this sentence did not improve readability.\n• Lines 416-417: I would start this paragraph with “Our melt lake parameterization suggests…” or something along those\nlines. Currently it is unclear if you are making a statement about albedo based on your work or you should be referencing\nliterature.\n• Line 422: Broken reference.\n• Line 477: That Github repository link is broken.\n(Remarks on code availability)\nThe link took me to a 404 error page.", "reviewer_2": "I appreciate the work that the authors put into the second version of this manuscript. The methodology is much clearer and I appreciate the addition of Figure 5. I agree with all of the changes made and only have a few minor additional comments/remarks:\nFigure 1 – I’m not really sure what I am supposed to get out of the Figure. I think some labels and a more descriptive caption would be helpful. What are the blue spots? Also, I don’t think this Figure is referenced in the text (please check this for other figures as well)\nIn Figure 5 – does the height of the cells in the cellular automaton model indicate their relative elevation? It might be helpful to label the ‘z’ axis in this figure.\nL272 – The sentence beginning with “As the depth of the melt lake…” does not make sense\nL317 – Include space after equation\nFigures, in general – Many of the figures have random horizontal and vertical lines around parts of the subfigures (eg. Figure 3, 11). This could be just in the rendering of the PDF but perhaps check to make the figures cleaner.\n(Remarks on code availability)"}, "author_response": "We sincerely appreciate the constructive feedback provided by the reviewers. We have\naddressed all suggested edits. Below you can find reviewer comments in blue and our\nresponses in black.\nReviewer #1 (Remarks to the Author):\nReview of “Predicting Mean Depth and Area Fraction of Antarctic Supraglacial Melt Lakes\nby Grau et al.\nsubmitted for consideration in Nature Communications\nSummary\nThe authors implemented relatively modest revisions in order to address comments made\nby two reviewers during the previous round of review. While many of the revisions improved\nthe manuscript, I found that a few of them were not beneficial. Below, I call attention to the\ndetrimental revisions and offered some recommendations for how they can be improved. A\nhandful of other revisions are also included below.\nMajor Comments\n1. You do not explicitly define the terms in Eqn. 2. For example, presumably h is the ICESat-\n2 elevation but that needs to be stated. Additionally, the bounds of the integration are listed\nas infinite but your sub-tracks are 10 km in length and it seems as though your frequencies\nare restricted to 10^-3.5-10^-9.2 (based on Fig. 2) but that I not defined or explained. Please\nprovide more detail in the formulation of this equation.\nThe bounds of integration have been changed to [0 to 10 km], which covers the full range of length\nscales in a Icesat-2 subtrack. Additionally, we have defined variables throughout the text prior to\nthe derivation of the Hurst exponent. Figure 2 has also been updated to display the power spectral\ndensity in a loglog plot with respect to the wavenumber. The slight curving towards the higher\nfrequencies in this loglog plot has also been explained further in the figure caption:\n“The slight deviation from the linear regression at higher wave numbers is due to the ATL06\nproduct processing aggregate photon counts in 40 meters intervals (Brunt et al. 2019), thus\nartificially smoothing roughness at the highest resolutions of the product”\n2. In my opinion, Figure 1 does not help convey self-affinity. I always imagine snowflakes\nwhen thinking about fractals (see https://personal.math.ubc.ca/~cass/courses/m308-\n03b/projects-03b/skinner/ex-dimension-koch_snowflake.htm). I like the multi-scale\napproach, with the zoom insets, but each zoom should show the same level of complexity\nbut at different scales. I imagine something in which large chevron-shaped lakes are visible\nat the largest scale and you zoom in once or twice to show that there are also smaller\nchevron-shaped lakes that are not as visible at the larger scales. That would be much more\neffective.\nFigure 1 has been updated substantially to emulate the self-repeating pattern at different horizontal\nlength scales referenced by the reviewer. We have drawn inspiration from prior illustration of self-\naffine geomorphic roughness in updating this figure.\n3. For Figures 8 and 9, the color scheme does not appear correct based on the legend. Based\non my interpretation of these figures, the darker blue color brackets the 25-50th percentiles\nand the lighter blue brackets the 5-95th percentiles. The legends suggest the opposite.\nFigure legends have been corrected.\n4. I appreciate the inclusion of two different predictions for the validation regions but I think\nthat the discussion of why either the area or depth parameterization is in error is still lacking.\nOne of the major discussion points that is missing is a discussion of the elevation data sub-\ntrack size and the frequency range used for solving the Hurst values. In Figure 2, you can see\nthat the slope of the line is not the same across the full range of frequencies under\nconsideration. If you focused on frequencies from ~10^-5 to ~10^-3 only, you would have a\nmuch shallower slope for that best-fit line. You mention the influence of the Hurst exponent\non results in lines 382-388 but that discussion tries to argue that Antarctica has little\nvariation in the exponent and does not mention how your interpretation of the data could\npotentially influence the agreement with observations.\nWe have done additional tests on the sensitivity of the predictions to the Hurst exponent and have\nfound that changing the Hurst exponent +/- 0.2 from the mean H = 0.4 doesn’t not alter the\npredictions significantly due to the weak dependence of area fraction on Hurst exponent. We have\nadded a sentence indicating this lack of sensitivity in Section 4.3. We note that in the example\ngiven, the main departure from a Hurst exponent of 0.4 is at very small horizontal length scales.\nSince most supraglacial melt ponds are hundreds to thousands of meters in horizontal dimension,\nsmall inaccuracies in the roughness at the smallest length scales do not play a significant role in\nsetting melt lake size. We have also added a note in the caption of Figure 2 explaining that this\ndeparture at the highest frequencies of the PSD is a result of how the ATL06 product is processed.\nLine-by-line Comments\n• Lines 20-21: The new sentence added at the end of the abstract is somewhat awkward and\ndoes not accomplish the goal of articulating the impact of the work. Please revise further to\ncircle back to the first two sentences that focus on albedo and ice shelf stability. Your\npenultimate sentence describes potential limits on pond size and depth. What are the\npotential implications related to albedo and hydrofracture?\nWe have modified the abstract as follows:\n“In this paper, we present a simple method for parameterizing supraglacial melt lakes in large\nclimate and ice-sheets model that has an intrinsic limit in potential future melt lake size. The\nimplementation of these parameterizations will improve the simulation of albedo and\nhydrofracture of water-filled crevasses in large-scale climate and ice sheet models. Ultimately,\nsuch parameterizations can be used to improve predictions of future ice sheet and ice shelf collapse\nin response to climate change.”\n• Lines 91-92: This sentence still reads as though you only select certain weak beams\nbecause you cannot expect the reader to know that sometimes data is not available from all\nthree beams. Please rephrase to indicate you use all available weak beams.\nLine has been changed to:\n“We select the first available weak beam from each data file in order not to double count data from\nboth tracks.”\n• Lines 122-129: Why not include the same level of detail about filtering as you described in\nthe response to the reviewers?\nParagraph has been altered to:\n“All sub-tracks with calculated Hurst exponents less than zero or greater than one are discarded,\nsince they cannot be reliably described as self-affine (i.e. outside the acceptable confines of the\nrange for the Hurst exponent). Typically, such cases are related to spurious elevation data causing\nlarge deviations in the power spectral density. Additionally, any linear regressions with R2 values\nless than 0.7 are discarded to ensure that the only sub-tracks considered are those in which the data\nstrongly indicate roughness is self-affine. Of the 1.8 million sub-tracks that were analyzed, 7.3%\nof the sub-tracks were discarded during the analysis process. 7.1% were discarded due to poor R2\nvalues when fitting, and the remaining 0.2% were discarded due to the calculated Hurst exponents\nbeing outside the acceptable range of [0,1]. As a result, over 92% of all the subtracks analyzed\nexhibit the qualities of a self-affine surface. Consequentially, we conclude that the Antarctic Ice\nSheet surface is largely self-affine.”\n• Line 203: You do not allow percolation into the subsurface so I would not use the term\n“percolation” to describe movement of water across the surface.\n“Percolation” has been changed to “Movement”.\n• Line 231-232: I agree that you need to explicitly state that you do not permit percolation\ninto the subsurface but this statement should come in the previous paragraph.\nWe are not certain whether this comment applies to these lines, but we have moved the statement\nto indicate that the conditioned walker model does not simulate subsurface percolation earlier in\nthe description of the model. Line 209-210 now reads:\n“This is under the assumption that the surface is non-porous with none of the precipitons\npenetrating the subsurface.”\n• Line 306 and onwards: Why is the “s” in Landsat capitalized here when it is typically not in\nliterature?\nAll of the “s”s in Landsat have been uncapitalized.\n• Lines 307-318: I like the additional detail that was added but this paragraph could be\ncondensed. Sentences 2-4 essentially repeat the same points but with slightly different\ndetails.\nThe paragraph has been reduced to four sentences. Some details have been kept due to other\nreviewer’s concern for understanding how the validation data set it is computing melt lake\ndepths.\n• Lines 336-337: Please rephrase. The revisions to this sentence did not improve readability.\nLine had been rephrased to the following:\n“This results in a distribution of monthly predictions calculated from equations 11 and 12, with\ninputs taken over the area where supraglacial lakes are observed in Landsat, and including: surface\nroughness parameters estimated as described in Section 2, and surface melt supply from RACMO.”\n• Lines 416-417: I would start this paragraph with “Our melt lake parameterization\nsuggests…” or something along those lines. Currently it is unclear if you are making a\nstatement about albedo based on your work or you should be referencing literature.\nLine has been altered to the following:\n“Our melt lake parameterization suggests that the albedo effect of melt lakes is not any greater on\nrougher ice surfaces than on smoother ones, since the area fraction of melt lakes does not depend\non the surface roughness amplitude.”\n• Line 422: Broken reference.\nCitation has been fixed in Line 422.\n• Line 477: That GitHub repository link is broken.\nReviewer #1 (Remarks on code availability):\nThe link took me to a 404-error page.\nLink has been adjusted and should take user to the corresponding repository.\nReviewer #2 (Remarks to the Author):\nI appreciate the work that the authors put into the second version of this manuscript. The\nmethodology is much clearer and I appreciate the addition of Figure 5. I agree with all of the\nchanges made and only have a few minor additional comments/remarks:\nFigure 1 – I’m not really sure what I am supposed to get out of the Figure. I think some labels\nand a more descriptive caption would be helpful. What are the blue spots? Also, I don’t think\nthis Figure is referenced in the text (please check this for other figures as well)\nWe have updated Figure 1 substantially to demonstrate an idealized self-affine glacier inspired by\nFigure 1 in Cael et al., 2017. The blue spots are supraglacial melt lakes on the glacier. We have\nadded more detailed explanation in the caption.\nIn Figure 5 – does the height of the cells in the cellular automaton model indicate their\nrelative elevation? It might be helpful to label the ‘z’ axis in this figure.\nThe height of the of the cells in the cellular automaton model indicates their elevation. Axes labels\nhave been added to this figure.\nL272 – \"The sentence beginning with “As the depth of the melt lake…” does not make sense\nSentence has been rewritten to:\n“On a self-affine surface, the average horizontal extent of lakes filling depressions does not depend\non the average roughness height (Cael et al. 2017). Therefore, we expect that the area fraction of\nsupraglacial melt lakes should have no dependence on σ.”\nL317 – Include space after equation\nSpace was added after equation.\nFigures, in general – Many of the figures have random horizontal and vertical lines around\nparts of the subfigures (eg. Figure 3, 11). This could be just in the rendering of the PDF but\nperhaps check to make the figures cleaner.\nWe have determined this is a rendering issue in Adobe Acrobat. We have remade the figures to\neliminate the issue."}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors appear to have addressed all comments from the last round of revisions.", "reviewer_2": "I appreciate the changes made to Figure 1 in the manuscript. Overall, I am happy with the modifications that the authors\nhave made and am happy to accept the article for publication."}, "author_response": ""}]

Decision letter