Thank you very much for reviewing our manuscript.  We are happy to learn your positive evaluation highlighted by the first sentence, “approach is relatively novel, a mathematical model used to represent this process is a good idea”.   Here, we provide responses to a few key points to enhance the information available to reviewers during the discussion round.  The first point is on the validation:

• Are the results sufficient to support the interpretations and conclusions?
  
  Not really… as noted above, the model derivation is nice, but comparing the model results to one measured lake is a little worrying.

There are only a very limited number of talik depth measurements under an isolated lake in a continuous permafrost – many of these examples are single drill points. The example at Peatball lake is, to our knowledge, the only quasi-3D dataset available in the Arctic.  We strongly believe that the TEM sounding survey of the Peatball Lake is the most comprehensive dataset, and therefore most appropriate for this comparison.

• Some minor issues noted below, and some jargon and unnecessarily complex language used to describe especially mathematical derivations.
• Eq 13 I think it may make more logical sense to present this in the opposite direction - the integral along the phase boundary (line) is not something that I can interpret easily or can be visualized, whereas something more like the flux across the phase change surface, or the volume integral of the total energy in the lake is more easily interpreted. I would start with the resulting equation and state which theorem (Stoke’s?) is used to get the first equation. Importantly providing a physical interpretation (in more simple terms) of what each expression (start and derived result) means and how it is useful and what it tells us about the system. This would greatly increase the utility of the work for those who are less interested in the mathematics and more interested in their application.

Unfortunately, as this study covers the different fields of studies, some jargon is unavoidable – indeed, we will work to improve this if given the opportunity to submit a revision.  We would like to highlight that the methodology itself is novel in many aspects.  We propose more explanation in the upcoming revision stage so that physicists, earth scientists as well as mathematicians can understand the content better.  One possible reason for confusion may be too many appearances of “Euler” and “Lagrange” in terminology across the related fields with slight variations.

For the sake of simplicity, the original manuscript focused on the mathematical technique, which appears as “Euler equation in the calculus of variation”.  However, we propose to enhance explanation with the physical context beyond the Newtonian mechanics, which hopefully helps readers can understand the methodological background of this study.  Additionally, the keywords are available in the Wikipedia, which provides fairly accurate explanations usually in plain language for readers who do not have a proper background.

https://en.wikipedia.org/wiki/Lagrangian_mechanics

“Euler–Lagrange equation” may replace “Euler equation in the calculus of variation”.

https://en.wikipedia.org/wiki/Euler%E2%80%93Lagrange_equation

We propose the following additional paragraph placed at the beginning of Chapter 2 Theory.

“This study uses the stationary action principle (the principle of least action) based on the Lagrangian mechanics, which generalizes the classical Newtonian mechanics.  The action is defined as the integral of the Lagrangian, which consists of kinetic and potential energies of the system.  In this permafrost application, the Lagrangian simply becomes the potential energy due to absence of kinetic energy.  The variational principle, the main tool in Lagrangian mechanics, can indeed derive the equations in the Newtonian mechanics. One of the related research topics using such a variational principle is a phase boundary propagation that can be analyzed by the phase field model or diffusion-interface model (Cahn and Hilliard, 1958; Cassel, 2013).  This model explains the diffuse phase boundary without surface tension, which appears in Newtonian interfacial physics between liquid and gas but irrelevant for liquid-solid interface. According to the second law of thermodynamics, monotonical decrease of the free energy is required for a non-negative entropy production (Singer-Loginova and Singer, 2008). This requires the time rate of change of the phase boundary to be expressed by the functional derivative of the free energy functional, which corresponds to the basin integrated energy flux in the permafrost thaw problem.  This study directly and analytically solves the Euler-Lagrange equation based on the stationary action principle rather than the entropy functional used in the phase field method.”

Additionally, Equation (15) uses the method of Lagrange multipliers (https://en.wikipedia.org/wiki/Lagrange_multiplier) which is a common tool in the machine learning field (e.g. maximum entropy principle) for optimization.  We plan to indicate the name of the method (method of Lagrange multiplier) for readers to understand the physical interpretation.

We will also address all points in the next stage.

References

Cahn, J. W., & Hilliard, J. E. (1958). Free energy of a nonuniform system. I. Interfacial free energy. The Journal of chemical physics, 28(2), 258-267.

Cassel, K. W. (2013). Variational methods with applications in science and engineering. Cambridge University Press.

Singer-Loginova, I., & Singer, H. M. (2008). The phase field technique for modeling multiphase materials. Reports on progress in physics, 71(10), 106501.

Thank you very much for reviewing our manuscript.  We prepare the final responses for each individual point.  The parts overlapping with our last responses to Reviewer 1 are still included but simplified.

Reviewer 1

• Approach is relatively novel, a mathematical model used to represent this process is a good idea.

In my view, though the model is nifty and mathematically neat (except some small areas where further details/assumptions should be stated to ensure applicability), there is essentially no field data that adequately corroborates the result. The lake used to ‘validate’ the model fits reasonably well along one depth problem, but the entire bathymetry of the lake is not represented, and it is stated that this lake may actually be composed of two different thaw features. It would be of great added value to compare the model to multiple lakes in an area where wind and uneven segregated ground ice distribution are not a large factor to see whether the lakes formed do indeed follow the model. Also to predict time series thaw in lakes as compared to model predictions. As it is, the comparison to data concerns me as it looks like validation, but the model would be more ‘trustworthy’ if stated in theoretical terms as opposed to in direct application.

Thank you for your positive overall evaluation.  Our point-by-point response will be presented below.

• For the most part, though some of the derivation needs a bit of clarification; likely the results hold, but as a reader I struggled to follow some of the steps. Sometimes the assumptions previously stated should be reiterated (especially in the discussion) to ensure the readers understand the limitations of this model.

We will include more complete assumptions for the theoretical equation in addition to the quasi-steady state approximation.  They are as follows:

The material of permafrost and talik is fully saturated with ice and water, respectively.  The thermal constants (thermal conductivity, latent heat, and thawing temperature) are constant and isotropic. Change in volume of water on thawing and freezing is negligible.

• Are the results sufficient to support the interpretations and conclusions? Not really… as noted above, the model derivation is nice, but comparing the model results to one measured lake is a little worrying.

There are only a very limited number of talik depth measurements under an isolated lake in a continuous permafrost – many of these examples are single drill points. The example at Peatball lake is, to our knowledge, the only quasi-3D dataset available in the Arctic.  We strongly believe that the TEM sounding survey of the Peatball Lake is the most comprehensive dataset, and therefore most appropriate for this comparison.  We will indicate the data availability in the revised manuscript.

• As noted above, some detail in derivations is missing. The manuscript is also quite long, and so I was unable to review the supporting information as I would like and cannot comment on the quality of the derivations therein.

We will try to improve the mathematical derivation in the next stage.

• The authors give proper credit to related work and clearly indicate their own new/original contribution? Yes, though the introduction may need a few additional references for some ‘obvious’ concepts which are clearly not the author’s own ideas, but generally accepted in the field.

Thank you.  We plan to improve the references in the introduction.

• Some minor issues noted below, and some jargon and unnecessarily complex language used to describe especially mathematical derivations.

We will reduce the number of jargons by choosing better terms. However, as this study covers the different field of studies, some jargons are unavoidable.  More explanation, which might include additional jargons to readers in other fields, was added in the revised manuscript so that broader range of readers can understand the content better.

• For the most part, some missing units and inconsistent use of symbols detailed below

We will provide the missing units.

• I am unsure of the comparison with field data - I think either this section should be expanded to include more sites, removed (which I am sure the authors agree would detract from the merit of this contribution), or perhaps re-phrased as an example application of this new method and not a test of the method proving its efficacy.

It is hard to obtain the talik depth measurements under an isolated lake in a continuous permafrost as stated above.

• L 23 … the Euler equation and the calculus of variations

The original manuscript focused on the mathematical technique, which appears as “Euler equation in the calculus of variation”.  We propose more complete explanation with the physical context beyond the Newtonian mechanics, which hopefully helps readers to understand the background of this idea in some extent.  “Euler–Lagrange equation” replaces “Euler equation in the calculus of variation” in the revised manuscript. 

We add the following paragraph at the beginning of Chapter 2, Theory.

The approach used in this study is based on the Lagrangian mechanics, which generalizes the classical Newtonian mechanics, using on the stationary action principle (the principle of least action).  The action is defined as the integral of the Lagrangian, which consists of kinetic and potential energy of the system.  In this application, the Lagrangian simply becomes the potential energy due to absence of kinetic energy.  The variational principle that is the main tool in Lagrangian mechanics can indeed derive the equations in the Newtonian mechanics. One of related research topics using variational principle to fluid mechanics is a phase boundary propagation, which can be analyzed by the phase field model or diffusion-interface model (Cassel, 2013).  This model explains the diffuse phase boundary without surface tension appeared in Newtonian interfacial physics between liquid and gas. According to the second law of thermodynamics, the free energy of the system must decrease monotonically to ensure a non-negative entropy production (Singer-Loginova and Singer, 2008). This needs that the time rate of change of phase boundary be expressed by the functional derivative of the free energy functional, which corresponds to the talik total energy flux in permafrost thaw problem.  This study directly and analytically solves the Euler-Lagrange equation based on the stationary action principle rather than the entropy functional used in the phase field method.

• L 24 an extremum of the functional -> a minimum of the energy associated with the functional description of the phase boundary (for clarity)

Thank you for the suggestion.  We would revise it, accordingly.  Equation (15) uses the method of Lagrange multipliers which is a common tool in the machine learning field (e.g. maximum entropy principle) for optimization.  We plan to indicate the name of the method (method of Lagrange multiplier) for readers to understand the physical interpretation.

• L 32 stabilizes thermokarst lakes -> stabilizes the size? shape? growth rate? be specific

It implies slowing the “growth rate”.  We will be more specific.

• L 45 … above an unfrozen water body (CITATION NEEDED)

We will provide an appropriate citation here.

• L 50 … bed consequently subsides (CITATION NEEDED)

Will be provided.

• L 84 advanced the technique by including (parameter of)<-remove advective heat transport

We will revise the part accordingly.

• L 141 k_L needs units!

Thank you.  It should be W/(m·֯C).

• L 151 should read: the letters in bold denote vectors.

We will revise the part accordingly.

• L 152 remove ‘also’

We will revise the part accordingly.

• L 181 in the horizontal direction -> the horizontal gradient

We will revise the part accordingly.

• L 182 this is not easily interpreted, especially for the general audience of this journal. Are you implying Stoke’s theorem? I think it would strengthen the derivation to begin from a more ‘certain’ or understandable place than eq 183

No, it is not Stoke’s theorem although it is related.  We quantify the free energy functional for the Euler-Lagrange Equation.  Please see the response above.

• Eq 13 I think it may make more logical sense to present this in the opposite direction - the integral along the phase boundary (line) is not something that I can interpret easily or can be visualized, whereas something more like the flux across the phase change surface, or the volume integral of the total energy in the lake is more easily interpreted. I would start with the resulting equation and state which theorem (Stoke’s?) is used to get the first equation. Importantly providing a physical interpretation (in more simple terms) of what each expression (start and derived result) means and how it is useful and what it tells us about the system. This would greatly increase the utility of the work for those who are less interested in the mathematics and more interested in their application.

We will add the paragraph presented above at the beginning of Theory chapter to clarify the approach before the derivation.  Also, Equation (15) uses the method of Lagrange multipliers which is a common approach in machine learning field, lately.  We hope the name of the method helps for readers to understand the physical interpretation.

• l 187 section heading is redundant: optimum phase boundary shape - is sufficient

That is right. It will be revised, accordingly.

• l 190 “here, we present the …” - present is not a good word choice

We will adjust it.

• l 193 “weighted phase boundary area” this is not easily visualized/interpreted. Can you describe more concretely what this key phase means and what quantity is weighted along the phase boundary before using it (it occurs several times throughout)

The weight is alpha. The area element without weight in Cartesian coordinate can be expressed as,

(equation here, see supplement).

The corresponding phase boundary area is,

(equation here, see supplement).

However, without the weights the optimum shape is sphere, which is irrelevant here.  As we consider the vertical thermal gradient, the phase boundary area must be weighted as

(equation here, see supplement).

The purpose of Section 2.1 is to show validity of these weights.

• l 195 it does not seem logically evident that the shape of a talik would preferentially minimize the total permafrost thaw given an amount of incoming energy. Please either provide a reference, explain the logic, or reconsider. The energy would simply Óow according to the thermal gradients, and more energy would be used where gradients are strongest, no? If there is data for the temperature gradients here or in other lakes that would be helpful.

The free energy of the system must decrease monotonically to ensure a non-negative entropy production (the second law of thermodynamics).  This law suffices for the minimization of talik total thaw.  Also, it is obvious from the energy conservation equation at the phase boundary (Equation 1) that the milder thermal gradient in solid permafrost results in faster thaw.  We prepare a diagram below for convenience.

(image here)

• eq 14 the volume is a double integral and not a triple integral - why?

It is because phi(x,y) describes the depth of talik (distance from the surface).

• l 224 (and future occurrences) “this 3D Stefan equation” should it not be the solution to the 3D Stefan equation? Where do you introduce the Stefan equation? eq 4? State in the text when this occurs otherwise this conclusion to the methodology here doesn’t really make sense

Stefan equation describes the phase boundary depth (active layer depth or frost depth) under the quasi-steady state heat balance.  The solution of the Euler-Lagrange equation describes the phase boundary depth beneath a seasonal water body.  Hence, we think the elliptic function is the 3D Stefan equation (not a solution of Stefan equation).

• eq 24 - define units for new parameters, and define r_{deg} in text. Where is this equation from? Appendix? if so please say where to find derivations, if not consider including this.

It is dimensionless (ratio; L/L). We could not find such a degradation ratio as a function of thawed permafrost depth in the literature. Therefore, we newly introduced it.  Equation 24 is the definition.

• l 245 unlikely -> not a

There is no evidence that it is impossible either.

• l 239-248 note if there was evidence of segregated ground ice in these sites

Please see the references, Kanevskiy et al., 2013, 2016.

• eq 25 r_{sub} what is this?? Is it the same as r_{deg}?? rename or deÒne

Sorry, that is a typography.  All r_sub should be r_deg.  We will fix it in the next revision phase.

• l 269 missing work : … not eh western shore determined form….

We will revise the part accordingly.

• lSSS 3.1 describe thaw rates a bit more (especially in reference to Figure) the range is wide, I assume the high thaw rates are observed in a similar location… also how complex is the shoreline shape? Does it vary around the lake?

Relatively flat lake bathymetry observed by Lenz et al. (2016) suggests uniform subsidence with deep talik.  However, we understand that it is complex near shoreline.  This result actually led to the discussion chapter.

• l 283 solution to … also remove ‘the’ based on 27 talk thickness …”

We will revise the part accordingly.

• l 301 model geometry not geometry model

We will revise the part accordingly.

• l 304 lake expansion is most rapid

We will revise the part accordingly.

• l 324 “wind erosion effect” not mentioned until now, please elaborate earlier or save this for the section reporting it exclusively

We will introduce the wind erosion effect earlier.  Thank you.

• l 341 is the assumption that the radial thermal gradient is zero accurate? Other publications report much more rapid than vertical thaw (though my focus is discontinuous PF) see McClymont et al. Devoie et al. work at Scotty Creek. Please cite something or report thermal gradients to support this.

Thank you.  We will.

• l 371 the preceding discussion all hinges on the zero lateral gradient assumption - please highlight this otherwise it seems unlikely

Thank you.  We will make this clear.

• l 395-396 this is the Òrst mention of anthropogenic processes, and none of the bullet points are direct human activities. Suggest to remove this idea unless there is an additional section on anthropogenic change

We will revise the part accordingly.

• l 408 what are horizontally oriented lakes? Please describe as this term is not clear

The horizontally orientation here is direction dependent elongation.  We will define the term.

• l 443 easier to understand would be: Wind-driven wave action make the water bodies round…. (remove asymptotically - this does not belong)

We will revise the part accordingly.

• l 493 why is it more rapid?

We will remove “more rapid talik”.

• l 510 for the horizontal stage (COMMA), A in Figure 6,

We will revise the part accordingly.

• l 527 solution to 3D Stefan equation is limited

It should be kept as is.  This is not a solution of 3D Stefan equation but the 3D Stefan equation.

• l 535 what about anisotropic thermal properties? Maybe also discuss these as well?

Sorry. We cannot catch it.  Thermal properties (e.g. thermal conductivity, latent heat, and thawing temperature) should be isotropic, constant, and uniform.  However, the thermal field is anisotropic (e.g. vertical temperature gradient).

• l 541 state the thermal effect on thermokarst morphology, I am not sure what this refers to in the manuscript as is written now, so it is either unclear or unsupported

Thermal effect is described by the 3D Stefan equation.

• l 557 wieighted phase boundary (again weighted according to what?)

It was explained about it earlier.  Please refer to it.  Thank you.

• Suggestion: due to the lack of evidence supporting the conclusions on shape, it seems that the argument is stronger for the phases of formation and evolution of thermokarst lakes, so it may be more relevant to report more strongly on this aspect in the abstract and conclusion? As an alternative/addition to my previous comments on data comparison

We plan to adjust the title and the manuscript to reflect the discussion section.

Reviewer 2

• The discussion paper analyzes geometry of taliks forming below permafrost thaw lakes, compares the resulting model to field data from one site, and provides a more qualitative discussion on permafrost thaw lake dynamics and influence of those dynamics on talik growth. The work provides insights into environmental controls on thaw lake taliks and implicitly on their evolution in a warming Arctic. The subject is of significant interest to the readers of TC and the approach used is relatively novel and appropriate. Overall, the discussion paper provides some interesting insights and is a welcome contribution. I do think the manuscript could be improved by paying some more attention to the presentation of the model in Section 2 and to implications of the work.

Thank you for the overall positive evaluation.

• The title of the discussion paper focuses on the analysis of Section 2, but doesn’t provide an adequate indication of Section 4, which provides a much richer picture of thaw lake and thaw lake talik dynamics. This paper might be more impactful if the authors can find a title that better reflects the breadth of the analyses.

It is a good idea as this title was created at the earlier mathematical derivation stage.  We suggest, “Evolution of Thermokarst Lake and Talik Geometry and the New Stefan Equation based on Lagrangian Mechanics” 

• Description around the derivation in Section 2 is often unclear and in some places imprecise and will make it difficult to follow for some. For example, it’s clear in the scalar equation 1 that q_f is the energy available to thaw permafrost per unit time and area. Generalizing to the vector equation 5, it’s a little easier to understand q_f as the velocity of the moving phase boundary scaled by the volumetric latent heat of fusion for water ice to make it an equivalent heat flux, as in equation 2. Referring to q_f as the ‘fusion heat vector’ is a bit obscure. This might be easier to follow and would avoid that jargon by doing the analyses in the velocity v instead (eg. By dividing both sides of eq 5 by \phi \rho L see Eq 2). An alternative approach would be to clearly describe the physical interpretation of the vector q_f around Eq. 5, give q_f a better name and stick to that name in the rest of the manuscript.

Thank you for the suggestion. We will try to improve the first introduction of heat equation. 

We plan to change the name to q_f from “heat for fusion or thawing” to simply “thaw energy” while this term can be refreezing with negative q_f.   Also, the vector may be named, ‘talik expansion vector”. 

• It would be helpful to summarize assumptions behind Eq 1 when Eq 1 is introduced. This is addressed somewhat in Section 4, but it would be helpful have that stated more explicitly.

The assumptions for Equation (1) will be included in the revised manuscript.  They are as follows:

The material of permafrost and talik is fully saturated with ice and water, respectively.  The thermal constants (thermal conductivity, latent heat, and thawing temperature) are constant and isotropic. Change in volume of water on thawing and freezing is negligible.

• The reader needs to know why the functional F is introduced this way in Equation 15 (i.e. you want to minimize the boundary area for a specified thaw volume, the symbol \lambda is a Lagrange multiplier, etc.)

As explained above for Reviewer 1, the second law of thermodynamics ensures a non-negative entropy production throughout the system. We will include this explanation as proposed above. 

• Also in line 197, it would be clearer to say “for a specified talik volume” instead “for the total talik expansion”. Similarly, the sentence starting on Line 194 could be clarified.

We will adjust the part on line 197 accordingly.  Thank you for the good suggestion.

• This paper contains several insights that could inform representations of thaw lake dynamics in Earth System Models. If possible, it would be useful if the authors could comment on that.

The analytical solution presented in this article can be potentially incorporated in the global or regional scale Earth system model to describe missing sub-grid scale processes such as lake dynamics with minimal additional computational resources.  We will improve the manuscript in this regard.

References

Cahn, J. W., & Hilliard, J. E. (1958). Free energy of a nonuniform system. I. Interfacial free energy. The Journal of chemical physics, 28(2), 258-267.

Cassel, K. W. (2013). Variational methods with applications in science and engineering. Cambridge University Press.

Singer-Loginova, I., & Singer, H. M. (2008). The phase field technique for modeling multiphase materials. Reports on progress in physics, 71(10), 106501.