The paper is focused on evaluating bulk algorithms by comparison of the computed heat flux  against the observed fluxes. The latter are obtained using the eddy covariance method over an ice-covered lake. The study is using a rather unique four-year dataset. Such an analysis has been rarely performed for rather small ice-covered lakes which represent complex conditions. Namely, the proximity of a coastline and an edge of a boreal forest can break the limits of applicability of the Monin-Obukhov similarity theory. The paper clearly shows that heat fluxes, especially, the sensible heat flux is underestimated by the standard bulk algorithms which might point to the fact that some physical mechanisms are not taken into account by the considered bulk algorithms. I find that the study certainly advances our knowledge of turbulent exchange processes over various types of surfaces. It also highlights the limitations of the existing models.

The main drawback of the study is the lack of a detailed analysis and a rather crude approach to using the existing models. At least, a better presentation of the possible sources of errors is needed. This concerns the following issues:

• the effect of stability on exchange coefficients and computed fluxes is studied by prescribing constant exchange coefficients and comparing the performance of such a model with other models where the stability-dependent transfer coefficients are used. However, all the models use different roughness lengths, i.e. the neutral transfer coefficients. (The authors should give the exact values of the neutral transfer coefficients of the algorithms, if possible) Thus, one cannot fully separate the effect of stability by comparing their results. Unfortunately, the used observations are limited to one level only and thus it is not possible to fully get rid of the uncertainty associated with the unknown roughness lengths for momentum and heat. The way to proceed would be to prescribe various values of roughness length, or find the best fit based on their data and evaluate the roughness length models.

• The authors use the value 1.8x10-3 for the neutral transfer coefficient at 1.7 m height in their most simple bulk model. Is this value representative for a smooth lake ice? What value would result from their own dataset? Why don’t the authors try to estimate the roughness length (or the neutral transfer coefficient) for momentum and heat/moisture based on their EC data? Moreover, roughness length might be dependent on wind direction. The authors do not compare the observed momentum fluxes with those obtained using bulk approach. The question is why? Such an analysis would be helpful in identifying the uncertainties related with roughness length and also stability functions. From their dataset and using best-fit roughness length it is possible to estimate the integral stability correction functions for different z/L (Psi-functions) and compare them with those prescribed in the bulk models. Such an analysis for the momentum Psi-functions would be free from the uncertainty associated with the uncertainty in the surface temperature. Also,  the presented results suggest that increasing roughness length for heat might improve the bulk model performance with respect to sensible heat flux, at least in March-April. The authors should include the sensitivity study to prescribing various roughness lengths.

• The EC fluxes are associated with a certain footprint area. What is the area of such a footprint for the considered site, what types of surface are expected to affect turbulence over the measurement site? Does coastline and forest frequently occur in the footprint area?

• Why do authors use a bulk model with the roughness length prescribed using Charnock formula for open water? What is their rationale behind that? Of course, for a snow-covered surface a Charnock-like formula for z0 was suggested by Andreas, but can it be applied to the considered lake? The authors do not study this issue.

Minor comments

Line 73 – all the models are to some extent semi-empirical. Better to avoid such a sentence without explicit explanation what is meant.

Line 119 – what is the sensitivity of the output of bulk algorithms to the used value of surface emissivity varied in the range of natural variability? For example, how much would the fluxes change if eps = 0.98 is used?

Line 200 T* is not a dimensionless temperature, but it is the temperature scale which has dimension [K]

Lines 230-235 it should be better explained how the increase of solar radiation results in negative heat flux and T0-Ta at daytime in spring. Obviously, this is only possible if somewhere around the lake the daytime heat flux and T0-Ta become positive.

The authors consider turbulence measurements over a lake in Finland. They compare different bulk paranetrizations of the fluxes with Eddy Covariance measurements, which had been obtained during four winter episodes. They find small differences between all parametrization models and conclude that the best agreement with the EC measurements is found for a method using a constant transfer coefficient without any stability correction.

The paper is well written and clearly organized. However, I think some improvement is necessary before its publication.

Major Revisions

1) Equation (14) is not correct. A logarithm is missing in the first term of the r.h.s. Please check your codes if they are correct (with logarithm).

2) The main conclusion (see abstract of the manuscript) is that the assumption of a constant transfer coefficient being independent on stratification is the best one. This needs much better explanation. Modellers might get the idea to ignore the stability correction in their runs in general. But this is against all previous experience over decades from observations, theory, and Large Eddy Simulation. So, if this is really the result, then the reader must be better convinced that it is not an artefact. Possible reasons might be conditions violating Monin Obukhov similarity or problems with the accuracy of measurements (e.g. due to influences of the boxes neat the small ‚tower‘?) and many others.

To better convince the reader I find it necessary to show results (fluxes) obtained by the EC method as a function of z/L or of the Richardson number (as in Grachev et al., 2007) or in many other papers (e.g. most recently Srivastava, Gryanik et al.). It would be helpful to show results of the phi- function (eqs. (9,10) and that behind the SHEBA equation (14) (see Grachev et al., 2007) as a function of z/L.

3) It is difficult to interprete the differences between all three schemes based only on the scatter plots (Figure 5).

4) When the final result remains unchanged, it needs to be explained more careful. It should be written that further research is necessary to test the robustness of this result. E.g. measurements over other lakes are necessary before a general suggestion to modellers can be given. Such results could depend on the lake size, where the flow over small lakes might be more inhomogeneous than the flow over large lakes and inhomogeneity might hide the stability dependence.

5) The underestimation of fluxes might be due to errors in the roughness lengths (especially the ratio between roughness length for momentum and for scalars is uncertain).

Minor Revisions

1) The introduction is interesting but the description of main goals and their explantion are rather short.

2) Line 66 with ‚the‘ bulk…

3) line 78: I would not say that the EC measurements really overcome heterogeneity problems.

4) line 82: present ‚a‘ unique

5) Considering figure 1 and description, I was not sure if data were selected with wind directed only along the fjord to minimize inhomogeneity impact.

6) equation (9): insert space between zeta and ‚for‘

7) Sometimes, the article is missing, e.g. line 247 (The surface energey …)

8) Line 247: better give equation for the energy balance

9) Figure caption 3. Sign of flux vice versa as compared with line 280 ?

10) Line 13: is the underestimation a problem due to the phi-function or due to the roughness assumption?

11) Line 295 and similar …. better: results over trhe melting surface

12) Line 310: patterns

13) paragraph starting at line 388: This comes very late although it is the main result (see abstract and above).

References

Srivastava et al. (2020) On stability correction functions over the Indian region under stable conditions. Meteor. Appl., 27, e1880

Gryanik, et al. (2020) New modified and extended stability functions for the stable boundary layer based on SHEBA and parametrizations of bulk transfer coefficients for climate models. Journal of the Atmospheric Sciences, 77, 2687–2716.