I greatly enjoyed reading this paper. The method for efficiently coupling the nonlinear surface energy balance to the linear subsurface heat conduction is a clever piece of matrix algebra, but it is not just that; it directly relates to a point of contention in the lively interactive discussions of Brun et al. (2022) and Potocki et al. (2022) concerning the mass balance of the Everest South Col Glacier. There are limitations, however. Many processes are generally handled sequentially in snow models (Clark et al. 2015), but this paper only couples two of them. Only idealized test cases are shown and not full model performance in real applications. From the test case results, I could take contrary conclusions that the added complexity of coupling is not needed and the standard skin-layer formulation is fine as long the timestep is not made too large, which is well known (“not too large” could still be prohibitively small for thin layers, though).

Specific comments:

Author list

“Brun Fanny” might like to have her name turned around.

Introduction

I don’t recommend writing a comprehensive review of SEB formulations, but only giving recent examples of applications of a skin layer and no examples using a finite surface layer in the introduction, rather than original model development papers, gives a distorted view. An uncoupled skin layer has been in use for snow models at least as far back as Yamazaki and Kondo (1990). There is a snow surface layer temperature in Anderson (1968).

22

There are many “numerical models” that are not snowpack/glacier models.

32

The surface energy balance is described as “profoundly non-linear”. Actually, this is a pretty benign nonlinearity in the field of nonlinear equations; it does not have multiple or chaotic solutions.

49

The “infinitely small horizontal layer” would be better described as infinitely thin.

65

While Eq. (1) is more generally applicable, it could already be emphasized that this is invariably implemented as a 1D model with T a function of z.

79

I think that there will be very few exceptions to this “usually” of allowing snow temperature to exceed the fusion point before calculating melt, but there are examples of models with phase changes over a temperature range in Albert (1983) and Dutra et al. (2010).

136

 “SNTHERM (Jordan, 1991), Crocus (Vionnet et al., 2012)”

146

Another step is required if the calculated melt exceeds the available snow mass.

234

LWout and H are given as examples of fluxes that are nonlinear in the surface temperature; L should also be mentioned as intrinsically nonlinear. H as defined by Eq. (B1) is only nonlinear if C_H is a function of surface temperature. It is, through Ri_b here, but models often neglect this nonlinearity because of the complexity of the resulting derivatives; it is not clear if that is done here.  A supplement giving the elements of the Jacobian might be a useful addition.

261

I understand the problem, but I don’t understand the benefit of returning the solution to the vicinity of the discontinuity.

The SEB should have a unique solution, but the Newton method is not guaranteed to find it. It can get trapped in a cycle of states around the solution. This situation can be diagnosed from the SEB, but I think that most models just give up and select the last iteration. Does the modified Newton method avoid this problem?

265

Another solution in use, with its own numerical errors, is to linearize the SEB and solve it in one step without iteration (e.g. Best et al. 2011). This is essentially the Penman-Monteith method.

Equation (11) and following

Be consistent in making diag, up and low superscripts or subscripts.

284

“The above equation” is Eq. (13).

286

“invert A_diag”

322

“cells which then become”

379

No refreeze in this test case.

421

“Concerning the glacier test-case, Fig. 3 shows”

424

“by about 0.50 K”

440

I have, indeed, seen time step oscillations like this in class 2 simulations. They are not the same as the well-known and catastrophic instability of the explicit Euler method with too large a timestep. Considering the wide use of class 2  models, a stability analysis to understand the origin of these oscillations (not necessarily for this paper) might be of interest. 

460

“only marginally worse”

490

Divergence of the class 2 model from the reference as the mesh is refined in the glacier test case (Fig. 10) is an odd result. I guess that this could happen if the time step in these mesh refinement tests is larger than in the reference. If so, this needs to be stated in the text.

Having said that, it is not apparent that the 225 cell simulation is worse than the one with 45 cells in Fig. 10c.

504

“the backward Backward Euler method” sounds like it goes forward. Just one “backward” required.

509

“mesh size too small”

6.4

Having found G from the SEB, the obvious thing to do in a class 2 model is to use it as a flux boundary condition for the internal temperature calculations. Can any real class 2 model be found that uses the surface temperature as a Dirichlet boundary condition? If not, section 6.4, Fig. 14 and the last sentence of the conclusion should be deleted. A note that this would be the wrong thing to do will suffice.

Search the text for “, that”. In all but one case, it should be “that” or “, which”.

Albert (1983): https://apps.dtic.mil/sti/pdfs/ADA134893.pdf

Anderson (1968): https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/WR004i001p00019

Best et al. (2011): https://gmd.copernicus.org/articles/4/677/2011/

Clark et al. (2015): https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2015WR017200

Dutra et al. (2010): https://doi.org/10.1175/2010JHM1249.1

Yamazaki and Kondo (1990): https://doi.org/10.1175/1520-0450(1990)029<0375:APMFSS>2.0.CO;2

Summary: 

This work proposes a methodological improvement to surface energy balance modeling over frozen ice surfaces by merging the benefits of two diverging current approaches to coupling air temperature and ice temperature.  The coupling approach appears effective and insightful and is an important contribution to the field. The paper presents two case studies, one over snow and one over a glacier (with highly idealized implementations) as demonstrations of the accuracy. There is a well motivated exploration of the implementation's dependence on time and spatial resolution, which are not only practically important for anyone wishing to implement this method, but also provide the opportunity to discuss numerical stability. 

General comments:

It is exciting to see this paper address both snowpack and glacier surface energy balance. It would be good to discuss (briefly) the physical similarities and differences (structure, air content) between the two.

There is a consistent overuse of commas in the setup ‘,that’ (many of which should be ‘which’ with no comma)

The manuscript is clearly structured in introducing a new method to approach temperature and melt numerical modeling and then applying that method to two test cases. However, the test cases are very specific and thus convey limited information about the broader application of the method – these limitations should be discussed, especially as a future goal would likely be to apply this numerical routine to more complicated cases.

Lastly, the finding that a coupled surface model can outperform other models at coarser grid sizes is implied here to be more computationally efficient due to the change in mesh size. However, this is not generally true when you are also changing the numerical scheme, so the assertion of computational cost savings which maintaining accuracy (as claimed here) should be backed up by either reports of the time taken for the computations and/or a clear statement that the numerical implementations are computationally identical by construction. This, if true, should also be mentioned in the conclusion as it is an important outcome! This is somewhat related to the discussion of numerical reduction (back to the same order of the original models) that you get from the Schur complement, but they are not discussed together and the data are not shown.

Specific comments:

L3-4: “This surface energy budget is the sum of the various surface energy fluxes, that depend on the input meteorological variables and surface temperature, and to which heat conduction towards the interior of the snow/ice and potential melting need to be added.”  the comma between ‘fluxes’ and ‘that’ is incorrect, as are similarly positioned commas throughout, and ‘that’ should be ‘which.’ 

L2-4: ‘and to which heat conduction towards the interior…” this sentence is unclear to me

L26: once the SEB acronym is introduced, it should be used consistently in the paper

L25-30:  There is a large focus on the nonlinearity of SEB processes, which is important but not hugely challenging in the modeling field, as many of the nonlinearities are easily solved. It would be good to mention this and discuss sources of nonlinearity in a more mechanistic sense. For example, the “regime change” mentioned is due to thermal energy being used for processes with different reaction coefficients in warming frozen ice vs. phase change. This will help build intuition to support the truncation method discussed later. Perhaps mention another example. 

L42: which domain? The ice domain?

L63: specify Fourier’s law of heat conduction

L90-95: specify the sign convention used for fluxes

Figure :1: clarify the meaning of the blue/orange colors of dots in the figures. Additional labels within the diagram would improve the clarity of the figure. It is also somewhat redundant to label Class 1 as a), class 2 as b) etc. since they are all in essentially the same panel. Consider just labeling the columns as class 1, class2, this paper.

L115: “We therefore do not treat the finite elements method, which is for instance used in the SNOPACK model.” -> “We therefore exclude implementations of the the finite elements method, such as in the SNOWPACK model.”

L179: omit the word “let’s” or use “we” instead

L275: the introduction of new terminology in representing equations 5 and 9 as a block-matrix system with decomposed components (Adiag, etc.) requires additional explanatio of the correspondence between terms in Equations 5 and 9 and their placement in Equation 11.

L350: the bottom boundary no heat flux assumption seems strong to me, or at least appropriate in a limited set of conditions. A citation or further discussion of this would help.

L353-355: these constants are also introduced on L 66 and 72-75, use the same symbols here to connect them. 

L354: thermal conductivity of ice is temperature sensitive! If making this assumption, please explain why it is warranted in this case (i.e., the temperature ranges reasonably experienced in this case are small enough that there is not meaningful variation?)

L385: albedo over what wavelength range? In most of the visible spectrum, this would be a quite low value in clean snow. Further, the longwave emissivity of snow is more density dependent due to the presence of air. It seems reasonable to use 1 for this approach, as the emissivity is still usually quite high

L435: the “lag” of one time step mentioned here is interesting and well explained. The impacts of this on interpreting a snow model output may be sensitive to the time step. If there is a long time step, this would be problematic as it may prevent modeling melt. A short time step may be more resilient to this lag.

L440: the observation that numerical instability is leading to differences between class 2 models and other models is interesting and seen clearly in Figure 4. The fact that this is not happening in the glacier model is only vaguely referenced. A direct comparison of the reasons for this – if there is an inherent numerical instability in class 2 models, why is there not an instability in the glacier model? Is all of the oscillation occurring in the meltwater percolation?

Figure 4: it is impossible to see the coupled surface line in panel b – consider adding markers or some other formatting choice which would allow us to see it clearly. Layering the coupled surface model on the front may help if markers are not working favorably.

Figure 6: right panel y axis would benefit from additional labels

L490-495: as worded, the phrase “the class 2 model exhibits the largest phase change rate errors for an initial number of cells of 225”  is ambiguous – is 225 the worst number of cells for C2 models or is C2 the worst option when working with 225 cells? From the graph, it is the second option, which is potentially less important than discussing the fact that for the other two model options, a larger number of cells generally confers better performance (within the parameter space explored here), but that is not the case for class 2 when moving from 90 to 225. Why might this be?

L504: implicit backward Euler method?

L506: “explicit” ?

L509: “too”

519: is the Dirichlet approach actually used? If not, it is not relevant to compare it here

Figure 13: this is a time series of temperature, not a graph of numerical instabilities and should be labeled as such. It seems that the goal is to point out the larger variance in the higher number of cells-driven runs, so I would recommend either adding the time average standard deviation, or converting this plot to a time series of deviation from some sort of rolling mean in order to focus more on the instability-driven variance. Or, add a second column that contains histograms of that variance for each case.

L560: the level of accuracy is similar but not identical

L613 “This approach is, for instance, used in the Crocus model” add commas

The authors present an approach to numerical modeling of snowpack or glacier interface with atmosphere using a finite volume method discretization of thermodynamic relations. The novelty of the approach lies in coupled computation of heat transfer through the ice/snow and the thermodynamic balance at the surface. The authors provide sufficient numerical experiments to support the agreement of their implementation with previously published results.

The only critical comment I would like to make is the relatively vague mathematical description of their approach, or the problem at hand. The authors discus the Fourier's law for the heat transfer in ice (Equation 1) and the balance of energy fluxes at the ice surface (Equation 3).  Then, they immediately follow on to numerical discretization, leaving the reader curious as to what assumptions and specific method choices they made. I would outline below a few of my concerns.

The authors start with the heat equation:

∂_th - div (λ grad(T)) = Q

where

h = c_p(T-T₀) + ρ_wLθ

which leads to 

c_p∂_tT + ρ_wL∂_tθ - div (λ grad(T)) = Q.             (1)

In the subsequent paragraph they discuss issues with representing the effects of phase changes on the temperature, but I believe they mean that they neglect the ρ_wL∂_tθ term in their model. Please state that clearly.

Moving on, the authors jump to Equation 5, where they present the discretized version of (1) using finite volumes. It would be useful to state the implicit assumptions here, that the three dimensional equation (1) is now considered as one-dimensional equation

c_p∂_tT  - ∂_z (λ ∂_zT) = Q,

which is then integrated over each "volume", which in this case is segment of length Δz_k. This integration, along with replacing the point variables with their volume averages (with abuse of notation: T_k = 1/Δz_k ∫Tdz), and using fundamental theorem of calculus (we are in one dimension now, no need for divergence theorem) gives

Δz_kc_p∂_tT  -  (λ ∂_zT)_{k+1/2 +} (λ ∂_zT)_k-1/2 = Δz_kQ, where subscripts k+1/2 and k-1/2 refer to the (top and bottom) endpoints of the cell Δz_k

At this point the authors introduce the approximation of the (λ ∂_zT)_k+1/2 term with Equation 6. I am curious, however, whether it is not better to leave the term (λ ∂_zT)_z=surf at the top of the first layer as is, and replace it with the term G from the surface energy balance equation (3)? I am not sure whether this is the way the authors achieve coupling, or whether they still discretize the temperature gradient at the ice surface using the surface temperature and half of the top layer size? The authors discuss in lines 103-105 that term G depends on surface temperature and temperature within ice, which indicates that this term is indeed discretized. 

I would urge the authors to provide a more detailed and careful mathematical description of their work, as it would improve the reproducibility of their result, not only for the finite volume method community, but also researchers working with other numerical discretizations.