We thank the reviewers for their thoughtful and insightful comments, which helped improve our manuscript. We outline the revisions we plan on making to the manuscript in response to the reviewer feedback and address their specific comments and evaluations below. Broadly, we will implement additional clarification in the form of additional text and figures regarding the choices made in our methodology in the random forest, training labels, and implications of constraints in the methods and available datasets. We will further contextualize our results in relation to findings from previous studies by performing additional analysis, such as incorporating bedfast and floating ice data, where required.

We thank the reviewers for their thoughtful and insightful comments, which helped improve our manuscript. We outline the revisions we plan on making to the manuscript in response to the reviewer feedback and address their specific comments and evaluations below. Broadly, we will implement additional clarification in the form of additional text and figures regarding the choices made in our methodology in the random forest, training labels, and implications of constraints in the methods and available datasets. We will further contextualize our results in relation to findings from previous studies by performing additional analysis, such as incorporating bedfast and floating ice data, where required.

We thank the reviewers for their thoughtful and insightful comments, which helped improve our manuscript. We outline the revisions we plan on making to the manuscript in response to the reviewer feedback and address their specific comments and evaluations below. Broadly, we will implement additional clarification in the form of additional text and figures regarding the choices made in our methodology in the random forest, training labels, and implications of constraints in the methods and available datasets. We will further contextualize our results in relation to findings from previous studies by performing additional analysis, such as incorporating bedfast and floating ice data, where required.  

Reviewer #3 comments: 

Nguyen and Lopez et al 

General Comments: The authors present ice phenology retrievals for over 30 000 lakes on the Alaskan Coastal Plain from the period of January 2019 to December 2023 using Sentinel-1 and Sentinel-2 to train and run a random forest machine learning model, which is validated using an Accumulated Degree Day (ADD) equation forced using ERA-5 data. The authors present clear, concise objectives that are statistically tested, with reasonable conclusions based on the assumptions of the input data. There are some issues with the Methodology, where statistical tests are presented in the results section but not described in the Methods. 

The main criticism which I raise with regards to the methodology is the question of ground truth or similar data. This study effectively produces an image stack of several Sentinel-1 polarimetric parameters (VV, VH, VV/VH) and retrieves ice phenology based on Sentinel-2 derived image labels (classifying as open water or ice from the Normalized Difference Snow Index (NDSI)). The resulting ice phenology retrieved from Sentinel-1 is then statistically compared against ADDT and ADDF, which are forced only with air temperature (if another variable is used, it is not described in an equation). There is no in-situ, or thermodynamic model that is used as validation in this study. The issue with using the thresholded NDSI values to discriminate open water and ice is that the NDSI does not perform well in the fringe seasons that the authors are looking to retrieve. Prior to snowfall, lakes that freeze with bare ice, or wind swept ice can be classified as open water. Concurrently, during the melt season, ice that has had the snowcover melt, exposing the clear ice beneath can be classified as open water, which can cause issues with the NDSI-derived ice flags. The authors need to show the reader with Sentinel-2 imagery the NDSI classification of freeze or melt to show its efficacy, but also quantify the accuracy of the ice labels derived from the NDSI. The other metric that is being used as validation is the ADDT and the ADDF, both of which are reliant on the air temperature at a 30km resolution. Firstly, there are higher resolution gridded air temperature reanalysis products (NLDAS is at 12.5km, and the NClimGrid is 5km) that cut reduce the variability in comparisons. Secondly, the ADDT and ADDF do not account for the thermal mixing that lakes undergo prior to freeze-up, introducing considerable error. I understand that with the sheer number of lakes in this analysis that a fully parameterized thermodynamic lake model may not be appropriate, but the authors should a) use reanalysis air temperature with a higher spatial resolution, and b) show that air temperature agrees with in-situ weather stations (e.g. Utqiagvik, Umiat, Inigok, Toolik Lake). This would better constrain the error associated with the ADDT/ADDF. 

We appreciate the reviewer’s general comments and criticisms of this work. We would like to clarify that the RF model used in this study is trained on the Sentinel-1 data itself (VV, VH, and the cross-polarized VV/VH ratio from both the lake interior and marginal landscape). The RF model uses training labels derived from Sentinel-2 NDSI as well as ERA5 reanalysis 2m air temperature data. Furthermore, ADDT and ADDF are only estimated from air temperature data alone, as there is no readily-available product for land surface temperature at this spatiotemporal resolution. 2m air temperature data is the most commonly-used alternative. ADDT/ADDF is an empirical index model that we use in place of in-situ or thermodynamically-modeled validation, and ADDT/ADDF is commonly used for lake ice studies (Stefan 1891; Richards 1964; Bilello 1980; Ashton 1985; Korhonen 2006; Arp et al. 2010; Kirchner and Hannam 2024). 

We agree that solely relying on NDSI to discriminate between open water and ice during the fringe season would not be advised, owing to misclassification of open water during freeze-up and melt from misidentified ice reflectance. For this reason, our strategy for training labels uses a combined NDSI threshold and temperature criterion check, as explained in the methods section 2.3, the response to reviewer #1’s second comment, and the response to reviewer #2’s third comment. Furthermore, the temperature data comprises 62% of training labels, as mentioned in the manuscript. This thresholding approach mitigates the confounding effects of ice phenology in the shoulder season mentioned by the reviewer. To make this clearer in the text, we will expand upon our explanation of the RF model training in the revision to the manuscript. Additionally, we agree that it would be illustrative to include a figure depicting Sentinel-2 imagery, corresponding NDSI value, and ice/water classification for randomly-selected lakes during both the ice-off and ice-on periods. We plan on creating this supplemental 6-panel figure to visually show the efficacy and accuracy of the RF model. 

While we appreciate the reviewer’s suggestions for the use of higher resolution climate reanalysis products, to our knowledge, and after an investigation in response to this comment, we note that the NClimGrid and NLDAS products suggested by the reviewer are available for the contiguous United States, not Alaska. NLDAS-3, which has planned Alaskan coverage, is currently in development, but at present only the meteorological model input (forcing) data is available, not the output data. For these reasons, we argue that the use of ERA5 is not only appropriate, but the best choice for this study. 

We appreciate the reviewer’s suggestion to compare ERA5 air temperature with in situ weather stations; the reviewer is correct that the 31 km spatial resolution of ERA5 cannot fully capture the spatial variability of air temperature or other climatic variables. However, a detailed validation of ERA5 reanalysis air temperature with in situ records is beyond the scope of this project. Instead, we will supplement the revision with a brief discussion of prior studies that have demonstrated excellent agreement between ERA5 and in situ weather stations across the Arctic (Clelland et al. 2024; Ibebuchi et al. 2024; Tarek et al. 2020; Pernov et al. 2024). For example, Clelland et al. 2024 has shown agreeance of ERA5 with in-situ data in Siberia, Ibebuchi et al. 2024 find ERA5 outperforms other reanalysis products for comparison with in-situ data across the contiguous United States, Tarek et al. 2020 found ERA5 hydrological modeling was comparable to in-situ data across North America, and Pernov et al. 2024 found a 0.94 R2  for temperature when comparing ERA5 to 17 weather stations across the Arctic. 

 The reviewer is completely correct in noting that the empirical ADDT/F model does not explicitly account for thermal mixing within the water column, potentially introducing error into the ADD model results. However, we agree with the reviewer’s remark that the sheer size and scale of this work precludes the implementation of a robust thermodynamic lake model. 

The authors present quite a large degree of variability in ice free period based on latitudinal gradient, 30 days per degree. There can be quite a bit of variance when modeling ice phenology using in-situ weather station data, including all the lakes within a 30km pixel to have one freeze-up and one melt date reduces this variability. This is likely why there is such a large spread shown in Figure 3. Additionally, based on the spatial clustering of high residuals (upwards of 60 days) in the lakes in and surrounding the Brooks mountain range, those lakes should not be included in the summary of ice free days per degree. The pattern of 30 days per degree is then discussed in Section 4.2. that it is six times that of the broader range of the Canadian Cordillera at 6 days per degree, “consistent with the heightened climate sensitivity of lakes near the northern limits of the seasonal ice zone”. Then the authors state “this sensitivity implies that continued Arctic warming could produce disproportionately large changes in ice phenology for these northern most lakes”. It appears that the baked-in error associated with large deviations in the southern region of the study site has influenced the 30 day per degree ice free days, and sweeping statements like this should omitted. 

In reference to the spatial clustering of residuals shown in supplemental figure 10, we agree that the lakes with high residuals (including up to 60 days) should not be included in the summary report of ice-free days per degree, as suggested by the reviewer. We will remove the relevant statements the reviewer highlighted. To address the necessary change in our summary report of ice-free days per degree, we will either remove lakes below a latitude threshold or process lakes below this threshold separately.  We will then report the updated summary report for ice free days per degree and make a clear distinction between lakes on the Arctic Coastal Plain and those in and around the Brooks Range. 

We will further remove the statement comparing the reported latitudinal/degree relationship to the Canadian Cordillera, and we will instead restrict our regional analysis to the subpopulation of lakes within the Arctic Coastal Plain and the Brooks Range foothills. 

The manuscript requires revision with respect to the validation (Sentinel-2, ADD) of ice covered regions both in data analysis and text to make the reader believe that the patterns being presented a function of ice dynamics in-situ. 

As mentioned above in our responses to this reviewer’s comments, we believe that our proposed inclusion of an additional figure showing randomly-selected imagery, an update to the summary report for ice free days per degree, and text surrounding these additions will address this revision request related to data analysis and text for validation of our method.  

Specific Comments: 

Page 2 Line 57: “the degree to which thermokarst lake morphology may influence ice phenology relative to estimates based on ADD metrics remains largely unknown” 

I don’t see how you could test how lake morphology affects ADD metrics. There are no inputs to the ADD that would be affected by morphometry. Other lake ice models have a mixing depth or bottom associated with them. 

This line was unclear and has been corrected. We fully agree that lake morphology would not affect ADD metrics which was the point of the sentence. While lake morphology may affect lake ice phenology, the ADD metrics would be unable to capture this behavior because they are independent of any morphometric properties: 

“As a result, the degree to which thermokarst lake morphology may influence ice phenology relative to estimates based on ADD metrics, which do not account for morphology, remains largely unknown.” 

Page 3 Lines 93-94: “we were able to retrieve phenology from 30,862 lakes (99.2%)… 

Why not all lakes in the dataset? What was the factor that inhibited the other 0.8% from being obtained? 

For the 31,108 lakes in the dataset, we did not retrieve phenology for 246 of them (0.8%).  Our determination of whether or not a lake had a successful detection required both ice-off and ice-on detection in any year within the five-year time period. The 246, or 0.8%, of lake either did not have an ice-on or ice-off detection. The reason that these lakes did not have either an ice-on, ice-off, or any detection may be due to there not being a minimum of two consecutive observations for ice or water for a given lake, the unlikely case in which ice-off would otherwise be detected outside of the window between April 1st to August 31st   or ice-on may have occurred prior to our beginning search date of September 1st. These time and date parameters set and called within the code may have excluded the ability to detect either or both of the ice phenology dates for these 246 lakes over the five-year period. We plan on further analyzing the reason why these lakes fell outside of the set constraints, their locations, and consistencies in morphometrics. We will include this information in the text, likely in the section 3.1 when addressing the ice phenology results. 

Page 4 Line 112: “the period January 2019 through December 2023” 

Why this time scale? Sentinel-1 and 2 have been in orbit since 2015, and after 2023. 

We selected this time scale based on the availability of the Harmonized Sentinel-2 MSI: MultiSpectral Instrument, Level-2A (SR) collection from Google Earth Engine, as described in the code. 2019 was selected as the start year, as the 2017-2018 data does not have global coverage in this collection. Furthermore, the start year of data for this collection is in 2017. Using 2023 as the final year for this data date range ensured that there were at least five years of data to retrieve lake ice phenology from, where we sought to define a nominal five-year time span over which to perform this analysis. We can include this information into the supplement if necessary. 

Page 5 Line 119: “-20 to -5 dB” 

For this location and others where you refer to the standardization of the symbologies – what justification is there for the max and min and the association made to the different scattering mechanisms associated with them? There needs to be a justification. 

The justification for the maximum and minimum for these ranges, as well as the different scattering mechanisms, are explained in Surdu et al. 2015 as well as Murfitt and Duguay 2020. We will add text explaining this in this paragraph and add these references as well. 

Page 5 Line 128: “we classified NDSI values > 0.4 as snow/ice”  

Does this mean that the input was then classified into a binary raster? 

This is correct. We use the 0.4 threshold on computed NDSI for the pixels within a lake-date pair to compute a binary classification. We then take the mean across all pixels within the lake-date pair of this binary classification to get an ice fraction (0-1). This ice fraction is then read into our RF label creation code that then labels an individual lake-date as either ice or water using a threshold of ice fraction > 80% being ice, and an ice fraction < 20% being water. We will update the text regarding the binary classification into section 2.2. Likewise, we will include information for the ice fraction thresholds to the first paragraph in section 2.3.  

Page 7 Line 194: This paragraph describes the ADD model 

Some equations are needed here to describe the ADD. Also, some justification here for using the ADDT and ADDF instead of a thermodynamic model, what the drawbacks are and the benefit of using it for such a large dataset. Where do you expect the air temperature to exhibit the highest error (e.g. mountainous regions). 

We agree with the reviewer that reporting the relevant equation for the ADD model would help motivate its use; we will update this paragraph to include the equations used in the ADD model. We will also incorporate a discussion of its prior use in similar studies of lake ice formation and decay, the strengths and limitations of the ADD model, and a justification for why we favor this method over a thermodynamic model for the regional scale study presented here. 

Page 7 Line 197, 198: “ranging from DOY 91 to DOY 242”…”DOY 244 to DOY 361” 

These are the endpoints of the DOY that the model has constrained as possible days for open water and freeze-up. The output having dates at the first or last potential should warrant a revision in the model to see why those dates are being flagged. 

This is true. The ranges reported here for ice-off and ice-on are exactly the windowed DOY we selected to search for ice-off and ice-on.  We originally chose these dates such that there are no instances in which ice-off is called after ice-on for a given year. To remedy this current issue, we plan to remove this search window and mark specific lakes in which ice-off is called after ice-on and will exclude them from the main reported summary statistics in this paper. However, we will note the number and percentages of lakes, as well as their spatial distribution and general morphometrics, that exhibit the physically impossible ice-off call after ice-on call. We will also report our ranges using either standard deviations from the mean or 5th and 95th percentiles to remove outliers. 

Page 8 Line 210: “After loss of Sentinel 1-B” 

This should be discussed in the data section – how did it affect your analysis? When was Sentinel 1-B lost? 

We will include information about the loss of Sentinel-1B in the methods for section 2.2, including the date in which the satellite was lost (December 23rd, 2021). The loss of Sentinel-1B affected our analysis in increasing the repeat time between acquisitions over our study area, thereby effectively decreasing the density of measurements after data transmission failure of the satellite. 

Figure 3: Figure 3 and in text does a good job at identifying that outside of observing a lake within a single Sentinel-1 revisit that the error increases substantially. There should also be a figure that shows the spatial pattern of comparison to the ADDT and ADDF estimates. Is this pattern consistent year to year, is the pattern consistent between freeze and melt? If consistent spatial patterns exist latitudinally, that may also drive the very large 30 day per degree reduction in ice free season. Figure S10 appears to show the spatial distribution of the residuals, showing that the largest residuals are in the southern most regions in the Brooks mountain range. 

We agree with the reviewer in that there should be a figure comparing the spatial patterns between the ADDT and ADDF model compared to the RF model calls of ice-off and ice-on. We will make a supplemental figure with 4-panel plots showing maps of ice-off and ice-on for both the ADD-derived and RF-derived ice phenology for each year. We will then include text describing the pattern consistency between years, ice-off and ice-on, and relation to the 30 days ice-free per degree change reported prior, although we also expect that the previously reported 30 days per degree relationship will change significantly after separating our analysis to lakes on the Coastal Plain and lakes in the Brooks Range foothills. 

Page 9 Lines 225-226: “we implemented a series of statistical analyses including …” 

These statistical tests should not be introduced in the Results section. They should be in a validation or statistical tests subheading in the methods section. 

  We will update the methods, specifically at line 194 (section 2.5) to include a brief overview of the statistical analyses employed in this study:  

“We employ the Kruskal-Wallis H-test, Mann-Whitney U-test, and Cliff's Delta (δ) to statistically evaluate any potential controls of latitude and lake morphometry on discrepancies with the RF model and ADD model estimates of ice-off and ice-off dates.” 

Figure 5: The X and Y axis need to be flipped. The independent variable is latitude, and the dependent variable is the ice-free duration. 

While it is true that the independent variable is latitude, and the dependent variable is the ice-free duration in this case, the choice of orientation for this figure was deliberate, such that latitude is displayed vertically as in a map-view of ice-free duration. However, we will include a statement in the revised text making explicit that the independent variable in this regression is indeed latitude and justify the choice of figure orientation.  

Page 15 Line 324: “ice phenology y difficulty” 

Typo 

Fixed. 

Page 15 Lines 327 – 329: “Likewise, the use of ADD to model thermal diffusion through the ice cover…” 

You should also mention that the ADD does not account for the thermal mixing that is required of lakes prior to initial skim ice formation. Larger/deeper lakes will take longer to mix and delay ice cover establishment. 

We appreciate this suggestion and will update the sentence to discuss how the empirical ADD model does not explicitly account for thermal mixing and what that may mean for ice phenology as a function of lake size.  

Page 15 Line 335: “models in figure 6” 

Should be “Fig. 6”. 

Fixed. 

Page 16 Lines 348-350 and Lines 355 – 358: There are two statements that are overly generalized here:  

“This sensitivity implies that continued Arctic warming could produce disproportionately large changes in ice phenology for these northernmost lakes, with consequently disproportionate effects on future green-house gas emissions (Walter Anthony et al., 2018).” 

“This indicates that the RF model detects a longer ice-free period for lower latitude lakes and a shorter ice-free period for higher latitudes relative to the ADD estimates, suggesting that the ADD model is more conservative in predicting ice-free duration at lower latitudes whereas at higher latitudes, the RF model is more conservative” 

The lakes in the lower latitudes are surrounding the Brooks mountain range, at higher elevations with variation in weather conditions. While this is acknowledged in the paper, the extrapolation that the 30 day per degree would potentially result in considerably larger greenhouse gas emissions does not take into account the different in elevation, situation, weather patterns, precipitation, etc. Latitude alone does not influence ice cover duration. This would be a good circumstance to test to see if the ERA5 air temperature agrees with the Toolik weather station located at the foothills of the Brooks range. 

The choice of ERA5 reanalysis air temperature was motivated in part because these reanalysis products do account for physical features and processes such as elevation, prevailing precipitation and weather patterns, etc. Previous studies have analyzed and validated the performance of different ERA5 temperature products for the Arctic region and have found strong performance of ERA5 data (Cao et al., 2020; Graham et al., 2019; Yu et al., 2021). However, a specific ERA5 air temperature analysis for the North Slope of Alaska, to our knowledge, remains outstanding. While we agree with the reviewer that this type of analysis would be beneficial, it is beyond the scope of this study to conduct these analyses, and such an analysis ought instead to be the subject of its own separate study.  

We do completely agree with the reviewer’s suggestion to conduct a subpopulation analysis of lakes within the Coastal Plain separately from lakes within the Brooks Range foothills, which we will do in the revision as discussed above. We will further modify the two specific statements identified here by the reviewer to discuss in greater detail the effect of elevation, situation, precipitation, etc. on ice phenology in addition to latitude.  

The statement comparing relative performance of the RF and ADD models as a function of latitude (i.e., lines 355-358) was intended to provide a generalized description of model performance over the study region. After performing the suggested changes to our analysis, we will append this discussion to make explicit that latitude is not the sole controlling factor that dictates ice phenology, and that systematic differences in elevation, precipitation regime, aspect, etc. can also affect ice phenology, and consequently might be expressed in the observed differences between ice-on/off date detections for the RF and ADD models discussed here. 

Page 17 Lines 387-388: “these larger and deeper lakes remain ice-free longer because they begin freezing later in the year likely due to their greater thermal inertia”. 

And the longer period associated with thermal mixing. 

We will modify this statement to the following: 

 “These larger and deeper lakes remain ice-free longer because they begin freezing later in the year due to a combination of longer periods of thermal mixing within the water column and their greater thermal inertia.” 

The authors mention the loss of Sentinel 1B – the density of measurements (i.e. repeat passes) must have been reduced as a result. I think that this should be mentioned and identify if it impacts the accuracy of the RF model before and after losing Sentinel 1-B. 

We agree that the reduction of density in the measurements from Sentinel-1 due to the loss of Sentinel-1B potentially affect the accuracy of our RF model. We will add text regarding the reduction of temporal visit time coarsening the granularity of the ice-on/off DOY estimates after its loss to the methods section. 

Supplemental Figures: Be sure that for the ice-off and ice-on dates in each of the figure captions you reference what analysis output the data is from. For instance, Figure S5: “Box plots of ice-off dates for phenology-detected lakes each year”. Is this the RF output? The ADD output? 

We appreciate the reviewer bringing this point to our attention. We will add additional text to the figure captions to explicitly state the specific analysis output data reported in each figure (e.g., S3, S5, S6, S7, S8, S10, S11, and S15).