Reviewer 1

Overview

This manuscript presents a study of thermokarst lakes in the Northern Territories, Canada, investigating hydrological, meteorological, and watershed controls on the lakes’ water balance. Twenty-five lakes were sampled five times over the spring-summer season in 2018 and analyzed for stable oxygen and hydrogen isotope ratios. The isotope data were used to calculate the evaporation-to-inflow (E/I) ratio and the isotopic composition of lake source water (δI). Four water balance phases were identified, associated with shifts in meteorological and/or hydrological changes. The E/I and δI were compared to meteorological and watershed data (e.g., precipitation amount and watershed area to lake area ratio (WA/LA)). Subsequently, the relationship between WA/LA and E/I was used to estimate the average E/I of 7340 lakes in this region, finding that the lakes are not vulnerable to desiccation, and that lakes with smaller WA/LA are likely to be more influenced by increased evaporation in a future warmer climate.

General comments

The manuscript presents an interesting approach to evaluate water balance changes and to estimate E/I, which is relevant for the scope of HESS. The manuscript is well written with an easy-to-follow structure and clear figures. However, there are some issues with the dataset size and associated assumptions that need to be more clearly outlined. I also have some questions concerning the chosen method to infer δI. I have listed my concerns and suggestions for improvement/clarifications below.  

• Thank you for your thorough review of our manuscript. We believe that we have addressed all your comments and have implemented most of the suggestions you made.
  
  
  
  

Precipitation dataset size: I have some concerns about the small size of the precipitation dataset used for the “isotopic framework”. Are 11 snow samples and 13 rain samples enough to produce a reliable LMWL and estimate δP, or would it be better to use a more conservative approach and use the GMWL? How well do the precipitation samples cover the annual range of precipitation isotopic variability in this region? And why not compare your data to GNIP data from Inuvik (from the 1980s)? Justify your approach to use local precipitation data (because they are most representative for the study area and cover the same period as your lake water dataset?). How does the uncertainty related to the precipitation data affect the uncertainty in the δ* and δI calculations? Some more information about the precipitation samples would be useful to clarify how well estimated δS, δR and δP are. Were the snow samples collected soon after snowfall, or can the “end of winter snowpack” have experienced post-depositional fractionation processes (e.g., sublimation) before sampling? Were the liquid samples event-based and/or representing all precipitation events from May to September 2018? Have the precipitation isotope data been amount-weighted? Please clarify.

• We have clarified the collection of snow samples. They were not collected during the winter as precipitation fell, but the end of winter snowpack was sampled by pushing a cylinder vertically through the snowpack at the end of winter, putting the sampled snow into a plastic bag to melt and then transferring the melted snow into a sampling bottle. In this way, we incorporate any fractionation processes that may have occurred as the snow was on the ground and measure the amount-weighted snow isotope composition by virtue of how we sampled the snowpack.
• We relied on the precipitation samples we collected to establish the LMWL because they most closely represent the input water to the lakes (compared to the GMWL), and we believe we have sufficiently sampled both rainfall and snowfall to determine the LMWL, at least for the year of 2018. Below we explain why we will remain with the LMWL we developed:
  • Unfortunately, there are only a few years of GNIP data available from Inuvik, from 1986-1989, with many months only having one year of data, and this is so long ago that we worry the values for δ_S and δ_R may have changed since the eighties as climate change alters global atmospheric circulation and reduces sea ice in the Beaufort Sea. Although, the equation of the GNIP LMWL line and our LMWL line are similar:
    • GNIP LMWL equation δ²H = = 7.3*δ¹⁸O - 3.5‰,
    • Our LMWL δ²H = = 7.1*δ¹⁸O - 10.0‰
  • By complete coincidence, another group has developed a LMWL for Inuvik using three years of data and 134 precipitation samples between 2015 and 2018 (Fritz et al., 2022). The LMWL, δ_R and δ_S values they calculate are very similar to ours:
    • Fritz et al. (2022) values for δ_R, δ_­S­, LMWL
      • Average, non-amount weighted δ_R: -16.9‰ (δ¹⁸O), -131.8‰ (δ²H)
      • Average, non-amount weighted δ_S: -23.8‰ (δ¹⁸O), -181.7‰ (δ²H) (Fritz et al., 2021).
      • LMWL: δ²H = 4*δ¹⁸O - 6.7‰
    • Our values
      • Average, non-amount weighted δ_R: -17.03‰ (δ¹⁸O), -129.54‰ (δ²H)
      • Average, non-amount weighted δ­_S: -24.61‰ (δ¹⁸O), -184.19‰ (δ²H)
      • LMWL: δ²H = 1*δ¹⁸O - 10.0‰
    • Given the close agreement in the LMWL, δ_R and δ_s values between our values and Fritz et al. (2022), we believe that redoing the analysis using the Fritz. et al (2022) data would not cause a meaningful difference in the results. We now discuss in the manuscript that our LMWL fits closely with the Fritz et al. (2022) and GNIP LMWL.
  • Post-depositional fractionation of snowfall was likely minimal (approx. 1‰ δ¹⁸O) – See Ala-aho et al. (2021), where they compare snowfall samples to corresponding snow layers at the end of winter for an Arctic tundra site (Figure 9a, 10a). Again, our agreement with Fritz et al. (2022) also suggests there was little change in snow isotope composition.
  • In regard to precipitation data affecting the uncertainty in δ^* and δ_I calculations: δ^* could impact δ_I through precipitation samples since we assume that atmospheric vapour (δ_As) is in equilibrium with the isotope composition of summertime precipitation (Equation A5). Since we seem to capture δ_R well when compared to Fritz et al. (2022), we think that we have likely calculated δ_I Also, lakes seem to plot around the LEL, which we calculate independently of lake water samples, suggesting that we have estimated δ^* accurately.
  • We have clarified the calculation of δ_S and δ_R and that they are the average of all samples without any weighting by precipitation amount.

Analytical uncertainty: Are duplicate measurements of every fifth sample enough to determine the analytical uncertainty? In many labs, each sample is injected multiple times (sometimes more than 10 times), the first replicate(s) discarded, and the rest used to calculate an average. How did you deal with drift and memory effects?

• By duplicate measurements, we do not mean duplicate injections. Each measurement consists of 8 injections of ~1000nL, of which the first two were discarded, with the remaining six injections averaged for the measurement.
• VSMOW and VSLAP standards provided by the instrument manufacturer were interspersed throughout the sample run to confirm the accuracy of the instrument.
• We have added this detail to the methods section.

 Isotope framework: I generally like the approach to present an “isotope framework”, but some clarifications are needed when comes to presentation of isotope data and the used terminology. One example is the inconsistent use of the delta notation (δ). Some examples of this are on line 44 where δ is missing when introducing water isotope analysis (which is the analysis of the stable isotope ratios δ18O and δ2H), and ‘δ-δ space’ which should be ‘δ2H-δ18O space’. There is also a parenthesis missing in equation 1, which should be written: δsample = (Rsample/RVSMOW-1) * 103. These fundamental things need to be stated correctly. Furthermore, it is not clear when δI is referring to δ18OI and when/if it refers to δ2HI as well?  

• We have corrected it to “δ²H-δ¹⁸O space”
• We now clarify that δ_I refers to both isotopes, and δ¹⁸O_I refers to just δ¹⁸ We use only δ¹⁸O_I values in the results and discussions since δ¹⁸O_I perfectly correlates with δ²H_I because δ_I values are bound to the LMWL.
• We have corrected the δ_Sample formula, however we have removed the “*10³” since a reviewer in another manuscript has pointed out to us that δ notation is not inherently in ‰ units (Coplen, 2011). This and the other manuscript will be a part of a PhD thesis, so we would like to keep the formulas consistent between the chapters.

Approach to calculate δI: Have you considered using the more recent ‘MWL source implementation’ method by Bowen et al. (2018) to model δI values? Using this approach, you supply a MWL equation (you could test the GMWL and your LMWL to see how much difference it makes) and a hypothesized LEL slope with confidence intervals, as well as your lake water δ2H and δ18O values with uncertainties.    

• We were vaugely aware of this approach but not considered using it. Our understanding of the benefits of the MWL source implementation approach is that it avoids biasing the δ_I of lakes when a regression through lake water samples is used to calculate the LEL, while there is correlation between evaporation-to-inflow ratios and δ_I. Since we calculate δ_I­ using a fixed point in δ²H-δ¹⁸O space that is calculated independently of lake water samples, and calculate a lake-specific evaporation line for each lake water sample, we believe that our method is not vulnerable to the errors that the Bowen et al. (2018) approach is meant to address. Our understanding is that the approach introduced by Bowen et al. (2018) is more geared towards large-scale studies covering larger areas than our region, and they actually encourage “model-based approaches” for site-specific studies like ours “We also believe that [a model-based] approach should be more frequently adopted in site-based studies as an alternative to regression…” (p1028, Bowen et al., 2018). Therefore, we will remain with our initial approach to calculating δ_I, but we will keep this other method in mind for our future research.

Specific comments

Lines 18-19: Rephrase the sentence starting “WA/LA strongly predicted average lake E/I ratio (R2 = 0.74)…” to “Lakes with smaller WA/LA tended to have higher E/I ratios (R2 = 0.74) because they received relatively less inflow. I think “strongly predicted” is an exaggeration. 

• We disagree that it is an exaggeration – a single variable explaining three quarters of the variability in a complex environmental system is uncommon. However, we changed the wording anyways.

Line 40: evaporation-to-inflow ratios. Use this wording throughout the paper (sometimes it says evaporation to inflow, sometimes evaporation/inflow ratio).

• We have changed the wording throughout the manuscript to be consistent.

Line 44: Change to “In several studies, stable water isotope (δ18O and δ2H) analysis…” and add a couple of references.

• Changed to: “In several studies, water isotope (δ18O and δ2H) analysis has been the primary method used to efficiently characterize the water balance of a large number of thermokarst lakes because the isotope composition provides an integrated measure of influential hydrological processes (Gibson, 2002; Edwards et al., 2004; Turner et al., 2014; Narancic et al., 2017).”

Line 60: Explain the concept “snowmelt bypass” here. Now it isn’t explained until in the methods (lines 88-89).

• Added: “Snowmelt bypass occurs when less dense (~0°C) freshet runoff flows underneath lake ice and passes through a lake without mixing with and replacing the deeper and denser (<4°C) lake water.”

Line 74: Remove “All lakes we selected were either headwater lakes or downstream of lakes that we sampled”, since this information is repeated on lines 76-77.

• Change implemented.

Line 87: How much did the first day of ice-free season vary between the lakes? 

• Added: “June 15 was chosen as it marked the first day of the ice-free season for most lakes and was intended to capture the influence of snowmelt bypass, however the southernmost lakes (Lakes 7-12) became ice-free on June 7, and some of the northernmost lakes (Lakes 49-55) became ice-free on June 17 or June 18.”

Lines 88-89: Move the definition of snowmelt bypass to the introduction and rephrase it. This sentence does not read well and needs to be clarified.

• Change implemented.

Line 97: Explain “end-of-winter snow”. Does this mean that you sampled the snow in the end of winter soon after it fell, or that you collected a core of snow accumulating over a longer period? If the latter, how do post-depositional processes impact the isotope values?

• There is likely minimal change in snow isotope compositions after snow is deposited, as we describe above. See Ala-aho et al. (2021).
• Also as we describe above, we have clarified that the snow samples are taken using a vertical core through the end-of-winter snowpack.

Line 98: Are δS and δR amount-weighted? Change to δSnow and δRain throughout the text to match the terms used in Figure 3 (or change to δS and δR in the figure).

• They are not amount weighted; we have made the terms similar between the text and figures.

Line 107: Equation misses parentheses around (Rsample/RVSMOW - 1)

• As described above, we have changed the formula to omit the x10³ _­now as δ notation is not inherently in ‰ units.

Lines 113-114: When introducing the “fundamental linear relationships”, describe the global relationship (i.e., the GMWL) as well. This could also be added to Figure 3.

• We now describe the GMWL and have added it to Figure 3.
• Added:The isotope composition of precipitation from anywhere in the world tends to plot closely to the Global Meteoric Water Line, a strong linear relation in δ¹⁸O-δ²H space which is defined as δ²H = 8.0 ∗ δ¹⁸O + 10 (Craig, 1961)."

Line 114: Change ‘δ-δ space’ to ‘δ2H-δ18O space’ (also on line 122).

• Change implemented.

Line 116: Is δP the average between the two δS and δR values, or the average of the full range of rain and snow values?

• It is the average of the two (δ_S and δ_R)values, otherwise the number of δ_S and δ_R samples would weight δ_P. We have clarified this now in the text.
• Added: “The average isotope composition of rainfall (δ_R) and snowfall (δ_S) were calculated by averaging the isotope composition of all the rainfall and snowfall samples we collected respectively, with δ_P representing the average of δ_S and δ_R.”

Lines 116-122: As it is explained now, it is not clear how the LEL was defined. Looking at Figure 3, the slope of the LEL is different between δP and δSSL compared to between δSSL and δ*? It also says that δSSL is located along the LEL. Which line is the LEL equation referring to? Why is it not a straight line from δP to δ*? Please clarify.

• We double checked our calculations and found an error: we had calculated δ­_SSLfor a lake fed with δ_R instead of δ_P, which caused the bend in the line. We didn’t initially notice this error because there is usually a slight bend in the line associated with uncertainties in the fractionation equations for δ²H, which influences δ­_SSL and δ^* The line has now been fixed by calculating δ_SSL using δ_P as the source water.
• We have now labelled the line referencing the LEL on Figure 3.
• Modified the description of the LEL for clarity: “Evaporated waters tend to plot along a Local Evaporation Line (LEL), which can be defined independent of measured lake
  
  water isotope compositions. Using this approach and for the case of a lake fed by waters of mean annual isotope composition of precipitation, the LEL was anchored at δ_P and the theoretical isotope composition of a lake at the point of total desiccation (δ^*), which is dependent on air temperature, relative humidity and the isotope composition of atmospheric moisture in the region (Gonfiantini, 1986). Along the LEL exists a useful reference point where the amount of evaporation a water body is experiencing is equal to the amount of water input, defined as δ_SSL (steady-state lake water isotope composition of a terminal basin).”

Line 122: What does SSL stand for? Add reference for δSSL definition.

• Steady state lake, we have added this now: “defined as δ_SSL (steady-state lake water isotope composition of a terminal basin).”

Lines 126-128: Describe the δI calculation here as well (not only in the appendix).

• We now split the paragraphs describing the calculation of δ_I and E/I, and more clearly describe how δ_I is calculated.

Line 131: “where δL is the isotope composition of the lake water and δE is the isotope composition of evaporated vapour from the lake (Gonfiantini, 1986)”

Line 136: cumulated rainfall?

• Total rainfall amounts within the periods of time being discussed.

Lines 150-151: What do you mean by “data were transformed if the distribution was non-uniform”?

• If the distribution of the data on a histogram was non-uniform, a mathematical function (e.g. ln) was applied to the data in order to transform the data into a uniform distribution.
• Changed to: “The distribution of each variable was plotted on a histogram and data were mathematically transformed if the distribution was non-uniform.”

Line 152: Rephrase this sentence (removing the “strong”) and add the R2-value.

• Change implemented.

Line 170: This is the first time you use δ18OI. Does δI mentioned throughout the paper refer to δ2HI and δ18OI collectively, or δ18OI? Please clarify.

• We now clarify this when we introduce δ­_I in the methods section
• Changed to: “Isotope compositions of the lakes were used to estimate the average isotope composition of source water (δ_I, or δ¹⁸O_I when referring to just δ¹⁸O)”

Line 170: “We observed distinct shifts in lake water isotope composition along the LEL…” – This is not easy to read from Figure 3. Do you mean for individual lakes, or the lakes in general? It could be interesting to indicate which lake is which, to be able to compare the different lakes’ responses during P1-P4. This could for example be done (in a supplementary figure?) by giving the lakes different colors, and assign each sampling date a different symbol, and/or by drawing lines between data points from the same lake.    

• We explored some of the lake-specific changes in water isotope composition in our data analysis, but it becomes very messy to plot all lakes at once, and most lakes followed the general pattern over time that we describe in the discussion (the four water balance phases). There are some interesting patterns for a couple lakes and for clusters of lakes where we have some upstream and downstream water samples, but overall, our goal with this paper was to generally describe the variability in lake water balances over space and time. Maybe in a future manuscript we will explore more lake-specific water balance evolutions, and the impact of upstream lakes (and their watershed properties) on downstream lakes.
• Changed to: “In general, we observed distinct shifts in lake water isotope composition…”

Lines 174 and 180: Remove “very”.

• Change implemented.

Lines 178 and 184: Change “compositions” to “values”

• Change implemented.

Lines 193-194: “On September 3, some lakes plotted close to the LMWL, indicating that their waters had experienced negligible amounts of evaporation (Figure 3).” Does this mean that the same lakes plotted on the LEL before, and that the high precipitation amounts during P4 “reset” the lake water to be closer to δI by removing the old evaporation signal? 

• Yes, their waters had no evaporation signal. These lakes likely had large WA/LA ratios and short water residence times as a result, making their fluctuations in E/I more sensitive to short-term changes in meteorological conditions.

Line 196: Two thirds of the lakes.

• Change implemented.

Line 221: Remove “results” before “appears”.

• Change implemented.

Lines 245 and 246: Change “cooler” to “lower”.

• Change implemented.

Line 278: Earlier you mention only 5 downstream lakes (e.g., in Table 1). Where did the 6th lake come from (also presented in Figure 5)?

• We had incorrectly labelled Lake 12 in our table as a headwater lake when it is in fact a downstream lake. We have corrected the table and text accordingly.

Line 320: evaporation

• Change implemented.

Line 321: led

• Change implemented.

Line 329: remove “as high as”

• Change implemented.

Line 357: Is δPs the value referred to as δP in Table A1, or where do you present that value?

• δ_Ps is summertime precipitation, while δ_P is all precipitation. This is clearer now in our revised Table A1, where we describe all the parameters in the isotope framework and how they are generated.

Figure 1: Add overview map (e.g., Canada), and add black triangles and red square to legend.

• Change implemented.

Figure 2: Add reference to meteorological data. Clarify that you show the cumulative precipitation amounts. The sample days are shown by vertical dashed, not dotted lines. Mean daily air temperatures are indicated by horizontal dashed lines. 

• Changes implemented.

Figure 3: Please explain all elements in the caption, e.g., that the δSnow, δRain (which are called δS and δR in the text) and δP are averages (it looks like δRain displays the median, but maybe the average and median are very close?), what the box and whisker plots show, and that the numbers refer to sampling dates. It would also be good to add the GMWL for reference. Which line does the LEL equation refer to? The line between δP and  δSSL and the line between δSSL and δ* have different slopes. Why? Change color for the snow samples or the May 1 samples, as it is confusing that both are displayed in grey. You could also use different symbols for precipitation samples and lake samples, to make them easier to differentiate.

• We have added descriptions for all elements shown in Figure 3. The colour of snow samples has been changed to purple, and the symbols for precipitation have been made different to the symbols for lakes. A legend has also been added for all the different isotope samples.
• As we mentioned above, we had made an error in calculating δ_SSL. The LEL is now a straight line. The calculation of δ_SSL is now clarified in Table A1.

Figure 4: It says four sampling dates in the caption but should be five. What do you mean by “lake-specific” change in E/I in (c)? Why not also present a panel with the measured δ18O values, to see how much these values differ from δ18OI?

• We have corrected the number of sampling dates.
• We have clarified what “lake-specific change in E/I” means, we calculated the difference in E/I between sampling points for each lake, and then show the distribution of those values in (c). We added the equation for this to the figure caption.
• We are not sure of the reason to compare δ¹⁸O to δ¹⁸O_I, we already show the distribution of δ¹⁸O values for each isotope sample type in Figure 3, and one can get an idea of δ_I values by observing how lakes plot above or below the LEL.

Figure 6: Are the temperatures and precipitation values in from the year of lake water sampling or the 1980-2020 values? And how were the curves in the three lowermost panels generated? The E/I values at the sampling dates are not the same as presented in Figure 4b? And how were the values between the sampling dates interpolated? Please clarify.

• We have now clarified in the figure caption that this is a conceptual diagram and not pure data, although the temperature and precipitation data do represent 1980-2020 averages. We have redrawn the E/I and δ_I lines to better represent our data.

Table 1: Change Polygon Extent to Ice-wedge Polygon Coverage.

• Change implemented.

 Table 2: Clarify in caption that you mean δ18OI. What do you mean by “adjusted R2”?

• Adjusted R² is used in cases with multiple explanatory variables, in an attempt to control for the natural tendency for R² to increase as additional variables are added. In R, the Wherry formula is used to calculate adjusted R²: R²_adj = 1 – (1 – R²)*((n-1)/(n-p-1)), where n is the sample size and p is the number of explanatory variables (i.e., predictors). When there is only one explanatory variable, R² ≈ R²_adj. We have added this to the caption.

Reference

Bowen, G.J., Putman, A., Brooks, J.R., Bowling, D.R., Oerter, E.J., Good, S.P., 2018. Inferring the source of evaporated waters using stable H and O isotopes. Oecologia 187 (4), 1025–1039. https://doi.org/10.1007/s00442-018-4192-5.

• References
  • Ala-aho, P., Welker, J. M., Bailey, H., Højlund Pedersen, S., Kopec, B., Klein, E., Mellat, M., Mustonen, K., Noor, K., and Marttila, H.: Arctic Snow Isotope Hydrology: A Comparative Snow-Water Vapor Study, 12, 150, https://doi.org/10.3390/atmos12020150, 2021.
  • Coplen, T. B.: Guidelines and recommended terms for expression of stable-isotope-ratio and gas-ratio measurement results, 25, 2538–2560, https://doi.org/10.1002/rcm.5129, 2011.
  • Fritz, Michael; Wetterich, Sebastian; McAlister, Joel; Meyer, Hanno (2021): Stable oxygen (δ¹⁸O) and hydrogen (δD) isotopes in precipitation from Inuvik in the western Canadian Arctic (2015-2018). PANGAEA, https://doi.org/10.1594/PANGAEA.935027 [data set]
  • Fritz, M., Wetterich, S., McAlister, J., and Meyer, H.: A new local meteoric water line for Inuvik (NT, Canada), 14, 57–63, https://doi.org/10.5194/essd-14-57-2022, 2022.

The manuscript employs the well known isotope mass balance formulations to compute the water balance, namely, the evaporative water loss to total inflow (E/I) ratio in 25 Thermokarst lakes in Canada. The authors subsequently compare the E/I to the lake and catchment geometry and cover [Catchment Area to Lake area ratio, Lake Depth, Vegetation, gradient, Lake area etc].  Interestingly, the E/I ratio of the 25 lakes correlate well (R2=0.74) with the Catchment/Lake area ratio, allowing the authors to compute the E/I ratio for 7340 Thermokarst lakes without needing the hydrologic or isotopic measurement in these lakes.

The 25 lakes have similar range in depth (0-4 m) dictating the volume of the water in the lakes is proportional to their area (because of similarity in lake depth). This would mean the underlying mechanism for the good correlation between E/I and Catchment/Lake area ratio is that the later is a good surrogate for Catchment/Lake Volume and thus to residence time (Volume/Inflow). The major question now is what would happen if the lakes (25 lakes) vertical dimension was significantly variable and if the 7340 lakes have variable depth properties? The validity of the assumption that the lakes have similar depth property must be clearly reflected in the manuscript or the uncertainties of using the approach or extending the approach beyond these lakes needs to be highlighted in the manuscript.

• Thank you for your comments, variability in lake depth, and therefore volume, is not something we had considered. We have done some simple statistics using lake depth data published in a previous study of 34 different lakes in this region by Pienitz et al. (1997). The median lake depth was 3.0 m, average of 5.0 m, and a range between 1.5 to 18.5 m. We now added in the discussion section of the manuscript that our E/I estimates may underestimate E/I for these deeper lakes, as increasing lake depth likely increases residence time and therefore the opportunity for lake water to evaporate. It would be interesting for future work to compare these deeper lakes to shallower lakes, with respect to E/I.

Other detailed comments (13 comments) are provide in the attached PDF.

• Replies to your pdf comments are in the attached pdf.

References

• Pienitz, R., Smol, J. P., and Lean, D. R.: Physical and chemical limnology of 59 lakes located between the southern Yukon and the Tuktoyaktuk Peninsula, Northwest Territories (Canada), 54, 330–346, https://doi.org/10.1139/f96-274, 1997.

• Replies to your comments in the PDF are in the attached PDF.
  
  
  
  

The authors investigated how lake and watershed properties and meteorological conditions influence the water balance of thermokarst lakes in the tundra uplands east of the Mackenzie Delta in the Northwest Territories of Canada. They sampled 25 lakes for isotope analysis five times in 2018, beginning before snowmelt on May 1 and ending on September 3 to calculate the ratio of evaporation-to-inflow (E/I) into the lakes and estimates the average isotope composition of lake source water (δI). They indicated the range of average E/I ratios compared very well with watershed area to lake area (WA/LA) parameter (R2 = 0.74). Furthermore, they used this relationship to predict the average E/I ratio of 7340 lakes in the region, finding that lakes are not vulnerable to desiccation in the region. The authors have a good data and other measurements that could compliment isotopic computations to establish a surrogate.

While the general idea for finding a surrogate parameter which is easier and less expensive to measure or used to estimate E/I ratios laudable, the competent application of the time-tested evaporation model (water -isotope -mass balance method) used by the authors to justify this surrogate build up is doubtful. There are several assumptions that were made or glossed over in the paper that are difficult to justify, and these make it difficult to accept the results as they are presented. A major update of the assumptions and re-evaluation of the E/I computations is suggested. 

• Thank you for your comments and suggestions on our manuscript. While we agree with many of your comments, we believe that you may have misinterpreted our approach to calculating E/I based on the presence of a singular LEL on Figure 3. The LEL on Figure 3 was placed there as a reference case only, and we have in fact calculated lake-specific evaporation lines for every lake water sample, as outlined in Appendix A. We agree with the suggestions you made under the assumption we had used a singular LEL approach, but we believe our actual methods address most of the concerns you had. We have attempted to clarify our isotope framework methodology in the revised manuscript by implementing your other suggestions (e.g., adding a table with all parameters used in the isotope framework).

Below are a few suggestions to improve the paper.

The authors indicated that the 25 lakes cover a range of watershed sizes, surface areas, among other characteristics (which include vegetation, land use type, rolling hills etc) and that the lakes were sampled as either a headwater lakes or a downstream lakes (line 70-75). The above information suggests the ff:

1) that the isotopic behaviour/ responses to evaporation from these lakes may vary widely and hence there is the need to individually and separately evaluate the LEL, h, delta P, delta A, and finally E/I and not lump them together as shown in Figure 3; I encourage the authors to study Bam and Ireson 2019 https://doi.org/10.1016/j.jhydrol.2018.12.032  since this work is similar.

• Although we display a LEL on Figure 3, we did not used a single evaporation line to calculate E/I for all lakes. We define the LEL (Local Evaporation Line), as being between δ_P and δ^*, defined in the Figure 3 caption and in lines 118-119 “… the LEL was anchored at δ_P and the theoretical isotope composition of a lake at the point of desiccation (δ*)…” The LEL is used as a reference component in Figure 3, to show whether lakes tend towards rainfall- (above the LEL) or snow-sourced (below the LEL) waters.
• We have, in fact, treated every lake water sample individually when calculating δ­_I and E/I, by generating a lake-specific evaporation line for every lake water sample, i.e. calculating δ_E and δ_I separately for each lake water sample, as outlined in Appendix A, Lines 368-370: “The isotope composition of the lake-specific input water (δ_I) was calculated following Yi et al. (2008), where δ_I is estimated as the intersection of the LMWL and the lake-specific LEL, defined as the line between the measured isotope composition of the lake (δ_L) and δ_E, which is the isotope composition of vapour evaporating from the lake.”.
• We have used single values for h, δ_P, and δ_As because these lakes were all sampled in the same year, and are situated in a small enough region that it is reasonable to assume that the atmospheric conditions they experience are similar enough to warrant using a single value for h, δ_P and δ_As­. We do see evidence there may be slight differences in δ_P latitudinally as the distance to the ocean changes, however this should not affect our isotope framework greatly.

2) there is a need to account for the inflow and outflow fluxes to these lakes, I did not see this accounting process or equations, neither was there a stated assumption that these fluxes have been considered negligible within the case and the period and why;

• It was not our goal to quantify individual lake water balance components (inflow, outflow, evaporation, precipitation) for each lake. We do not assume inflow and outflow to be negligible, but we assume that lakes remain at a relatively steady volume (inputs ≈ outputs) over the period that lake water isotope samples integrate a hydrological signal (water residence time), see Line 375: “…assuming that lakes are well mixed and in hydrological and isotopic steady state.” We believe it is reasonable to assume that lakes remain at roughly the same volume given that no lakes appear to experience desiccation and therefore are always fed enough water to remain above their sill level, maintaining a relatively stable water volume.

3) If even we assume that these lakes are well-mixed throughout the period (this assumption was not made though), precipitation addition and output fluxes are negligible and hence E/I ratios could be computed; we know a range of other factors that do significantly affect evaporation rates from lakes and hence lake isotope composition and ultimately the E/I. Among these are the exposed surface area, depth, the geometry/ basin bathymetry, vegetation cover etc., are all these the same of the 25 lakes? I believe the authors have sufficient data to account for these parameters, since the ultimate goal is to find a surrogate.

• As stated above, we do assume lakes to be well-mixed, and we do not assume outflow to be negligible, but we assume that the lakes are in a steady-state condition.
• Precipitation inputs are also integrated into the “inflow” portion of the Evaporation-to-inflow ratio, by nature of the isotope methodology.
• We have accounted for other parameters that have the ability to impact E/I and δ_I. As described in lines 139-151, we tested many lake and watershed attributes (all listed in Table 1) for significant correlation with E/I and δ_I¸ however we only showed statistically significant relationships in Table 2. We clarify this point now in the Table 2 caption. “Only significant (p<0.001) relationships are shown in the table.”

4) the choice or search for delta P (initial lake water isotope composition) using the intersection point of the LMWL and LEL is not clear. We see in Figure 3 that both May (and Sept ) samples plotted on or above the LMWL and if for anything the isotopic compositions of the May /June samples could be used as the starting point for those lakes to estimate the E/I between May and Sept or August. Extrapolating the LEL to the bottom to meet the LMWL to find initial source water has been proven to be unrealiable (sese Bennettin et al., 2018, https://doi.org/10.5194/hess-22-2881-2018) .

• We did not use δ_P as the initial lake water isotope composition in our methods. We define δ­_P as the average isotope composition of precipitation, from which we anchor the LEL which we provide as a reference only. We calculated an “initial lake water isotope composition” (i.e. lake waters before evaporation took place) as δ_I, as described in our response to comment #1. See Lines 368-370: “The isotope composition of the lake-specific input water (δ_I) was calculated following Yi et al. (2008), where δ_I is estimated as the intersection of the LMWL and the lake-specific LEL, defined as the line between the measured isotope composition of the lake (δ_L) and δ_E, which is the isotope composition of vapour evaporating from the lake.”. Since we are not computing a linear regression through lake water samples to define the lake evaporation line, our approach to calculating δ_I is not vulnerable to the biases described by Bennettin et al. (2018).

I would recommend that since five (5 ) samples were collected from each of the 25 lakes from May to Sept, it will be useful an LEL is constructed for each of the lakes and used together with the rainfall isotopes and temperature measurements to estimate the dA for each of these lakes (see Bam and Ireson 2019; Bennett et al., 2008; Gibson et al., 2008a,b,2016; Gibson and Reid, 2014; Skrzypek et al., 2015). Use the May/June sample isotope composition as the initial lake water isotope to estimate the E/I for the period (May/June to Sept where feasible).

• The studies you have referenced seem to use variable δ_A values because they use multiple years of water sampling, or because of large distances between lakes where one would expect differences in the atmospheric vapour source. Since all our samples are collected in the same year, and we assume that lakes experienced reasonably similar atmospheric conditions given their proximity to one another, we believe it is reasonable to use a single value for the isotope concentration of atmospheric vapour (δ_As).

There cannot be any E/I estimate outside the ice-free period, water fluxes lost during the dry winter periods originate from the sublimation of frozen ice over and snow, which covers on the surface of the lakes. The lakes are closed to evaporative losses prior to the open-water season. 

• We do not claim that the water samples taken from lakes on May 1 represent the water balance conditions on that exact date, see L83-85 “The May 1 samples also represent the hydrological status of lakes during the freeze-up period of 2017, since virtually no hydrological activity occurs during the winter months due to the complete freezing of the soil and lake surface.”. We use samples from May 1 because it allows us to evaluate how the freshet has influenced lake waters.

References

Skrzypek, G., Mydłowski, A., Dogramaci, S., Hedley, P., Gibson, J. J., and Grierson, P. F.: Estimation of evaporative loss based on the stable isotope composition of water using Hydrocalculator, 523, 781–789, https://doi.org/10.1016/j.jhydrol.2015.02.010, 2015.